RISK ASSESSMENT WITH THE USE OF THE MONTE-CARLO METHOD

1*Dep. «Hydraulics and Water Supply», Dnipro National University of Railway Transport named after Academician V. Lazaryan, Lazaryana St., 2, Dnipro, Ukraine, 49010, tel. +38 (056) 273 15 09, e-mail water.supply.treatment@gmail.com, ORCID 0000-0002-8525-7096 2*Dep. «Hydraulics and Water Supply», Dnipro National University of Railway Transport named after Academician V. Lazaryan, Lazaryana St., 2, Dnipro, Ukraine, 49010, tel. +38 (056) 273 15 09, e-mail water.supply.treatment@gmail.com, ORCID 0000-0002-1531-7882 3*Dep. «Life Safety», Prydniprovska State Academy of Civil Engineering and Architecture, Chernyshevskoho St., 24а, 49600, tel. +38 (056) 756 34 57, e-mail berlov@pgasa.dp.ua, ORCID 0000-0002-7442-0548 4*Dep. «Life Safety», Prydniprovska State Academy of Civil Engineering and Architecture, Chernyshevskoho St., 24а, 49600, tel. +38 (056) 756 34 57, e-mail cherednichenko@pgasa.dp.ua, ORCID 0000-0002-1457-9282


Introduction
The emission of chemically hazardous substances during industrial accidents poses a threat to the lives of employees of these enterprises and the population at all. In this regard, an extremely important problem is the injury risk assessment in the event of such man-made accidents. To solve the problems of this class, a number of parameters are essentially plausible, such as the intensity of the chemical substance emission. This forces to use the models representing the probability of one or another parameter, in addition to deterministic mathematical models. Therefore, the scientific direction of the mathematical models development for predicting the environmental pollution level in case of emergency emission of chemically dangerous substances is of practical interest. This makes it possible in a certain way to take into account the probability of a number of parameters that affect the formation of chemical contamination zones.
To assess the risk of human injury at chemically hazardous objects [1,2,4,[6][7][8][9], as a rule two approaches are used. They are a Gaussian model or a normative technique used in the State Emergency Service of Ukraine (SES). Based on these approaches, one can quickly determine the extent of chemical contamination zones, but they have a number of significant disadvantages, for example, they do not take into account the influence of buildings. In connection with this creation of mathematical models that allow quick determining the dimensions of chemical contamination zones and the risk of damage to humans is an urgent task [10][11][12][13]. The chemical damage risk is defined as an area where the concentration of a chemically dangerous substance exceeds the maximum allowable concentration (MAC).

Purpose
The main purpose of the work is the development of a numerical model for quick prediction of chemical contamination zones in case of accidental ammonia emission at the pumping station, as well as the development of a model for the damage risk assessment and the depth of wound in the event of fragments scattering formed during the pipeline explosion.

Methodology
Emission of chemically hazardous substances can lead to extremely negative consequencesdeath of people (staff at the industrial site, population). The risk of lethal injury to humans depends on many factors, among which the mass of a chemically dangerous substance entering the atmosphere is determinative. As it is known, this value is of probabalistic nature. However, emission limits for certain objects can be set based on the analysis of available statistical data. One can use the Monte-Carlo method to determine the mass of a chemically dangerous substance at an industrial site (for example, an ammonia pumping station). In this case, the methodology for assessing the risk of injury will be as follows: 1. To set the limits for possible release of a chemically dangerous substance at industrial site known from expert judgment (M1 is the minimum known mass of the substance; M2 is the maximum known mass of the substance).
2. To determine the most likely mass of a chemically dangerous substance 0 Q that can get into atmospheric air in the event of a possible emergency.
3. To determine the air pollution zone in the event of the release of a chemically dangerous substance in the amount 0 Q . In this zone, there is a subzone where the concentration of a chemically dangerous substance exceeds the limit value (for example, a lethal concentration).
4. To determine the number of people N who were in a subzone where a lethal concentration of a chemically dangerous substance was predicted.
Thus, solving a problem consists of two main steps: -the first is the calculation of the possible emission intensity of a chemically dangerous substance by the Monte-Carlo method; -the second is the calculation of the zones of chemical contamination with the definition of subzones of lethal injury.
To estimate the level of chemical pollution of the atmospheric air we will use the threedimensional mass transfer equation [2,3,5]: here Cis ammonia concentration; ,, u v w -vector components of the wind flow velocity; x y z -ammonia emission source coordinates;  -coefficient taking into account the chemical decay of impurity, precipitation scavenging; Q -ammonia emission rate; g w -the rate of gravity sedimentation of the impurity; t -time. Boundary conditions for the mass transfer equation are considered in the work [5].
During the calculations, we will take into account unevenness of the vertical diffusion coefficient and the air velocity in height: 1 1 where 0.15 p  ; 1 m  ; 1 0.2 k  ; 0 0.1 1 k  . For numerical integration of equation (1) we will use finite-difference methods [5]. We perform preliminary splitting of equation (1) into the sequence of solving the following equations: For numerical integration of equations (3), we use an alternating-triangular difference scheme [5] and the Euler method.
For mathematical modeling of the wind flow field at an industrial site, we use the model of potential flow. The calculation is based on equation [5]: where Р -is the velocity potential. The air flow velocity components are defined as follows: The boundary conditions for equation (4) are as follows: on impermeable boundaries and on the upper surface of the calculation area; 2) at the boundary where the flow enters the calculation area, n V -known air velocity; 3) P = const -at the outflow boundary of the calculated area.
For the numerical solution of equation (4), we use the method of total approximation, so we reduce this equation to the form: where t -time dummies. The calculated dependence for determining the velocity potential is written in the following form in two steps of splitting: The calculation according to this formula is terminated if the following condition is fulfilled: where n -iteration index (number of steps in time); ε -small number.
We determine the components of the velocity vector on the edges of the difference cells as follows: Based on the difference equations considered, a computer program was created that includes several subprograms: 1) subprogram for calculation of velocity potential field; 2) subprogram for calculation of air velocity field in the conditions of building; 3) subprograms for calculation of impurity concentration in the air for different time points after emergency emission.
For ease of use, the initial data file has been separated. The user adds to this file data on the location of buildings at the pumping station, the location of emergency emission, the intensity of the emission and other parameters.
FORTRAN was used to encode the difference equations.

Findings
The developed numerical model was used to calculate the zone of chemical contamination in case of accidental ammonia emission at the pumping station located near the Bashmachka settlement ( Fig. 1). Based on expert data analysis, it has been determined that ammonia emission can range from 200 to 500 kg. Based on the Monte-Carlo method, it is estimated that, within this range, the highest probability of accident ammonia emission (19%) corresponds to a 300 kg emission. Therefore, this emission was taken for the computational experiment. Since there is some inertia at the pump stop, it is clear that the emission will occur within a short time, so it is assumed that it will last 3 sec., that is, we have a semi-continuous emission.
The calculation is based on two approaches. The first is a calculation based only on the kinematic model (1). The second is a calculation based on two models: aerodynamics (4) and mass transfer models (1). That is, for the first calculation, the location of buildings at the pumping station is not taken into account. For the second calculation, the presence of buildings at the pumping station is taken into account.  Fig. 2 shows a chemical contamination zone of atmospheric air at the pumping station where ammonia emission takes place; the calculation is performed only based on the kinematic model. We see that chemical contamination zone has the form of "drop", no deformation of this zone is found. Fig. 3 shows the predicted zone of chemical contamination, to determine which the location of buildings at the pumping station is used.   Fig. 4 shows the dynamics of changes in the ammonia concentration at the industrial site near the building (Fig. 1, position of the receptor 2).
As we can see from Fig. 4, in the event of emergency ammonia emission, its concentration at the industrial site, near the industrial building, will be substantially higher than the MAC (MAC = 20 mg/m 3 ), i.e. there will be a risk of toxic damage to people at the site. Let us note that the computation time is 5 sec.
At the second stage of the study, the risk of injury to humans in case of fragment scattering due to explosion at the ammonia pumping station was assessed. For example, such an explosion may occur on the pipelines located at the pumping station (Fig. 5). The impact zone is defined by the range of the pipeline fragments. To simulate the process of fragment scattering in the air, Newton's second law of motion dynamics of material point was used. In vector form for a fragment having mass m, motion dynamics is modeled by the equation: where V -is the vector of fragment scattering speed; The fragment formed during explosion has a complex geometric shape. To apply model (9), we perform the following procedure. First, we determine the volume of fragment. Next, suppose that a "reduced sphere" with a radius R has such a volume. Then knowing that the volume of fragment is equal to W , we find the radius of the "reduced sphere" using the expression Further, knowing the mass of the fragment material and its initial velocity, we make calculations based on the model (9). Preliminary vector equation (9) is written in the projections on the X, Y axis (the Y axis is directed vertically up): The weight of fragment is determined as follows: ст mW    , where ст  -is the density of the fragment material.
Euler method was used for numerical integration of equations (10) - (11). The midsection is calculated as follows: The angle of fragment scattering is the initial data (Fig. 6, angle  ).
Below the figure shows the area of possible human injury in case of ejection of steel fragment having a reduced diameter of 1 cm. The fragment ejection is assumed to be at 2 m height. The initial fragment velocity is taken at 150 m/sec. An important parameter for assessment of the impact zones in case of fragment scattering is to determine the angle at which the fragment leaves the explosion zone (Fig. 6, angle  ) Fig. 7 presents the results of the calculation for the scattering angle 30  . Fig. 7. Impact zone in case of fragment scattering after the explosion at the ammonia pumping station As we can see from Fig. 7, this fragment may reach the boundaries of Bashmachka and Kalynivka settlements. That is, there is a threat of damage to people at a sufficiently large distance from pumping station.
In addition to this model, a model of fragment movement in the body of the animal, which appeared to be in the impact zone, was also developed (Fig. 8) The velocity of fragment before the "obstacle," that is, the body, is determined based on solving equations (10) - (11). The fragment "meets" the body of the animal (human) at an angle  (Fig. 8), which can be calculated in the process of solving equations (10) - (11), based on the determined values of the movement velocity components u, v and for a specific height H of fragment above the earth surface. We set the height H. Then the movement of fragment in the body begins. We model this process with the following equation of material point motion (Newton's second law): where V -the velocity motion vector of the fragment in the body of an animal (human); R Fresistance power: m -fragment weight.
Choosing the coordinate axis ОХ in the direction of fragment movement in the body, one can write the equation of movement (12) in projection on this axis: . We numerically solve this equation by the Euler method. The calculated dependence for determining the fragment velocity value in the body (at a new time step 1 n  ), is determined as follows by the Euler method: Using this difference dependency, we determine the fragment velocity in the animal body at each time step. The depth of fragment penetration into the body   xt (Fig. 8), at each time step, is determined as follows: where 0 0 x  -corresponds to the starting point of the animal body, i.e. the place where the fragment enters the body (Fig. 8);   xt -is the new position of the fragment in the body as a result of its movement.
It should be noted that dependence (15) makes it possible to estimate the wound size in the body that is to determine the damage severity to the animal or person in the first approximation. Table 1 shows the data for determining the depth of fragments penetration into the body of the animal for the damage area of 1239 m. To calculate the wound size in the body the height 1.6 Hm  is taken as the calculated one. That is, the data on the fragment penetration and its flight speed at this height were initial to calculate the fragment movement in the body. Let us note that time 0 t  corresponds to the moment of fragment impact on the body.  Table 1 we see that in about 0.03 sec., the animal's body will be wounded through.