CC BY
17
УДК 502.55
A. I. Kuptsov, R. R. Akberov, F. M. Gimranov
DEVELOPING OF COMPUTER PROGRAM USING CFD RESULTS FOR PREDICTING PROPAGATION OF HAZARDOUS GASES RELEASED THROUGH VENT STACKS
TO THE ATMOSPHERE
Keywords: gas release, gas vent stack, numerical modeling, computer program.
A physico-mathematical model allowing predictions of propagation of clouds of hazardous gases in the atmosphere appearing after gas releases through gas vent stacks is developed. Results of calculations carried out with the use of the developed model are compared versus available experimental data. For allowing urgent predictions of gas propagations, computer program is developed using results of the carried out numerical calculations.
Ключевые слова: сброс газа, свеча рассеивания, численное моделирование, компьютерная программа.
Разработана математическая модель, которая позволяет спрогнозировать характер распространения облаков опасных газов в атмосфере в результате сброса через свечу. Приведено сравнение результатов расчетов по разработанной модели с экспериментальными данными. Для срочного прогнозирования распространения газа предложена компьютерная программа на основе численных расчетов.
At chemical enterprises, for avoiding incidents in emergency situations the systems of urgent emptying of the processing equipment are installed. In case of light gases (i.e. gases having density smaller than 0.8 times density of air such as hydrogen, methane, ethylene, ammonia, etc.), it is allowed to discharge the gases through gas vent stacks (vertical pipes for gas releasing), i.e. major blowouts of the gases directly to the atmosphere are permitted. In addition, such discharges of gases from the processing equipment as well as from pipelines directly to the atmosphere can be also a part of scheduled works carried out prior to regular maintenance procedures.
However, analysis of the incidents having occurred at a number of the Russian enterprises [1, 2] and experimental data for natural gas emissions [3] have shown that such methods of equipment emptying are insecure, as in this case at the terrestrial surface toxic and/or explosive concentrations of gases can take place. Since the hazardous gases are to be released to the atmosphere, it is necessary to develop a method of accurate and fast assessments of spatiotemporal propagation of the resulting air-gas clouds, which would account for meteorological conditions (i.e. a class of atmospheric stability, wind speed, etc.) as well as for design of the gas vent stack.
At present time, physico-mathematical models based on using the CFD (computational fluid dynamics) technology together with modern CFD software packages found their way to predictions of consequences of gas propagations in the atmosphere. One of such CFD software packages is Fluent [4], which is a multi-purpose CFD package. This study also made use of the Fluent software. Such packages as Fluent are extremely expensive and it is well known that highly qualified experts are needed to properly run the packages. In addition, computational time required for calculating one scenario is too time demanding, which casts doubts on their applicability for urgent assessment of consequences of an incident and administering measures on their localization and full elimination.
An alternative to full CFD studies, from our viewpoint, is the use the exhausting number of CFD calculations obtained in advance for various cases of atmospheric conditions, gas vent stack heights, gas compositions, etc. to develop a simple and reliable tool for quick estimations of the consequences of gas discharges. One of the ways to use results of the CFD studies is constructing computer programs. The programs allow an easy-to-use method to get the results in several easy steps, as will be described in detail below. To summarize the above, in the present paper, we will show our physico-mathematical model based on using the modern CFD package Fluent, compare our CFD results obtained with the use of this model with experimental data [5] and show computer program and using it for taking urgent measures during hazardous gas discharges to the atmosphere.
Physico-Mathematical Model
Prior to presenting the mathematical equations governing the process of propagation (dispersion) of a hazardous gas in the atmosphere it is necessary to discuss the main physical processes occurring in the atmosphere together with the main assumptions with regard to the physical phenomena. Next to the surface of the Earth, there is a relatively thin layer of the atmospheric air called the atmospheric boundary layer (1.5 km deep), in which the parameters of the air are strongly dependent on the momentum, heat and mass transfer from the terrestrial surface. As a rule, the atmosphere is turbulent in the boundary layer. The main approach to dealing with turbulent (chaotic) flows is replacing the instanteous values of wind speed, temperature, pressure, density, gas content, etc. at every spatial point with sums of time-mean values averaged over some periods of time and fluctuations. Substitution of the sums to the original equations of continuity, momentum (Navier-Stokes), energy, gas transport and state gives the equations, which will be similar to the original equations. The main discrepancy between the original and final equations is the additional term in the
equations of momentum, energy and gas transport. For example, in the momentum equation the additional term is called the Reynolds (turbulent) stress term, which is similar in form to the existing viscous (molecular) stress term. Thus, similarly to the viscous stress term, the Reynolds stress term can be replaced with the product of turbulent viscosity and the rate of deformation of the time-averaged wind speed. For determining the new parameter, turbulent viscosity (which rather defines properties of the flow than properties of the fluid), several approaches are available up to now. One of the most widely used approaches for determining turbulent viscosity is using the standard k-e turbulence model [6], which will be described below.
Next, we will introduce the terms: unstable, neutral and stable atmosphere. As is known, temperature usually decreases with height in the entire troposphere (the first ~11 km of the atmosphere) including the boundary layer (~1.5 km deep). This is called the lapse rate. Depending on the lapse rate value (vertical profile of temperature), the atmosphere can be unstable (most turbulent), neutral and stable (least turbulent). Usually, instead of using the lapse rate value for determining whether the atmosphere is unstable, neutral or stable, the Monin-Obukhov length is used. The Monin-Obukhov length is negative for unstable atmosphere, infinite for neutral atmosphere and positive for stable atmosphere. Hazardous gases will propagate in the atmosphere differently depending on these stability scenarios. Here we will study the processes only within the atmospheric boundary layer.
The flow is assumed to be turbulent, compressible and steady-state (stationary). Despite the fact that the fluid flow is considered stationary, the time derivative will be preserved in all the equations, which is required for obtaining the numerical solution by the time marching numerical scheme. The initial values of all parameters (at the initial time moment) have no physical meaning and are chosen such that the converged solution is obtained at a low computational cost.
At the beginning of calculations, there will be no gas emission from a gas vent stack to the atmosphere; therefore, the medium will be pure air. After that the gas emission will take place, and the medium will be an airgas mixture. Several gases will be considered in this study. At first, the hazardous gas will be methane, so that a comparison of our CFD results with our own experimental data published earlier in [3] can be made. After that, we will consider ethylene, which is another hazardous gas.
The resulting physico-mathematical model is based upon simultaneous solution of the following equations (1)-(10) [7]:
Continuity equation:
dp d(put)
- = 0.
(1)
dt dxi
where p is time-averaged density of air (for simulation of the atmosphere with no gas emission) or of an air-gas mixture (for simulation of the atmosphere with gas emission), kg/m3; xt are Cartesian coordinates in meters with i ranging from 1 to 3 (xi =x is the coordinate in the wind direction; x2=y is the coordinate in the direction normal to the wind direction; z is the vertical
coordinate); ui are components of the time-averaged wind speed in m/s (uj=u; u2=v; u3=w). Momentum equation:
d(pur ) + d(puiuj )
dt
dp d ._£_ +—
dxr dxr
dxr
(
(1 + /it )
du, du
+
dx, dx,
- ^
duk dxk
(2)
+ Pgr
where p is time-averaged pressure, Pa; gi is acceleration due to gravity, m/s2; 5j is Kronecker's delta; ^ is molecular viscosity of air or of an air-gas mixture, kg/(m-s); is turbulent viscosity of air or of an air-gas mixture, kg/(m-s).
Energy equation:
d(ph) + d(purh) = d dt dxr dxr
Л + c,
Pr,
t
dT dx
—, (3)
where h=cpT is enthalpy, kJ/kg; cp is specific heat at constant pressure of air or of an air-gas mixture, kJ/(kg-K); T is time-averaged air temperature, K; X is thermal conductivity of air or of an air-gas mixture, W/(m-K); Prt is the turbulent Prandtl number taken as Prt = 0.85 in line with [8].
Gas transport equation:
d
dxi
d(pYs ) + d(pu,Ys )
dt
dx,
pD +
Ht
Sc
dYs
t J
dx,
+ Sy , (4)
where Ys is time-averaged mass fraction of component s of the air-gas mixture; Sct is the turbulent Schmidt number [9, 10]; D is molecular diffusion coefficient for the air-gas mixture, which is dependent on mixture composition, m2/s; SY is a source term, which determines generation or loss of the gas admixture, kg/(m3-s) in the interior of the computational domain and it is taken equal to zero in this study so that the only source of the gas is the vent stack.
Transport equation for k (turbulent kinetic energy):
d(p) + d(pkur) = d
dt dxr dxr
Г HЛ H+— ak J
grad к
(5)
+ 2HtEjEj + Gb -ps + Sk
where Gb is generation (or suppression) of turbulence by buoyancy forces, kg/m-s3; Sk is a source term, kg/(m-s3) defined as [7, 11] and
Sk =-Pg,
Ht
■ Pr,
(6)
-p
where ft is thermal expansion coefficient taken equal to 0.00367, 1/K [11].
Transport equation for e (dissipation rate of turbulent kinetic energy):
d(pe) , d(Peui) =_d_
dt dxi dxi
H + |grad e
(7)
+ C,EB- (lMtEpEp + C3G) " C^y + Se
where Se is a source term, kg/(m-s4), which is taken equal to zero in this study; tensor Ey is defined using expressions of the standard k-e model [6]; C1e, C2e, ck, cE are coefficients of the turbulence model; C3e is a coefficient determining buoyancy forces. Equation for turbulent viscosity:
+
Vt = CMP — ' (8)
s
where C^ is a turbulent coefficient.
Equation of state for pressure:
p = pRT, (9)
where R is specific gas constant for air or for an air-gas mixture, kJ/(kg-K).
Turbulence coefficients (also known as turbulence constants) are taken from [12, 13]:
CM = 0.0333 ; ak = 1.0; ae
= 1.3; C1s = 1.176 ;
C 2s = 1.92.
(10)
Geometry of the computational domain used in the CFD model is depicted in Figure 1. The domain has a regular orthogonal shape with dimensions 2000x250x250 m (lengthxwidthxheight), and the domain's spatial discretization was carried out using three-dimensional elements of a regular hexahedral shape. Dimensions of the computational domain were not chosen arbitrarily. Instead, two circumstances were accounted for here. One of the circumstances was the fact that one of the boundary conditions had to be zero concentration of the gas, whereas the second circumstance was limitation for the domain's length required to reduce the computational costs. Figure 2 shows a gas vent stack placed inside the computational domain at the distance 150 m from the inlet boundary plane.
Fig. 1 - Computational domain and computational grid for modeling of propagation of gas in the atmosphere
The boundary conditions imposed at the boundary planes shown in Figure 1 follow the conventional notations adopted in the CFD community. Justification of these boundary conditions is given in [7].
Left and Right planes: "Symmetry". Inlet plane: "Velocity Inlet", and the wind is horizontal and normal to the plane. Top plane: "Velocity Inlet", and the wind is horizontal and tangential to the plane. Outlet plane: "Pressure outlet". Bottom plane: "Wall". Upper surface of the gas vent stack (see Figure 2): "Wall" - for modeling the atmosphere with no gas emission, and "Pressure outlet" for modeling the atmosphere with gas emission.
The main feature of this model is the use of the modified values of turbulence coefficients as well as the use of the source term in the transport equation for k [7], which allows considering various options of atmospheric stability on the basis of Monin-Obukhov's similarity theory [14]. Modeling is carried out in two steps. At the first step, the atmospheric boundary layer is calculated with no gas emission into the computational domain in order to produce atmospheric parameters (vertical profiles of the wind speed, turbulent kinetic energy, vertical temperature stratification, etc.) and maintain them constant throughout the entire domain length. These parameters essentially influence spatiotemporal propagation of the released gas in the atmosphere. In other words, the atmospheric boundary layer at the beginning of calculation is modeled using such simplification that the boundary condition at the exit from the gas vent stack is set to be an impermeable solid wall ("wall"). After the first step, the atmosphere is modeled with inclusion of the gas emission; at that, the boundary condition at the exit from the gas vent stack is changed to be gas pressure or gas mass flowrate through connecting the Fluent software with our own user-defined functions (UDF).
For assessing applicability of the proposed physico-mathematical model, the use was made of the field data gathered by us in the vicinity of the working equipment and published in an earlier paper [3]. In study [3], concentration of methane was measured at the height of 1.5-2 m at distances from the gas vent stack 300, 520 and 1000 m during major blowouts carried out from the main gas pipeline (internal diameter is 720 mm; length is 0.334 km) through a gas vent stack of height 2.7 m and diameter 150 mm. During the time of measurements, the wind speed varied in the range from 1.3 to 4.3 m/s. The air temperature was 17.75°C. During the measurements, the initial pressure in the gas pipeline was 4.4 MPa and the mass flowrate through a gas vent stack was 2.73 kg/s. The gas volume, vented out per one operation, was 7119.78 m3.
Results and Discussions
Comparison of the methane concentration determined via using the above described CFD model versus experimentally measured concentration at gas releases from a gas vent stack [3] is shown in Table 1.
Fig. 2 - Computational grid on the surface (z=0 m) near location of the gas vent stack
Table 1 - Average values of calculated data obtained via using the CFD model of this study and experimentally determined data published in [3]
It can be seen that the data obtained with the use of the CFD model satisfactorily agree with experimental data. For example, at the distance of 300 m from the gas vent stack the discrepancy makes 15-17%, increasing to 50% at the point corresponding to the distance from the gas vent stack of 520 m. However, at the distance from the stack of 1000 m, discrepancy decreases and makes 27.4%. These discrepancies can be attributed to incompleteness of data on meteorological conditions present during measurements and insufficient number of measurement points in the vertical direction (measurements were made only at one height). Determination of complete vertical profiles of atmospheric parameters was not conducted during field measurements.
One of the positive features of the proposed physico-mathematical model is a possibility of establishing dependences of atmospheric parameters on time [7], and capability of modeling of non-stationary gas discharges [15, 16]. This feature becomes very crucial in modeling emissions of long durations. Besides, application of the model allows carrying out assessment of influence of the stack's design on process of effluence of gas from the gas vent stack, which is necessary for achieving success in designing of new stacks and optimization of operational parameters when operating the already existing stacks.
However, such CFD packages as Fluent are extremely expensive and the use of CFD technologies requires highly-qualified experts. In addition, computational time required for calculating one of the variants can demand time ranging from several hours to several weeks, which is unacceptable from the point of view of urgent assessment of consequences of an incident and administering measures on their localization and full elimination.
Therefore, for ensuring safety by means of discharging a hazardous gas through a gas vent stack we propose to use computer program (Figure 3) for urgent predictions of propagation of hazardous gases in the atmospheric air. For development of computer program, it is necessary to solve several independent from each other problems:
1) choosing some control points for subsequent calculations (vertical thermal stratification, wind speed value, height and diameter of the gas vent stack, pressure and temperature of a hazardous gas contained inside the processing equipment or inside the pipeline) based on the objective reality;
2) numerical modeling of discharge of a hazardous gas using the proposed physico-mathematical model and search for all possible scenarios in the entire set of all control points
3) representation of results (distances, at which hazardous concentrations occur) in a form convenient for engineers.
In database program collected 486 separate numerical calculations of discharge of ethylene through a stack.
Fig. 3 - Computer program using CFD results for predicting propagation of hazardous gases released through vent stacks to the atmosphere
Conclusions
When a gas is released through technological vent stacks under certain meteorological conditions, on the territory of an industrial site, explosive and/or toxic concentrations of the gas can occur. The proposed physico-mathematical model allows predicting propagations of hazardous gases in case of their releases through gas vent stacks. The use of computer program, developed on the basis of CFD calculations carried out in advance, allows for reasonably quick and simple determination of optimum design parameters of a gas vent stack, operational parameters, favorable meteorological conditions and secure distances in case of releases of hazardous gases through gas vent stacks.
References
1. Информационная бюллетень Ростехнадзора. Выпуск №32. 2007 год. URL:ib.safety.ru/assets/pdf/Bull_32 /Bull_32_32-41.pdf Обращение: 05.03.2015.
2. Завгороднев А. В., Мельников А.В., Сафонов В.С. Проблема обеспечения безопасности сброса газа в атмосферу на объектах транспортирования и хранения природного газа. Безопасность труда в промышленности. 2011. №. 11. С. 66-71.
3. Завгороднев А. В., Акопова Г.С., Толстова Н.С., Мельников А.В. Результаты исследований рассеивания в атмосфере организованных нестационарных выбросов газа на объектах газотранспортных предприятий. ТЕРРИТОРИЯ НЕФТЕГАЗ. 2011, № 12. С. 90-97.
4. Fluent Inc. Fluent 6.3. User's Guide, Lebanon, 2006.
5. Купцов А.И., Акберов Р.Р., Гимранов Ф.М. Влияние метеоусловий на динамику рассеивания опасного газа, сбрасываемого через технологические свечи. / А.И. Купцов, Р.Р. Акберов, Ф.М. Гимранов // Проблемы
Distance from the gas vent stack, m Data for methane concentration, mg/m3
Physico-mathematical model of this study Experimental data [3]
300 5.52 4.75
520 4.90 3.25
1000 9.25 7.26
program iriinij ÇF-Çr rrajltf tor predicting prrjMflavçn о? nazatjoi/!
Ambient air: Design parameters of the stadc Gas parameters (stack):
wind speodj m/s; diameter, m: pressure, MPa:
□ mV» m м □ Ф4
0 S ml* □ 0-2 □ Û.4
□ IO nV» Ü M 1-й
dlmobpherfc btrdtifiMtwn: height, m : temperature, K;
1 1 wistatir □ > № »3
□ nwitnil Й io □ 273
Й statte □ is □ МЭ
Distances within which the MAC exists: 17 to 10933 rti
сбора, подготовки и транспорта нефти и нефтепродуктов. - 2015. - № 4 (102). - С. 171-177.
6. B.E. Launder and D.B. Spalding, Mathematical Models of Turbulence, Academic Press, New York, 1972. http://dx.doi.org/10.1002/zamm.19730530619
7. Купцов А.И., Акберов Р.Р., Исламхузин Д.Я., Гимранов Ф.М. Численное моделирование пограничного слоя атмосферы с учетом ее стратификации. / А.И. Купцов, Р.Р. Акберов, Д.Я. Исламхузин, Ф.М. Гимранов // Фундаментальные исследования. - 2014. - № 9-7. - С. 1452-1460.
8. B. Weigand, Analytical Methods for Heat Transfer and Fluid Flow Problems, Springer, New York, 2004.
9. S. Hassid, Turbulent Schmidt number for diffusion models in the neutral boundary layer, Atmospheric Environment, 17 (1983), 523 - 527.
10. I. Yimer, I. Campbell and L.-Y. Jiang, Estimation of the turbulent Schmidt number from experimental profiles of axial velocity and concentration for high-Reynolds-number jet flows, Canadian Aeronautics and Space Journal, 48 (2002), 195 - 200.
11. J.E. Pieterse, CFD Investigation of the Atmospheric Boundary Layer under Different Thermal Stability Conditions", M.S. Thesis, Stellenbosch University, 2013.
12. C. Alinot and C. Masson, Aerodynamic simulations of wind turbines operating in atmospheric boundary layer with various thermal stratifications, ASME 2002 Wind Energy Symposium, Reno, Nevada, USA, January 14—17, (2002).
13. C. Alinot and C. Masson, k-s model for the atmospheric boundary layer under various thermal stratifications, Journal of Solar Energy Engineering, 127 (2005), 438 - 443.
14. A.M. Obukhov, Turbulence and Dynamics of Atmosphere, Gidrometeoizdat, Leningrad, 1988.
15. Купцов А.И., Акберов Р.Р., Исламхузин Д.Я., Гимранов Ф.М. Проблемы расчета рассеивания легких газов в атмосфере при их выбросах со свечи с учетом рельефа и застройки местности и атмосферной устойчивости / А.И. Купцов, Р.Р. Акберов, Д.Я. Исламхузин, Ф.М. Гимранов // Вестник Казанского технологического университета. - 2014. - № 6. - С. 284286.
16. Adel I. Kuptsov, Roald R. Akberov, Fidais M. Gimranov / Calculation of gas parameters at the exit from a gas vent stack by means of calculating duration of emptying of the processing equipment // Contemporary Engineering Sciences, Vol. 9, 2016, no. 3, 103-111.
© А. И. Купцов - аспирант каф. промышленной безопасности КНИТУ, artpb@yandex.ru; Р. Р. Акберов - канд. техн. наук, доц. каф. промышленной безопасности КНИТУ; Ф. М. Гимранов - д-р. техн. наук, проф., зав. каф. промышленной безопасности КНИТУ.
© A. I. Kuptsov - postgraduate department of industrial safety KNRTU; artpb@yandex.ru; R. R. Akberov - candidate of technical sciences, associate professor department of industrial safety KNRTU; F. M. Gimranov - doctor of technical sciences, professor, head of industrial safety chair KNRTU, expert92@mail.ru.
