Advanced multi-perspective computer simulation as a tool for reliable consequence analysis

Process Safety and Environmental Protection 9 0 ( 2 0 1 2 ) 129–140

Contents lists available at ScienceDirect

Process Safety and Environmental Protection journal homepage: www.elsevier.com/locate/psep

Advanced multi-perspective computer simulation as a tool for reliable consequence analysis Spyros Sklavounos a,b,∗ , Fotis Rigas a a b

National Technical University of Athens, School of Chemical Engineering, 15700 Athens, Greece Hellenic Army Academy, Department of Mathematics and Engineering Sciences, 16673 Vari, Greece

a b s t r a c t Major accidents involving hazardous materials are a crucial issue for the chemical and process industries. Many accidental events taken place in the past showed that dangerous substances may pose a severe threat for people and property. Aiming at loss prevention, a series of actions have been instituted through international regulations concerning hazardous installations safety preparedness. These actions involve efficient land-use planning, safety studies execution, as well as emergency response planning drawing up. A key factor for the substantial consideration of the above is the effective prediction of possible accident forms and their consequences, for the estimation of which, a number of empirical models have been developed so far. However, (semi-)empirical models present certain deficiencies and obey to certain assumptions, thus leading to results of reduced accuracy. Another approach that could be used for this purpose and it is discussed in this work, is the utilization of advanced computational fluid dynamics (CFD) techniques in certain accident forms modeling. In particular, composite CFD-based models were developed for the simulation of several characteristic accident forms involving isothermal and non-isothermal heavy gas dispersion, confined and unconfined explosion in environment of complex geometry, as well as flammable cloud fire. The simulation cases were referred to real-scale trials allowing us to conclude about the validity of the quantitative results. Comparisons of the computational predictions with the experimental observations showed that obtained results were in good agreement with the experimental ones, whereas the evaluation of statistical performance measures proved the simulations to be statistically valid. © 2011 The Institution of Chemical Engineers. Published by Elsevier B.V. All rights reserved. Keywords: Loss prevention; Consequence analysis; Environmental impact; Gas dispersion; Explosion; Cloud fire

1.

Introduction

Huge amounts of chemicals are produced and used daily, while new chemicals are constantly being introduced in the world market with an average rate of 400 new substances to be registered every day in the Chemical Abstracts Service (CAS REGISTRY) (Binetti et al., 2008). Many of them are classified as hazardous materials, due to the fact that they present toxic, flammable and/or explosive properties, able to be involved in a major accident (toxic dispersion, fire and/or explosion), usually severely impairing people and property. The tragic chemical accident that happened in Seveso in 1976 was the cause for the Council of European Communi-

ties to establish in 1982 the 82/501 Directive aiming at the prevention of major accidents involving dangerous substances. In 1996, the Directive was revised and extended introducing new requirements relating to safety management systems, emergency and land-use planning. As a result, the Council Directive 96/82/EC on the control of major-accident hazards (the socalled Seveso II Directive) was adopted. Recently, the Directive 2003/105/EC of the European Parliament and of the Council of 16 December 2003 amended the Council Directive 96/82/EC on the control of major-accidents. In addition to accident prevention, Seveso II and its recent amendment aim at the limitation of the consequences of such accidents, not only to people (safety and health aspects), but also to the environment (environmental aspect) (EEC, 2003, 1996).

∗ Corresponding author at: National Technical University of Athens, School of Chemical Engineering, 15700 Athens, Greece. Tel.: +30 2102824365; fax: +30 2102816147. E-mail addresses: [email protected], [email protected] (S. Sklavounos). Received 3 November 2010; Received in revised form 13 April 2011; Accepted 15 June 2011 0957-5820/$ – see front matter © 2011 The Institution of Chemical Engineers. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.psep.2011.06.008

130

Process Safety and Environmental Protection 9 0 ( 2 0 1 2 ) 129–140

Nomenclature Cε1

k–ε turbulence model unitless constant equal to 1.44 Cε2 k–ε turbulence model unitless constant equal to 1.92 C k–␧ turbulence model unitless constant equal to 0.09 CH4 .molf molar fraction of methane in air (dimensionless) Di kinematic diffusivity of component i (m2 s−1 ) total black body emission (W m−2 ) Ebv energy release rate (J s−1 ) Er Gv spectral incident radiation (W m−2 ) h specific static (thermodynamic) enthalpy (m2 s−2 ) htot specific total enthalpy (m2 s−2 ) molar concentration of the reactant compo[I] nent I (mol m−3 ) k turbulence kinetic energy per unit mass (m2 s−2 ) or the elementary reaction K˛v absorption coefficient (m−1 ) scattering coefficient (m−1 ) Ksv LFL lower flammable limit (dimensionless) number of components in the reacting mixture Nc (dimensionless) p static (thermodynamic) pressure (kg m−1 s−2 ) modified pressure (kg m−1 s−2 ) p P total pressure (kg m−1 s−2 ) Pk production rate of turbulence (kg m−1 s−3 ) buoyancy production term (kg m−1 s−3 ) Pkb Prt turbulent Prandtl number (dimensionless) spectral radiative heat flux (W m−2 ) qrv Rk the elementary reaction rate of progress for the k reaction SE source term due to the chemical reaction rate involving component I (kg m−3 s−1 ) energy source (kg m−1 s−3 ) Si SI mass source of fluid component i (kg s−1 ) time (s) t u fluctuating velocity component in turbulent flow (m s−1 ) velocity vector (m s−1 ) U Ut velocity tangent to the wall (m s−1 ) static thermodynamic temperature (K) T Tlow atmospheric temperature in a specific trial (K) constant temperature value equal to 5000 K Thigh the extinction temperature of combustion reacText tion (K) uj fluctuating velocity component in x, y, or z direction (m s−1 ) upper flammable limit (v/v) UFL V volume of integrated cell (m3 ) molecular mass of component I (kg) WI X, Y horizontal coordinates (m) xj x, y or z variables mass fraction of the component I of the reacting YI mixture (dimensionless) Z vertical coordinate (m)

Greek letters A coefficient (dimensionless) coefficient (dimensionless) B molecular diffusion coefficient of component I i (kg m−1 s−1 ) the identity matrix or Kronecker delta function ı (dimensionless) ε turbulence eddy dissipation rate (m2 s−3 ) E turbulence dissipation rate (m2 s−3 ) the von Karman constant   dimensionless parameter eff effective viscosity accounting for turbulence,  + t (kg m−1 s−1 )  molecular (dynamic) viscosity (kg m−1 s−1 ) turbulence dynamic viscosity (kg m−1 s−1 ) t stoichiometric coefficient for component I in kI the elementary reaction k P density (kg m−3 ) density (kg m−3 )  i density of component i (kg m−3 ) dispersion parameter (m)  dimensionless turbulence model constant for the k equation equal to 1.0 dimensionless constant in the k–ε turbulence ε model equal to 1.3 Superscripts T transpose of matrix Subscripts A, B, C,. . . components in the reacting mixture I a component of the reacting mixture the elementary reaction k k component i in the fluid mixture i j x, y or z direction P product components in the elementary reaction k

Instead of using (semi-)empirical models, an alterative is to attempt a computational fluid dynamics (CFD) based approach which implies the development of composite numerical models (Sklavounos and Rigas, 2010). In this paper, the development and application of such advanced methods based on numerical modeling techniques is discussed, as far as the computational simulation of certain accident types is regarded. In particular, recent successful applications of CFD models in consequence analysis are presented. The issues addressed concern the effective simulation of isothermal and non-isothermal heavy gas dispersion plus total dose estimation; confined and unconfined explosions and resulting shock wave propagation in geometrically complex environment, in addition to specific impulse estimations; flammable cloud fires in open space and estimation of the thermal load and developing overpressure as well. Moreover, the results of an innovative multi-parametric study associated with the application of a validated numerical model in tunnel safety systems design is presented. On the basis of existing validation studies and the arisen range of application, some conclusions were made about the advantages and disadvantages of numerical modeling in relation to empirical models and actual experimentation. In short, comparisons of the com-

131

Process Safety and Environmental Protection 9 0 ( 2 0 1 2 ) 129–140

putational predictions with the experimental observations showed that the obtained results were in good agreement with the experimental ones being sometimes more accurate than those provided by semi-empirical models, whereas the evaluation of statistical performance measures proved the simulations to be statistically valid.

2.

Theoretical background of CFD modeling

CFD codes are divided into three parts: pre-processor, solver and post processor for problem definition, problem solution and results processing, respectively. The platform of CFX 5.6 and 5.7 versions provided by ANSYS were utilized for the research purposes of the project. The code consists of four sequential parts:

• Outflow boundary, where the fluid (or a multi-component fluid mixture) exits the domain. It is also one-way flow boundary condition and it is defined by the mass flow rate (or the velocity, or the static pressure) of the outlet. • Opening boundary, where the direction of the fluid flow is a priori unknown, for example due to the generation of recirculation eddies. Fluid temperature and relative pressure should be known for the definition of this type of boundary. • Wall boundary, which can be fixed or slipping in relation to the motion of the fluid. Walls are impermeable to fluid flow and allow the permeation of heat into and out of the domain. Once a boundary is denoted as wall, the relative velocity of the fluid on it is set to zero according to the non-slip condition.

2.1. (1) CFX-Build, in which the computational domain geometry is constructed as a unified set of parametric surfaces and it is divided into a number of smaller control volumes (cells) through the mesh generation process. (2) CFX-Physics Processor, where boundary conditions, physical models, domain fluids and solution scheme are defined. (3) CFX-Solver, where the code performs calculations to the final solution using iterative methods, until the predefined statistical error level has been achieved. (4) CFX-Post, in which the elaboration of the results takes place, supported by outstanding graphics capabilities (visualization of the geometry and control volumes, vector plots showing the direction and magnitude of the flow, visualization of the variation of scalar variables with time through the domain). Quantitative calculations are also performed in this stage.

Turbulence modeling

The equations that describe the fundamental transport phenomena of a particular problem are the Reynolds Averaged Navier–Stokes (RANS) equations that express the conservation of mass, momentum and energy: Fundamental transport equations: Mass balance

The characteristic physical process that fluid flow in realconditions presents is the turbulent motion. One of the most prominent turbulence prediction models implemented in many general purpose CFD codes is the k–ε eddy-viscosity model. The family of eddy-viscosity models suggests that turbulence consists of small eddies which are continuously forming and dissipating, and in which the Reynolds stresses are assumed to be proportional to mean velocity gradients. Both the velocity and length scale are solved using separate transport equations, hence they are named two equation models. The k–ε model has proven to be stable and numerically robust having a well established predictive capability. The k–ε model uses the scalable wall-function approach instead of standard wall functions, improving robustness and accuracy when the near-wall mesh is very fine. Eddy-dissipation models such as the Shear Stress Transport (SST) model and Reynolds stress models such as SSG model (Park and Kwon, 2004; Menter, 1994) are also available in the code. The role of the k–ε model is to introduce two new variables into the transport equations: the turbulence kinetic energy (k) and the turbulence eddy dissipation (ε), leading to the following forms: Continuity equation

∂ + ∇(U) = 0 ∂t

∂ + ∇(U) = 0 ∂t

(2.1)

Heat balance ∂(U) T + ∇(U ⊗ U) = ∇[−pı + (∇U + (∇U) )] + SM ∂t

(2.4)

Energy equation (2.2)



∂P ∂(htot ) − + ∇(Uhtot ) = ∇ ∇(T) + ∂t ∂t

  t

Prt



∇(h) + SE

(2.5)

Momentum balance ∂(htot ) ∂p − + ∇(Uhtot ) = ∇(∇(T)) + SE ∂t ∂t

The mean total enthalpy (htot ) is given by: (2.3)

The equations relating to the flow (fundamental and submodels) may be numerically defined by the specification of certain conditions on the external boundaries of the domain. The following boundary conditions may be imposed on a boundary in the CFX code:

htot = h +

(2.6)

where the contribution from the turbulent kinetic energy (k) is provided by: k=

• Inflow boundary, where the fluid (or a multi-component fluid mixture) enters the domain of interest. It is strictly one-way flow boundary condition and is defined by the temperature and the mass flow rate (or the velocity, or the static pressure) of the inlet.

1 2 (U ) + k 2

1 2 (u ) 2

(2.7)

Momentum equation ∂U T + ∇(U ⊗ U) − ∇(eff ∇U) = ∇p + ∇(eff ∇U) + B ∂t

(2.8)

132

Process Safety and Environmental Protection 9 0 ( 2 0 1 2 ) 129–140

where (B) is the sum of body forces and (p ) is the modified pressure given by p = p +

2 (k) 3

(2.9)

and eff =  + t

(2.10)

The k–ε model assumes that the turbulence viscosity is linked to the turbulence kinetic energy and dissipation via the relation: t =

C K2 ε

(2.11)

The values of k and ε are directly calculated from the differential transport equations for the turbulence kinetic energy and turbulence dissipation rate: ∂(k) + ∇(Uk) = ∇ ∂t ∂(ε) + ∇(Uε) = ∇ ∂t

  +   r

k

  +   r

ε

2.2.

Combustion modeling

In accident scenarios that include fire (i.e. flash fire, fireball), the need for combustion reaction and radiation transport modeling emerges. The equation of transport for component I with mass fraction YI is described by the following mass balance equation: ∂(uj YI ) ∂(I (∂YI /∂xj )) ∂(YI ) + = + Sl ∂t ∂xj ∂xj

where the source term SI is due to the chemical reaction rate involving component I. Chemical reactions can be described in terms of K elementary reactions involving NC components that can be written as: NC 

 kI I



∇ε +

ε (Cε1 PK − Cε2 ε) k

2 PK = t ∇U(∇U + ∇UT ) − (∇U)(3t ∇U + k) + Pkb 3

(2.12)

(2.13)

(2.15)

The value of the wall shear stress ( ω ) is obtained from:

K 

(



KI

 − KI )Rk

where Rk is the elementary reaction rate of progress for reaction k, which can be calculated using the eddy dissipation model. Several techniques have been developed for use in combustion modeling including the finite rate chemistry and eddy dissipation methods. One of the basic approaches is the eddy dissipation model, which is based on the concept that chemical reaction is fast relative to the transport processes in the flow. This is valid for diffusion flames, transient development of fires and flares (Magnussen, 1989). The model assumes that the reaction rate may be related directly to the time required to mix reactants at the molecular level. In turbulent flows, this mixing is dominated by the eddy properties, and therefore, the rate is proportional to a mixing time by the turbulent kinetic energy k and dissipation ε:

(2.16)

where u∗ ∇y 

(2.17)

ε k

(2.22)

In the eddy dissipation model, the rate of progress of the elementary reaction k is determined by the smaller of the two following expressions:

and u = (C )

1/4

(2.21)

k=1

rate ∝

∗

(2.20)

I=A,B,C,...

where kI is the stoichiometric coefficient for component I in the elementary reaction k. The rate of production/consumption SI for component I can be computed as the sum of the rate of progress for all the elementary reactions, in which component I participates in:

(2.14)

The k–ε model uses the gradient diffusion hypothesis to relate the Reynolds stresses to the mean velocity gradients and the turbulent viscosity. Specifically in the CFX code, it utilizes the scalable wall-function approach to improve robustness and accuracy when the near-wall mesh is very fine. The wall-function approach employs empirical formulas that provide near-wall boundary conditions for the mean flow and turbulence transport equations. Thus, the following logarithmic correlation applies:

y∗ =

 KI I



∇k + PK − ε

SI = WI

ω = uT u∗

NC 

⇔

I=A,B,C,...

where PK is the turbulence production due to viscous and buoyancy forces, which is modeled using:

1 Ut = log(y∗ ) + C uT K

(2.19)

1/2

(k)

(2.18)

In the scalable wall-function approach, the (y*) value used in the logarithmic correlation (2.15) is limited by a lower value of y∗l = max(y∗ , 11.06), where 11.06 represents the intersection between the logarithmic and the linear near-wall profile. As a result, the computed y∗l value is not allowed to fall below this limit and therefore, all mesh points are outside the viscous sublayer. The scalable wall functions allow solution on arbitrarily fine near wall grids, which is a significant improvement over standard wall functions.

Reactants limiter :

Rk = A

ε k

 × min

[I]  VKI

 (2.23)

where [I] is the molar concentration of component I which includes the reactant components.

Products limiter :

 ε   [I]WI  p Rk = AB ×  k

v W P KI I

(2.24)

where P loops over all product components in the elementary reaction k.

Process Safety and Environmental Protection 9 0 ( 2 0 1 2 ) 129–140

2.3.

Radiation modeling

Available techniques for thermal radiation transport computing are the Differential Approximation or P1 model, the Discrete Transfer Radiation Model (DTRM), the Rosseland model and the Monte Carlo model. The approach to thermal radiation transport due to fuel gas combustion has been accomplished through the P1 model, which actually is a simplification of the radiation transport equation (Seaid et al., 2004). The P1 model assumes that the radiation intensity is isotropic in space. The radiative heat flux in the diffusion limit for an emitting, absorbing and linearly scattering medium can be computed as: P1 model :

qrv = −

1 ∇Gv 3(K˛v − Ksv ) − AKsv

(2.25)

The equation for the spectral incident radiation arises by substituting the above terms into the radiation transport equation: −∇ ·



1 ∇Gv 3(K˛v − Ksv ) − AKsv



= Kav (Ebv − Gv )

(2.26)

where A is the linear anisotropy coefficient.

2.4. Discretisation and solution of the governing equations There are several numerical methods used to solve fluid flow equations. The main methods applied to solving a problem that involves partial differential equations are the finite element method (FEM), the finite difference method (FDM) and the finite volume method (FVM). The FVM was originally developed as a special finite difference variation, now developing to the most well-established and thoroughly validated general purpose CFD codes. It is a numerical method for solving partial differential equations calculating the values of the conserved variables averaged across the volume. One advantage of the FVM method over FDM is that it does not require a structured mesh, although a structured mesh can also be used. In FVM, the values of the conserved variables are located within the volume element, and not at nodes or surfaces. In contrast with FEM, FVM utilizes only one node per cell and this is important especially for large meshes and transient calculations, where FEM leads to extreme memory and capacity requirements. Moreover, FVM is especially powerful on coarse non-uniform grids and in calculations that use adaptive grid refinement, where the mesh moves to track interfaces, shocks, or even reacting fronts. Once the composite simulation model (Navier–Stokes equations plus any additional physical sub-models) has been constructed, the entire system of partial differential equations is discretised (namely, it is converted to a system of algebraic equations), which is solved via iterative solution algorithms (Chung, 2002). The solution begins from an initial approximation to the flow solution and gradually iterates to a final result that satisfies the imposed boundary conditions, ensuring that mass, momentum and energy are conserved locally for each grid cell and thus globally over the whole computational domain. The solution process requires the discretisation of the three-dimensional spatial domain into finite control volumes using the mesh constructed in the pre-processing stage. In the FVM numerical technique discussed above, the govern-

133

ing equations ((2.4), (2.5), and (2.8)) are integrated over each control volume, so that the relevant quantity (mass, momentum or energy) is conserved in a discrete sense for each cell. Each control volume surrounds a unique node, in which all solution variables and fluid properties are stored. Afterwards, the integral equations obtained are converted to a system of algebraic equations via the substitution of a variety of finitedifference-type approximations for the terms in the integrated equation representing flow processes such as convection, diffusion and sources. The latter are solved iteratively at nodal points inside each cell aiming at minimization of the RMS (root mean square) residuals (iteration loop error) until the desired convergence is satisfied. For non-steady state, the transient terms may be approximated using the first and second Backward Euler schemes, applicable for constant and variable step sizes. The second order Backward Euler method is a robust, implicit timestepping scheme and, in contrast with that of first order, is second order accurate in time. However, it is not monotonic and therefore inappropriate for some quantities that must remain bounded, such as turbulence quantities. Consequently, the first order Backward Euler scheme applies in turbulence equations. It is robust, as well as the second order scheme, fully implicit, bounded, conservative in time, and does not create a time-step limitation. The transient term has no bearing on the steady state solution, but it is only first order accurate in time and therefore may induce numerical diffusion in time. The first order Backward Euler scheme approximates the transient term as: ∂





v

ϕdv

∂t

=

V(ϕ − ϕ0 ) t

(2.27)

and the second order scheme as: ∂





v

ϕdv

∂t

=

V (1.5ϕ − 2ϕ0 + 0.5ϕ00 ) t

(2.28)

Common numerical techniques utilized for the completion of the advection term discretisation are the first order Upwind Differencing Scheme (UDS), the central differencing scheme (CDS) and the high resolution scheme (HRS) (Chung, 2002). Difference schemes in their majority are based on series expansion approximations (such as the Taylor series) for continuous functions. UDS is very robust, namely numerically stable, and is guaranteed to not introduce non-physical overshoots or undershoots in the solution. In contrast with central differencing scheme, UDS takes into account the flow direction when determining the value at a cell face being able to produce reasonable results even though strongly convective flows occur in the domain. HRS scheme is considered more accurate than the other two schemes. However, it is possible to meet problems in which a particular scheme does not essentially lead to convergence. Therefore, the choice between different numerical schemes requires prior test-cases in order to conclude about their appropriateness to a certain application. The linear set of equations that arise by applying the finite volume method to all elements in the domain are discrete conservation equations. The system of equations can be written in the form:

 nbi

anb i = bi i

(2.29)

134

Process Safety and Environmental Protection 9 0 ( 2 0 1 2 ) 129–140

where ϕ is the solution, b the right hand side, a the coefficients of the equation, i is the identifying number of the finite volume or node in question, and nb means “neighbor”, but also includes the central coefficient multiplying the solution at the ith location. The node may have any number of such neighbors, so that the method is equally applicable to both structured and unstructured meshes. The set of these, for all finite volumes constitutes the whole linear equation system. Several solvers employ a solution strategy according to which the momentum equations are first solved for a guessed pressure and an equation for a pressure correction is obtained leading to a guess-correction process that typically requires a large number of iterations. Other codes like CFX employ a coupled solver which solves the hydrodynamic equations for pressure and velocity components as a single system. This solution approach uses a fully implicit discretisation of the equations at any given time-step and, as a result, reduces the number of iterations required for calculating the transient solutions in a time-dependent analysis. The linearized system of discrete equations described above can be written in the general matrix form: [A][ϕ] = [b]

(2.30)

where [A] is the coefficient matrix, [ϕ] the solution vector and [b] the right hand side. The above equation can be solved iteratively by starting with an approximate solution, ϕn , that is to be improved by a correction, ϕ , to yield a better solution, ϕn+1 , i.e. ϕn+1 = ϕn + ϕ

with rn , the residual, obtained from (2.33)

Repeated application of this algorithm will normally yield a solution of the predefined error level.

3.

Turbulence modeling

The characteristic physical process that fluid flow in realconditions presents is the turbulent motion. One of the most prominent turbulence prediction models implemented in many general purpose CFD codes is the k–ε eddy-viscosity model. The k–ε model has proven to be stable and numerically robust having a well established predictive capability. The k–ε model uses the scalable wall-function approach instead of standard wall functions, improving robustness and accuracy when the near-wall mesh is very fine. Eddy-dissipation models such as the Shear Stress Transport (SST) model and Reynolds stress models such as SSG model (Park and Kwon, 2004; Menter, 1994) are also available in the CFX code.

3.2.

(2.32)

rn = b − Aϕn

3.1.

Combustion and radiation transport modeling

(2.31)

where ϕ is a solution of ˙  = rn Aϕ

tion of certain conditions on the external boundaries of the domain. (7) CFX-Solver, where the code performs calculations to the final solution using iterative methods, until the predefined statistical error level has been achieved. (8) CFX-Post, in which the elaboration of the results takes place, supported by outstanding graphics capabilities (visualization of the geometry and control volumes, vector plots showing the direction and magnitude of the flow, visualization of the variation of scalar variables with time through the domain). Quantitative calculations are also performed in this stage.

CFD code structure

CFD codes are normally divided into three parts: preprocessor, solver and post processor for problem definition, problem solution and results processing. The CFX code, which was utilized for the research purposes of the project, consists of four sequential parts: (5) CFX-Build, in which the computational domain geometry is constructed as a unified set of parametric surfaces and it is divided into a number of smaller control volumes (cells) through the mesh generation process. (6) CFX-Physics Processor, where boundary conditions, physical models, domain fluids and solution scheme are defined. The equations used for describing the fundamental transport phenomena of a particular problem were the Reynolds Averaged Navier–Stokes (RANS) equations that express the conservation of mass, momentum and energy. The equations relating to the flow (fundamental and sub-models) are numerically defined by the specifica-

In accident scenarios that include fire (for instance, flash fire or fireball), the need for combustion reaction and radiation transport modeling emerges. The equation of transport for any component of the reacting mixture is described by its individual mass balance equation, in which the chemical reaction rate of production/consumption of the component is included. The rate of progress for a reaction may be calculated using the eddy dissipation model, which is based on the concept that chemical reaction is fast relative to the transport processes in the flow. This is valid for diffusion flames, transient development of fires and flares (Magnussen, 1989). Available techniques for thermal radiation transport computing are the Differential Approximation or P1 model, the Discrete Transfer Radiation Model (DTRM), the Rosseland model and the Monte Carlo model. In this paper, the approach to thermal radiation transport due to fuel gas combustion has been accomplished through the P1 model, which actually is a simplification of the Radiation Transport equation (Seaid et al., 2004).

4.

CFD simulation in consequence analysis

4.1.

Isothermal releases

Isothermal releases involve release and dispersion of a gas at its physical state. What determines the formulation of the dispersing cloud is not the heat transfer between the cloud and the environment, but the gas density. Consequently, it is justified to consider heavy-gas releases as the most important from the safety viewpoint, because the forming cloud of a heavy gas remains at low height until to be diluted from the atmospheric air.

Process Safety and Environmental Protection 9 0 ( 2 0 1 2 ) 129–140

135

Fig. 1 – Computed vs. experimental concentration histories received on the back face of the obstacle for the downwind dispersion Thorney trial no. 26.

In this work, the CFD simulation of the isothermal Thorney Island trials – phase II (McQuaid and Roebuck, 1985) is presented, which concerned the dispersion of heavy gas in obstructed terrain. CFD results showed that gas concentration histories may be efficiently approximated with reasonable predictions of peak concentration and cloud passage duration, even when predictions are made upwind (i.e. opposite to wind blow direction) (Sklavounos and Rigas, 2004a). Comparisons between simulation results of different turbulence models (k–ε, k–␻, SST and SSG Reynolds) revealed that all models reproduced concentration histories in good agreement with the experimental data. However, k–ε model appeared to be more accurate. In Fig. 1, computational concentration histories for downwind dispersion trial no. 26 are compared vs. the experimental data received at the back face of the obstacle.

4.2.

Non-isothermal releases

A different type of dispersion is that resulting from cryogenic releases. What dominates after the spill is the violent vaporization of the liquefied substance, which leads to the formation of a cold, denser-than-air, vapor cloud (Nielsen, 2002). In such extreme conditions, the cloud absorbs heat from the ground and atmosphere dispersing as heavy rather than light gas for considerable distances before it lifts-off, although the released gas may be much lighter than air at normal conditions (Sklavounos and Rigas, 2005). Consequently, in case of releases of liquefied flammable gases such as hydrogen or natural gas, the risk of cloud ignition is much higher due to its long-lasting stay at low height. A qualitative view of the

latter is given in Fig. 2, where the dispersing cloud profile is compared considering pressurized (non cryogenic) and liquefied (cryogenic) hydrogen release (Rigas and Sklavounos, 2005a). In 2010, Health and Safety Laboratory (HSL) in collaboration with the US Fire Protection Research Foundation (FPRF) accomplished a research project concerning the evaluation of source term dispersion models especially for liquefied natural gas releases (HSE, 2010; FPRF, 2007). Some of them have been used in earlier risk analysis studies by the authors (Rigas and Sklavounos, 2002, 2004), as well as in validation studies (Sklavounos and Rigas, 2006a; Carissimo et al., 2001). A number of experimental works has been implemented focusing on the study of cryogenic releases as described by Thyer (2003). Among the most comprehensive experimental datasets in this category are those obtained from liquefied hydrogen experiments (Witcofski and Chirivella, 1984) and liquefied natural gas trials in Burro Series and Coyote series experiments (LLNL, 1984, 1983, 1982). In Fig. 3, the experimental gas concentration histories vs. the computational ones for LH2 trial no. 6 are shown. Apparently, computational concentration histories are found in good agreement with the experimental ones. The discrepancies are mainly associated with the peak overpressures, which seem to appear instantly during the experiment. Through the application of particular statistical performance measures: Fractional Bias (FB), Geometric Mean Bias (MG), Geometric Mean Variance, Mean Relative Square Error and Normalized Mean Square, it was inferred that simulation results were in good agreement with the experimental ones and within a

136

Process Safety and Environmental Protection 9 0 ( 2 0 1 2 ) 129–140

Fig. 2 – Downwind transverse cross-section 30 m from source demonstrating concentration distribution within the cloud for pressurized hydrogen release scenario after 9.6 s (a) and cryogenic hydrogen release scenario after 14.2 s (b). factor of two of the observed values (Sklavounos and Rigas, 2005). A comparative evaluation of the reliability between the CFD model and the heavy gas dispersion models SLAB and

DEGADIS was implemented, using the dispersion data from Coyote trials nos. 3, 5, 6 and 7 and calculating statistical performance measures MG and VG (Fig. 4). Actually, the model built in CFX code demonstrates statistically insignificant tendency for overprediction of gas concentration, since the MG value lies close to unity (Sklavounos and Rigas, 2006a).

4.3.

Fig. 3 – Computed and experimental hydrogen concentration history at the position (9.1, 0, 1) for the LH2 spill trial no. 6.

Cloud fires

Flash fires and fireballs are considered as cloud fires. Flash fires may be initiated either by ignition of a flammable vapor cloud formed from an instantaneous release, or by delayed ignition of a cloud from a continuous release. They pose severe risk due to the emitted thermal radiation and the relatively large duration (of the order of seconds) of the resulting fire (HSE, 1996). The instantaneous or continuous releases considered in risk studies would physically correspond to a spreading transient puff or a long steady-state plume (Pula et al., 2005). In Coyote series trials, finite duration LNG releases were performed in the open-air in order to study experimentally the evolution of cloud fire after its ignition. CFD simulations of the trials showed that cloud burning and resulting thermal emissions may be efficiently approached by means of computational fluid dynamics utilizing the P1 model for radiation transport and the eddy dissipation model for combustion computations (Rigas and Sklavounos, 2006).

137

Process Safety and Environmental Protection 9 0 ( 2 0 1 2 ) 129–140

Table 1 – Combined values of the statistical performance measures for the computed predictions of thermal impulse in Coyote trials nos. 3, 5, 6 and 7. Statistical measures Ideal value Combined values

FB

MG

0 −0.3

1 0.74

Fig. 4 – Geometric mean (MG) vs. geometric variance (VG) plot: relative ranking of CFX, DEGADIS and SLAB performance combined from dispersion predictions of Coyote trials nos. 3, 5, 6 and 7.

The application of statistical performance measures to computational predictions and experimental measures of thermal impulse gave the results shown in Table 1. The negative value of fractional bias tabulated in the second column clearly states that the thermal impulses were overestimated, so that computational predictions may be viewed as conservative. In order to check the reliability of numerical modeling predictions, the statistical measures of geometric mean bias (MB) and geometric mean variance (VG) may be applied in a manner similar to that in dispersion modeling performance assessment. The calculation of the statistical quantities of Mean Relative Square Error (MRSE) and Normalized Mean Square Error (NMSE) show that their values are close to zero, further confirming the statistical validity of the computational results. Research on modeling fireball phenomenon is also currently curried out by the authors with reasonable quantitative estimations of the thermal load received in the vicinity, whereas a good qualitative approach is accomplished for the formation and the evolution of the fireball after burning initiation (Fig. 5).

4.4.

Confined and unconfined explosions

The more severe hazard of an explosion is the development of a shock wave (for condensed explosives and under specific conditions for gas explosions) or a pressure wave (for the majority of gas and dust explosions) (Baker and Tang, 1991). Shock waves generate noticeably higher overpressures, thus leading to more extensive damage. The most common approach to estimating the consequences of explosion events is the use of empirical data from explosion records and obser-

VG 1 1.26

MRSE

NMSE

0 0.22

0 0.29

vations, which correlate certain overpressure values with a particular degree of damage. Pertinent quantitative criteria may be found in the literature (CCPS, 2000). TNO multi-energy, Baker–Strehlow, and TNT equivalency are well-established models on estimating the effects of potential gas explosions. Actually, only the latter can be applied as far as dense explosives are considered (Philips, 1994). This method assumes that the energy being released is equated to a mass of TNT, which would give an equivalent amount of damage in case of explosion. However, the assumption of an equivalent mass of TNT does not take into account the individual explosive properties of the material, which in fact, play an important role in blast wave characteristics; for example, higher explosive density entails higher detonation velocity and hence greater rate of energy release, which results in increased overpressures (Mader, 1996). Recent approaches of explosive phenomena using numerical techniques have shown good performance in simulating the overpressure waves produced from an explosion (Miura et al., 2004). CFD simulations were performed for the experimental trials of unconfined field-scale explosion in obstructed terrain (HSE, 2001) and confined lab-scale explosion in branched tunnel (Binggeli et al., 1993). Comparing the computed estimations for overpressure histories obtained from the 3-D numerical model built in CFX code, with the experimental ones, it was inferred a reasonably good agreement for peak overpressures and specific impulse. Through the calculation of the specific impulses, which are equal to the integral of the overpressure curve (Table 2), it was revealed that specific impulses were always overestimated with an average error of about +12% and +15%, respectively, namely the predictions may be viewed as conservative (Rigas and Sklavounos, 2005b; Sklavounos and Rigas, 2004b). Furthermore, the numerical model developed for the simulation of explosion phenomena was utilized for the design study of explosion protection systems against confined explosions. In particular, the quantitative estimation of the protective effect of explosion relief vents in case of confined explosions inside tunnels was addressed. A series of virtual experiments performed by computer simulation, revealed how the number of vents, the ratio between the main tunnel diameter (Dc /Dt ) and the branch diameter, as well as the angle between the vents and the tunnel (ϕ), influence the blast wave attenuation. In fact, a number of virtual experiments showed that the use of either horizontal, or vertical branched vents, provides an effective method for the attenuation of a shock wave following an explosion (Fig. 6), whereas the statistical elaboration of the results revealed that the attenuation is mainly affected by the number of vents and their diameter. Strangely, the angle between the side vents and the main tunnel appeared to slightly affect the explosion wave pressure decrease. Eventually, the influence of the above parameters was quantitatively illustrated in functional diagrams, so that the total attenuation effect may be promptly estimated, if the design variables are known. In addition, two regression models with reasonable fitting to the calculated data were constructed, which express

138

Process Safety and Environmental Protection 9 0 ( 2 0 1 2 ) 129–140

Fig. 5 – . Fireball development after release and ignition of 1708 kg of propane. The burning cloud moves upwards obtaining positive buoyancy due to high temperatures, whereas air inflow from the left (inlet boundary condition) causes horizontal fireball shift to the right (a: 1.2 s after ignition, b: 3.2 s after ignition).

Table 2 – Relative errors between computed and experimental specific impulses for confined and unconfined explosion test cases. Gauge

Experimental value (kPa s)

Unconfined (open terrain) explosion Nо. 1 Nо. 2 Nо. 3 Nо. 4 Nо. 5 Average Confined (branched tunnel) explosion Nо. 1 Nо. 2 Nо. 3 Average

Computational value (kPa s)

Relative error (%)

46 285 30 408 19 386 17 298 16 819

50 543 34 392 21 518 19 114 19 931

+9.2 +13.1 +11.0 +10.5 +18.5 +12.5

1335 1726 721

1273 1553 558

+4.9 +11.1 +29.2 +15%

Process Safety and Environmental Protection 9 0 ( 2 0 1 2 ) 129–140

139

Fig. 6 – . Gradual attenuation of the shock wave while crossing the branches (case: vertical vents, ϕ = 45◦ and Dc /Dt = 0.75). the attenuation effect as a dependent variable of the design variables and their interactions (Sklavounos and Rigas, 2006b).

5.

Concluding remarks

A series of CFD-based simulation tests were presented in this paper. The development and validation of numerical models focused on the simulation of specific types of major accidents and the quantitative estimation of physical quantities used in consequence analysis. The overall research work implemented so far showed that a series of characteristic accident types may be simulated and the associated consequences be estimated through CFD modeling, in particular: • Dispersion simulations for isothermal and non-isothermal releases provided concentration predictions within a factorof-two being more accurate than that provided from semiempirical, heavy gas dispersion models, as revealed through comparison tests. • Explosion simulations for open, obstructed field and confined spaces of complex geometry, provided specific impulse estimations with an average error into the range between +12 and +15%. • Cloud fire simulations involving flash fires and fireballs provided radiative heat load estimations within a factor of two of the observed values. In addition, the application of validated CFD explosion model in multi-parametric safety design studies proved that CFD approach may provide an alternative mean for avoiding the risk and the high cost of real experiments. Indeed, the quantitative estimation of the protective effect of sided vents in tunnel systems was accomplished, which led to the proposal of simple, one-way equations for prompt use in practice.

Acknowledgements The authors gratefully acknowledge the European Federation of Chemical Engineering that awarded our research work with the triennial EFCE Excellence Award in Process Safety 2010.

References Baker, W.E., Tang, M.J., 1991. Gas, Dust and Hybrid Explosions. Fundamental Studies in Engineering 13. Elsevier, Amsterdam. Binetti, R., Costamagna, M., Marcello, I., 2008. Exponential growth of new chemicals and evolution of information relevant to risk control. Annali dell’Istituto Superiore di Sanita 44 (1), 13–15.

Binggeli, E., Binggeli, F., Schlapfer, D., Bucher, K.M., 1993. Airblast predictions in tunnel/entrance configurations due to he-detonations near the tunnel portal. In: Proceedings of the 13th International Symposium on Military Applications of Blast Simulation , The Hague, Netherlands. Carissimo, B., Jagger, s., Daish, N., Halford, A., Selmer-Olsen, S., Riikonen, K., Perroux, J., Würtz, J., Bartzis, J., Duijm, N., Ham, K., Schatzmann, M., Hall, R., 2001. The SMEDIS database and validation exercise. Journal of Environment and Pollution 16, 614–629. CCPS, 2000. Guidelines for Chemical Process Quantitative Risk Analysis. AIChE, New York. Chung, T.J., 2002. Computational Fluid Dynamics. Cambridge University Press, Cambridge. EEC, 2003. Directive 2003/105/EC. European Community, Brussels. EEC, 1996, Directive 96/82. European Community, Brussels. FPRF, 2007. Evaluating vapor dispersion models for safety analysis of LNG facilities. Report MSU/2007/04. FPRF, Massachusetts. HSE, 2010. LNG source term models for hazard analysis. Report 789. Health and Safety Executive, Sheffield. HSE, 2001. Explosion hazard assessment: a study of the feasibility and benefits of extending current HSE methodology to take account of blast sheltering. Report HSL/2001/04. Health and Safety Executive, Sheffield. HSE, 1996. Review of flash fire modeling. Report 94/1996. Health and Safety Executive, Sheffield. LLNL, 1982. Burro Series Data Report Vol. 2, LLNL/NWC Report UCID: 19075. LLNL, Livermore. LLNL, 1983. Coyote Series Data Report Vol. 1, LLNL/NWC Report UCID: 19953. LLNL, Livermore. LLNL, 1984. Vapor Burn Analysis for the COYOTE Series LNG Spill Experiments, LLNL/NWC Report UCID: 53530. LLNL, Livermore. Mader, C.L., 1996. Numerical Modeling of Explosives and Propellants. CRC press, New York. Magnussen, B.F., 1989. Modeling of NOx and soot formation by the eddy dissipation concept. In: Proceedings of the First Topic Oriented Technical Meeting , International Flame Research Foundation, Amsterdam. McQuaid, J., Roebuck, B., 1985. Large scale field trials on dense vapour dispersion. Report EUR 10029 (EN), Health and Safety Executive, Sheffield. Menter, F.R., 1994. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA Journal 32, 1598–1605. Miura, A., Mizukaki, T., Shiraishi, T., Matsuo, A., Takayama, K., Nojiri, I., 2004. Spread behavior of explosion in closed spaces. Journal of Loss Prevention in the Process Industries 17, 81–86. Nielsen, M., 2002. Spreading of cold dense clouds. In: Fingas, M. (Ed.), The Handbook of Hazardous Materials Spills Technology. McGraw-Hill, New York, USA, pp. 18.1–18.15. Park, S.H., Kwon, J.H., 2004. Implementation of k–ω turbulence models in an implicit multigrid method. AIAA Journal 42, 1348–1357. Philips, H., 1994. Explosions in the Process Industries. IChemE, Warwickshire.

140

Process Safety and Environmental Protection 9 0 ( 2 0 1 2 ) 129–140

Pula, R., Khan, F., Veitch, B., Amyotte, P., 2005. Revised fire consequence models for offshore quantitative risk assessment. Journal of Loss Prevention in the Process Industries 18, 443–454. Rigas, F., Sklavounos, S., 2006. Simulation of Coyote series trials – Part II: A computational approach to ignition and combustion of flammable vapor clouds. Chemical Engineering Science 61, 1444–1452. Rigas, F., Sklavounos, S., 2005a. Evaluation of hazards associated with hydrogen storage facilities. International Journal of Hydrogen Energy 30, 1501–1510. Rigas, F., Sklavounos, S., 2005b. Experimentally validated 3-D simulation of shock waves generated by dense explosives in confined complex geometries. Journal of Hazardous Materials 121, 23–30. Rigas, F., Sklavounos, S., 2004. Major hazards analysis for populations adjacent to chemical storage facilities. Process Safety and Environmental Protection 82, 341–351. Rigas, F., Sklavounos, S., 2002. Risk and consequence analysis of hazardous chemicals in marshalling yards and warehouses at Ikonio/Piraeus harbor, Greece. Journal of Loss Prevention in the Process Industries 15, 531–544. Seaid, M., Frank, M., Klar, A., Pinnau, R., Thommes, G., 2004. Efficient numerical methods for radiation in gas turbines. Journal of Computational and Applied Mathematics 170, 217–239. Sklavounos, S., Rigas, F., 2010. Advanced multi-perspective computer simulation as a tool for reliable consequence

analysis. In: Proceedings in the 13th International Symposium on Loss Prevention and Safety Promotion in the Process Industries , Brugge. Sklavounos, S., Rigas, F., 2006a. Simulation of Coyote series trials – Part I: CFD estimation of non-isothermal LNG releases and comparison with box-model predictions. Chemical Engineering Science 61, 1434–1443. Sklavounos, S., Rigas, F., 2006b. Computer aided modeling of the protective effect of explosion relief vents in tunnel structures. Journal of Loss Prevention in the Process Industries 19, 621–629. Sklavounos, S., Rigas, F., 2005. Fuel gas dispersion under cryogenic release conditions. Energy and Fuels 19, 2535–2544. Sklavounos, S., Rigas, F., 2004a. Validation of turbulence models in heavy gas dispersion over obstacles. Journal of Hazardous Materials 108, 9–20. Sklavounos, S., Rigas, F., 2004b. Computer simulation of shock waves transmission in obstructed terrains. Journal of Loss Prevention in the Process Industries 17, 407–417. Thyer, A.M., 2003. A review of data on spreading and vaporization of cryogenic liquid spills. Journal of Hazardous Materials A99, 31–40. Witcofski, R.D., Chirivella, J.E., 1984. Experimental and analytical analyses of the mechanisms governing the dispersion of flammable clouds formed by liquid hydrogen spills. International Journal of Hydrogen Energy 9, 425–435.