Abstract: Gas flaring, a long-established but unacceptable practice in the Nigerian Petroleum industry, has deleterious effects on the environment. Research efforts have concentrated on using empirical formulae derived from experiments in predicting heat fluxes around a gas flare. However numerical method, which is capable of eliminating the severe inaccuracies in calculation from empirical formulae, is yet to be explored. The flare was considered as a three-dimensional turbulent jet issuing into a continuous cross flow of air. The momentum and scalar fluxes were approximated by the k-ε turbulence model and the resultant conservation equations were solved using the finite volume technique. The laminar flamelet concept was used to characterize the local thermochemical state of the combusting mixture, while the discrete transfer method was used to compute the radiative heat flux. A model for an accurate prediction of the heat flux at any point within the vicinity of a gas flare was developed. Thermal radiation profiles predicted at ground level exhibited Gaussian behaviour. Deviation of the predicted flux result at ground level estimated by the empirical formulae were between 15.0 and 35.0%. Based on the existing standard for maximum allowable radiation flux of 6.3 kW m^-2, the model minimum safest distance from the flare was 1.1 m, whilst the empirical formulae gave minimum safest distances ranging between 0.8 and 1.4 m. This model had superior prediction capabilities than the existing empirical formulae. Safe distances could be deduced from heat flux standpoint for both humans and habitat using the model.

INTRODUCTION

Nigeria is a country richly blessed with natural gas and it ranks 9th in the world and 2nd in Africa in its gas reserves (Abdulkareem and Odigure, 2006). Her oil wealth has been exploited for >45 years, but while oil companies have gained tremendously from the resources, communities living around it have suffered from the daily pollution caused by the by-products of non-stop gas flaring. Nigeria has an estimated 145 million cubic meter of proven gas but due to the absence of a natural gas exploration policy and the fact that most of the gas encountered are associated gas in oil exploration, the country flares 75% of the gas it produces and re-injects only 12% needed to enhance oil recovery (Abdulkareem and Odigure, 2006). This high percentage of gas flared is mainly due to lack of an efficient regulatory framework, poor access to local and international energy markets and financing constraints for projects aimed at reducing gas flaring.

Nigeria is known to flare more gas than anywhere else in the world (Fig. 1). It is estimated that about 1.6 million cubic meter of gas is being flared in Nigeria (Abdulkareem and Odigure, 2006), the highest in any member nation of OPEC. In Western Europe 99% of associated gas is used or re-injected into the oil reservoir.

Fig. 1:	Amount of gas flared/vented worldwide in 2004

But in Nigeria, despite regulations introduced >20 years ago to outlaw the practice, most associated gas is flared, causing local pollution and contributing to climate change.

Flaring in Nigeria was made illegal under regulations in 1984 and only allowed in specific circumstances on a field-by-field basis pursuant to a ministerial certificate. The Nigerian government has not enforced environmental regulations effectively because of the overlapping and conflicting jurisdiction of separate governmental agencies governing petroleum and the environment as well as because of non-transparent governance mechanisms (Kaldany, 2001). Neither the Federal Environmental Protection Agency (FEPA), nor the Department of Petroleum Resources (DPR) has implemented antiflaring policies for natural gas waste from oil production, nor have they monitored the emissions to ensure compliance with the regulation.

The large volume of gas being flared makes the process of solution gas flaring an environmental, ecological, economic and social concern for the public, oil and gas producers and regulatory agencies. The principal problem in flaring is one of safety. Personnel and sensitive structures must be protected from the intense thermal radiation emitted by the flare. There is a limit to which fluxes released can be tolerated. Typical levels by Bjorge and Bratseth (1996) are 1.6 kW m^-2 for continuous exposure, 4.7 kW m^-2 for exposure of limited duration and 6.3 kW m^-2 for short time (<1 min) exposure. Brzustowski and Sommer (1973) and Schwartz (1977) suggested that people with appropriate clothing can tolerate a radiant heat intensity of about 4.73 kW m^-2 for several minutes. The several minute time spans is satisfactory for situations, where a worker is infrequently in the flare area or must go into the flare area briefly to take some action. Oenbring and Sifferman (1980a) recommended allowable fluxes of 600 Btu h^-1 sq ft (1.89 kW m^-2) for continuous exposure time, 750 Btu h^-1 sq ft (2.37 kW m^-2) for 2 h, 900 Btu h^-1 sq ft (2.84 kW m^-2) for half to 1 h, 1,100 Btu h^-1 sq ft (3.47 kW m^-2) for 5-10 min, 1,400 Btu h^-1 sq ft (4.42 kW m^-2) for 2-5 min and 2,000 Btu h^-1 sq ft (6.31 kW m^-2) for 15 sec (escape only). The Nigerian Department of Petroleum Resources (DPR) has set a maximum limit of 6.31 kW m^-2 for a limited duration.

Over the years, the heat flux from flares has been obtained either using empirical formulae developed from laboratory experiment or by direct measurement from field during flaring. Although, empirical formulae are easy and convenient for performing repetitive calculations, they provide little or no information about the physical and chemical processes that govern the radiative output of a flare (Fairweather et al., 1991). They provide different results for the same location to the extent that some overestimated up to 50% of the actual value. This is very dangerous and could result in wasteful resources because it is one of the major factors used to determine flare stack. The effect of overestimating results in a taller than required stack and increase cost, whilst underestimating results in a shorter-than-required stack, which exposes personnel and equipment to potentially dangerous radiation levels. Thus, it is very important to predict radiation from flares as accurately as possible.

Cook et al. (1987a) performed fifty-seven field scale experiments and compared their results with alternative data derived from laboratory scale experiments and also with results obtained from empirical prediction methods described in American Petroleum Institute recommendation for the design of flaring systems. Comparisons of flame length and trajectory revealed significant differences between the field scale data and results derived from recommended prediction methods although, on average, recommendations for the global radiative characteristics of flares, quantified in terms of fraction of heat radiated were in reasonable accord with experimental data.

Several researchers have derived mathematical models based on solutions of the fluid dynamics equations in order to provide predictions of the structure of and radiation associated with a turbulent reacting jet in a cross-flow.

Botros and Brzustowski (1979) based structure predictions on finite-difference solutions of the three dimensional elliptic equations of fluid flow, describing turbulence by the k-ε model and the non-premixed combustion process by a flame sheet model, which assumed chemical reaction to proceed via a fast single irreversible step. It was however, found necessary to assume entrained oxygen to be composed of a reacting and an inert fraction in order to obtain agreement between model predictions and experimental data.

Galant et al. (1984) produced a complete model of the radiation received around a flare, based on solutions of a pseudo-stream function formulation of the fluid flow equations. The overall model embodied a k-ε turbulence model coupled with a combustion model based on the rate of turbulent mixing approach of Magnussen (1981). Mechanisms of soot formation and oxidation were included in order to permit application to large-scale flare and radiative heat transfer within a flame modeled using an adaptation of the six-flux approach of Gosman and Lockwood (1973). Predictions of the radiation received around a flare, obtained with an equivalent emitting body representation of the flame, were found to be in good agreement with data obtained from field scale experiments performed under low wind speed conditions.

Fairweather et al. (1991) described a model based on finite-difference solutions of the fluid dynamic equations of three-dimensional elliptic flow, with a two-equation k-ε model to calculate the distribution of Reynolds stresses. He described the turbulent nonpremixed combustion process using the fast chemistry assumption together with a specified probability density function (pdf) for a single conserved scalar variable (mixture fraction) coupled to solutions of modeled transport equations for the mean and variance of that variable. In addition, the laminar flamelet concept was used to determine the instantaneous thermochemical state of the combustion mixture as a function of the conserved scalar with flamelet prescriptions obtained from laminar counterflow diffusion flame calculations that employed a global reaction scheme for hydrocarbon combustion. Predictions of the complete model for flame trajectory and length were found to be in reasonable qualitative agreement with experimental data on natural gas flames from wind tunnel studies of Birch et al. (1989) over a range of cross-flow to velocity ratios.

Fairweather et al. (1992) described an extension of the model of Fairweather et al. (1991) to cover large-scale soot fires through the incorporation of the models for the formation and oxidation and for the radiative fluxes received by targets external to the fire. The usefulness of the complete model for predicting the consequences associated with atmospheric venting and flaring of flammable gases at sonic velocities was assessed and compared with experimental data on natural gas flame for both laboratory scale experiment of Birch et al. (1989) and field scale experiment of Cook et al. (1987a).

Cook et al. (1987a) performed different experiments to determine the size, shape and radiative characteristics of a flare in addition to measurements of thermal radiation incident about it. Cook et al. (1987b), based on the experiments performed by Cook et al. (1987a), derived a mathematical model in predicting levels of thermal radiation received about a flare using hybrid point source model in which the flare as an emitter of radiation was represented by a series of point sources, which was uniformly distributed along the flame locus. Fairweather et al. (1991) developed a model based on finite-difference solution of fluid dynamic equations for predicting the flame trajectory and length of a flare. Fairweather et al. (1992) extended the research of Fairweather et al. (1991) by developing a model for predicting the structure and received thermal radiation around a turbulent reacting jet discharging into a cross flow.

The present study is also based on the solution of fluid dynamic equations but using the finite volume method. The resulting algebraic equations were solved using the SIMPLER method of Patankar (1980). In using this method, the continuity and momentum equations were used to derive an exact equation for the pressure. With this equation, the pressure field is directly known during each iteration and this enhances convergence and reduces computer time.

The flare was considered as a three-dimensional turbulent jet issuing into a continuous cross flow of air. The flow in and around the jet as shown in Fig. 2 can be predicted by solving the full three-dimensional form of the partial differential conservation equations for mass, momentum and scalar transport. The five main features that are modeled in this procedure are the fluid dynamics conservation equations, the turbulence characteristics of the flow, the chemical reaction, soot and finally the heat transfer by radiation.

Fig. 2:	Turbulent jet in a cross flow

Governing equations: In general, flows can be described by means of the equations of motion, describing transport of mass, momentum and energy. This also holds for the jet flows considered in this study. The steady three dimensional flow field in and around the deflected jet is governed by the continuity equation, the momentum equations and the scalar transport equation. With the density-weighted averaging, the equations of the continuity, momentum and scalar transport may be written, for high-turbulence Reynolds numbers, respectively as:

(1)

(2)

(3)

Where the turbulent scalar flux is:

(4)

Where:

φ	=	A generalized dependent scalar variable
	=	Turbulent momentum flux, which for this study is approximated via the two equation
k-ε	=	Turbulence model

The k-ε model is the most widely used turbulence model, especially in the case of combustion. It has a fairly good performance in industrial flows (Luppes, 2000). The turbulent shear stress is related to the mean velocity gradients via a turbulent viscosity defined as:

(5)

Where:

	=	The turbulent viscosity
Cμ	=	0.09 (Jones and Launder, 1972)

The respective modeled transport equations for the turbulence kinetic energy, k and the rate of dissipation of the turbulent kinetic energy, ε are based upon the transport equation:

(6a)

(6b)

This model contains four empirical constants. The following values of these constants are used, which are recommended by Jones and Launder (1972):

The modeling approach used in this study is by describing the turbulent nonpremixed combustion process using the fast chemistry assumption (that is sufficiently fast for all reactions to go to completion or equilibrium as soon as the reactants are mixed), together with a specified probability density function (pdf) for the single conserved scalar variable (mixture fraction) coupled to solution of the mean and variance of the variable. The laminar flamelet concept is then used to determine the instantaneous density and temperature of the combusting mixture as a function of the conserved scalar.

In addition to the continuity, momentum equations and equations describing turbulence quantities there is the need for the balance equations for the mixture fraction and the mixture fraction variance, respectively. The conservation equation are:

(7)

(8)

Where C_D and σ_t are constants, respectively. Fluctuations of the mixture fraction are taken into account by introducing a presumed probability density function (pdf). A generally accepted form for the mixture fraction pdf is the beta function pdf:

(9)

Where the parameters α and β are defined in terms of the mean mixture fraction and varianc;

(10a)

(10b)

and the mean density ρ obtained from:

(11)

Predictions of soot level within the flame were incorporated in the model through the solution of balance equations for soot mass fraction, Y_c(s) and soot particle number density, N, using the following model transport equations:

(12)

(13)

Where the last terms in both Eq. 12 and 13 are the instantaneous rates of formation of the soot mass fraction and soot particle number density, which are given in Fairweather et al. (1992). Thermal radiation received at ground level around the flare was computed from converged flow field calculations using the Discrete Transfer Radiation Model (DTRM). DTRM is based on solving Radiative Transfer Equation (RTE) for some representative rays fired from the domain boundaries. The method is used in calculating radiative heat transfer to a receiver positioned within the finite difference solution domain by integrating the equation of radiative transfer along a number of representative rays that converge on the receiver. These rays are considered to emanate from a single point in a specific direction chosen in order to divide the circular field of view of a receiver into constant polar and azimuthal angles.

The path of each ray from that point to a boundary of the finite-difference solution domain is then tracked through the various control volumes of the domain. In this way, each ray can be divided into a number of optical path lengths, with each length corresponding to the distance traveled through any particular control volume.

The temperature and mixture composition associated with each path length is also known from results obtained from the flow field calculation for each of the homogenous control volumes. The information is used to derive the intensity of radiation arriving at the receiver along each ray, by integrating the transfer equation from the boundary of the solution domain, with the net radiative heat flux to the receiver over a prescribed field of view determined by summing the incoming intensities resolved to the unit vector central to that field of view.

Solution procedure: The governing equations for the turbulent jet flows considered in this study were given. Most of the equations are conservation equations with convection, diffusion and source terms. Some of the equations are algebraic expressions without partial derivatives. Each equation contains a number of unknown variables.

All equations together, in combination with the relevant boundary conditions, form the mathematical model that describes the ensemble of relevant physical processes and phenomena in the jet flow. Exact analytical solutions to the governing equations are extremely difficult to obtain. Only upon the introduction of gross assumptions and simplifications can analytical solutions be obtained. Such analytical solutions would only be an approximation of the exact solution. Another approach is to solve the equations numerically using a computer to obtain a sufficiently accurate solution. This is the approach employed in this research.

The equations are formulated in integral conservation form for an arbitrary control volume. The computational domain is subdivided in control volumes or cells. For each control volume, surface integrals over the boundary faces of the convection and diffusion fluxes are evaluated and the fluxes are approximated and written in terms of the discrete values of variables in neighboring control volumes. In this study, a grid in three dimensional Cartesian coordinate is employed.

The governing partial differential equations for the conservation of mass (Eq. 1), momentum (Eq. 2) and scalar quantities (Eq. 3) under steady-state can be represented in tensor notation form as:

(14)

Where:

The transport Eq. 14 is reduced to its finite-difference form by integrating over the computational cells into which the domain is divided. The resulting algebraic equations can be written in the following form:

(15)

Where:

For the three dimensional Cartesian coordinate system used the neighbours are six and will be represented by subscript N, S, E, W, T, B,. Where, N, S, E, W, T, B denote North, South, East, West, Top and Bottom, respectively of the nodal point p being considered. The a_p’s are the convection and diffusion coefficients of Eq. 13. The exact form of the a_p’s depends on the nature of the governing equations (elliptic, parabolic or hyperbolic) and on the form of finite difference approximation used (upwind, central, etc.). The source term S_φ is defined as:

(16)

Where:

The transport equations are discretized by integrating the respective equations over a control volume and using an interpolation procedure to evaluate the values of variables on the faces of each control volume.

Boundary conditions: Due to the elliptic nature of the equations, boundary conditions are prescribed for all the variables at all boundaries of the solution domain. The boundary conditions for the various faces of the calculation domain were specified as follows. The x-y plane through the origin was treated as a plane of symmetry; thus the velocity component was set equal to zero there, while the normal gradients of all other variables were taken to be zero (Neumann conditions). The other x-y boundary, the top x-z plane and the upstream, y-z plane were considered to be sufficiently far away from the jet and the velocity field there was characterized by U = U_∞, V = 0, W = 0. The bottom x-z plane coincided with the wall, where all the velocity components were set equal to zero, except that v was given the value V_j for the grid points representing the jet. At the downstream y-z plane, the normal gradients of the variables were assumed to be zero. For the boundary conditions of the turbulence quantities, the wall was treated by the wall-function method. The treatment amounts to assuming that the universal logarithmic velocity profile prevails near the wall that the diffusion of the turbulence kinetic energy is zero in that region and that the dissipation rate is inversely proportional to the distance from the wall.

Solution convergence: The model requires the solution of the steady flow, three-dimensional forms of ten-coupled partial differential equation. Global convergence is obtained when the pressure correction in the SIMPLIER method used is no longer needed, which means that the residual mass source S_c is equal to zero and hence the continuity equation is satisfied. For conservation of momentum, the residual sources of the momentum eequations should be zero. The residual sources of the other dependent equations should also be zero. For this study, the algorithm is considered to have converged when the normalized residual sources have reached a value smaller than a specified small value (10^-4).

RESULTS AND DISCUSSION

The usefulness of the model developed is to be able to predict accurately the heat flux received by personnel during gas flaring. The model equations were simulated using FORTRAN 95 programme. The heat flux obtained by the model was validated using results obtained from experiment performed and also compared with other empirical formulae.

Verification of the model: The model was validated using data from the experimental research of Birch et al. (1989). The simulation was carried out by considering a release of natural gas containing 94% methane by volume. The predicted and experimental heat flux results for the symmetry and cross-stream plane are shown in Fig. 3 and 4, respectively.

Fig. 3:	Variation of radiative heat flux on the symmetry plane (i.e., on the ground)

Fig. 4:	Variation of radiative heat flux on the cross stream plane (i.e., on the ground)

Table 1:	Properties used for comparison

The predicted heat fluxes are similar to the results obtained from experiment as shown in Fig. 3 and 4. Downwind of the maximum heat flux, it falls off more slowly because of the influence of the wind, since the wind is in that direction and the place is largely filled with either visible flame or hot combustion products dispersed in the direction of the wind.

Comparison of the model with some empirical formulae: The model was compared with four different existing empirical formulae obtained from the literature using data (Table 1) for the properties of the flared gas in the flare stack. The results (Fig. 5) show similar gaussian thermal radiation profiles at the ground level for all the models, but different value for the maximum and for distances downwind from the source. For this comparison, the heat flux predicted and that from empirical formula were plotted against the distance from the base of the flare stack.

Fig. 5:	Comparison of the model with four empirical formulae

The differences in the value of heat flux at the near field (that is close to the flare stack base) vary largely from each other and have nearly close values at the far field (at some distance away from the flare stack base).

At the near field, the maximum value predicted by the model is 15.5% higher than that of Kent (1964), 36.5% lower than that of Oenbring and Sifferman (1980b), 41.9% higher than that of De Faveri et al. (1985) and 22.5% lower than that of De Faveri et al. (1985) (surface model).

At the maximum limit of 6.31 kW m^-2 for a limited duration set by the Nigerian Department of Petroleum Resources (DPR), the model predicted a safe distance of above 75 m from the flare stack base, while Kent (1964) predicted a distance of above 60 m and that of Oenbring and Sifferman (1980b) and De Faveri et al. (1985) (surface model) a safe distance of above 105 and 80 m, respectively while the whole area is safe as for De Faveri et al. (1985).

The correct value of the heat flux is very important because in many cases, it forms the basis for determining the flare stack height and location. If underestimated it puts the workers at great risk of exposure to dangerous heat fluxes and if overestimated result to taller than required stack, which results in waste of resources.

CONCLUSION

Gas flaring is a long-established but unacceptable practice in the Nigerian Petroleum industry, with monumental negative social, economic, environmental and ecological consequences to Nigeria over the past five decades. The practice resulted in wastage of valuable resources much needed for economic development and has caused environmental degradation and health risks.

Consequently the reduction and elimination of gas flaring would benefit the nation, the local oil producing communities, the oil producer and the environment in a win-win situation.

One of the effects addressed in this study is the thermal pollution, since there is a limit to which the human body can tolerate the fluxes released during gas flaring. Furthermore, both habitat and structural buildings nearby also have heat thresholds. This impact of radiant heat from flare is a key element in petroleum industries, which requires accurate and efficient prediction methods which should not be subjected to large errors, which may result from a number of assumptions made in the empirical formulae predictions. Thus, it is very important to predict radiation from flares as accurately as possible.

The usefulness of the model is its ability to enable accurate prediction of the heat flux and temperature at any point in the vicinity of a flare. Such information is also useful in siting and determining the height of a flare stack.

The heat fluxes predicted were plotted against the distance from the base of the flare stack and were compared with that from the empirical formulae. The differences in the value of heat flux at the near field (that is close to the flare stack base) vary largely from each other and have nearly close values at the far field (at some distance away from the flare stack base).