ISSN 2307–3489 (Print), ІSSN 2307–6666 (Online)

Наука та прогрес транспорту. Вісник Дніпропетровського
національного університету залізничного транспорту, 2018, № 3 (75)

ЕКОЛОГІЯ НА ТРАНСПОРТІ

Введение

Технологические процессы предприятий авиационной промышленности включают использование, транспортировку и хранение отравляющих химических веществ (ОХВ) в сжиженном состоянии [10]. Нарушение правил эксплуатации оборудования приводит к его отказам, которые сопровождаются выбросом в атмосферу ОХВ с образованием токсичных облаков [13]. Одним из наиболее опасных видов техногенной аварии является разрушение емкости хранения или транспортировки сжиженного газа (СГ) с образованием пятна пролива [15] (рис. 1). Массовая концентрация токсичного газообразного вещества в газовоздушной смеси характеризует негативное отклонение от нормального химического состава воздуха и наряду с экспозицией является опасным параметром для обслуживающего персонала, оказавшегося в пределах зоны техногенной аварии.

Экспозиция обслуживающего персонала промышленного объекта определенным концентрациям ОХВ формирует поражающий фактор – ингаляционную токсодозу. Превышение пороговых значений токсодозы приводит к социальным последствиям – отравлению различной степени тяжести и человеческим жертвам. Поэтому определение риска предприятия для такого рода аварии является важной и актуальной инженерно-практической задачей.

Рис. 1. Развитие техногенной аварии

Fig. 1. Development of technogenic accident

Оценка последствий техногенной аварии включает в себя определение вероятности поражения обслуживающего персонала, который может быть подвергнут воздействию ОХВ, на основе математического моделирования рассеяния токсичной примеси в атмосфере [8].

Математическое моделирование физических процессов выброса и рассеяния вредной примеси в приземном слое атмосферы позволяет осуществить прогноз полей массовой концентрации ОХВ, определить ингаляционную токсодозу и вероятность поражения обслуживающего персонала [5, 17, 18].

Рис. 2. Схема вероятностной оценки
последствий аварии

Fig. 2. Scheme of probabilistic assessment
of accident consequences

Цель

Учитывая вышеизложенное, целью данной работы является разработка адекватной математической модели процесса испарения токсичной примеси с поверхности пятна пролива в результате техногенной аварии, формирования токсичного облака и его дальнейшего распространения в приземном слое атмосферы для получения пространственно-временных полей поражающего фактора ингаляционной токсодозы и определения полей вероятности поражения человека на основе пробит-анализа.

Анализ публикаций. Адекватное описание физических процессов смешения нереагирующих газов с воздухом и дальнейшего распространения смеси при истечении струи в открытое пространство или замкнутое помещение с принудительной (или естественной) вентиляцией возможно только с использованием системы нестационарных уравнений Навье-Стокса для сжимаемого газа. Ограниченные возможности современных компьютеров не позволяют эффективно осуществлять прямое численное решение этих уравнений. В настоящее время численное моделирование турбулентных течений осуществляют путем решения осредненных по Рейнольдсу-Фавру уравнений Навье-Стокса, дополненных моделью турбулентности [11]. Однако большинство моделей турбулентности не описывают с одинаковой степенью адекватности различные типы течений. Особенно это касается течений с интенсивными отрывами потока и/или большими градиентами давления температуры.

В работе [2] указано, что современные методики прогноза последствий аварий на химически опасных объектах и транспорте типа «Токси» [4, 6], «Аммиак», «SLAB» [9] реализуют модель Гаусса или аналитическое решение уравнения массопереноса и не учитывают застройку расчетного пространства зданиями. Применение численных кинематических моделей [19] для оценки территориального риска также ограничено случаями рассеивания ОХВ над ровной поверхностью. В некоторых работах учитывают сложный рельеф местности в процессе решения уравнения массопереноса конечноразностным методом [2, 9], но не берут во внимание либо трехмерный характер обтекания зданий [2], либо эффект сжимаемости течения, что не позволяет использовать эти математические модели для расчета последствий воздействия других поражающих факторов (взрывной ударной волны, теплового излучения), которые могут присутствовать одновременно при техногенных авариях на транспорте.

Кроме того, современные методики оценки загрязнения в основном базируются на детерминированном подходе [1], а при вероятностной оценке последствий поражения обслуживающего персонала на основе пробит-анализа используют зависимость вероятности от пробит-функции в табличном виде для экспертного анализа [4, 16, 18]. Это не позволяет применить данный подход автоматически в компьютерной системе для получения нестационарных полей поражающих факторов и вероятности поражения, поэтому требует усовершенствования вычислительной технологии.

В связи с этим существует необходимость построения новых математических моделей и расчетных схем для численного моделирования трехмерных течений многокомпонентных газовых смесей с учетом фактора загромождения пространства постройками, эффектов сжимаемости и химического взаимодействия. Они позволяют определять полный набор опасных параметров потока для различных сценариев развития техногенной аварии, рассчитывать поражающие факторы (в том числе и токсодозу) и строить пространственно-временные поля условной вероятности поражения обслуживающего персонала, необходимые для оценки индивидуального риска.

Методика

Постановка задачи рассеяния ОХВ. Рассмотрим формирование и движение газовой смеси на открытой промышленной площадке, на которой произошло аварийное разрушение емкости хранения сжиженного газа (рис. 3).

Рис. 3. Схема техногенной аварии:
1 – пятно пролива; 2 – поток воздуха;
3 – примесь; 4 – газовоздушное облако

Fig. 3. Scheme of technogenic accident:
1 – spillage spot; 2 – air flow; 3 – admixture;
4 – air-gas cloud

Расчетная область представляет собой параллелепипед с прямолинейными образующими, расположенными в правой декартовой системе координат (X, Y, Z) с основанием в плоскости XOZ (ось Y ориентирована в направлении, противоположном действию сил тяжести Земли). Эту область разбиваем на пространственные ячейки, причем размеры граней подбираем в соответствии с характерным размером особенностей области (шероховатостью обтекаемой поверхности, размерностью обтекаемых объектов). Под влиянием окружающей среды СГ испаряется с пятна пролива и поступает в приземный слой атмосферы с суммарной интенсивностью G_∑. Свежий воздух со скоростью ветра поступает через входную грань расчетной области, перемешивается с примесью, образуя газовоздушное облако с массовой концентрацией Q.

Для упрощения математической модели и ускорения компьютерных вычислений принимаем предположение, что физические процессы перехода выбрасываемого в приземный слой атмосферы вещества (вскипание) из жидкого состояния в газообразную фазу проходят мгновенно и в бесконечно-тонком по высоте слое. Это допущение несколько огрубляет расчеты, но в общем достаточно адекватно позволяет описать процесс попадания газовой примеси в воздух окружающего пространства.

Течение газовой смеси в расчетной области будем определять параметрами окружающей атмосферы, площадью пятна пролива, параметрами газовой примеси, поступающей в результате испарения в атмосферу. В какой-то момент времени испарение может прекратиться, и поступления примеси в область не будет.

Основные уравнения математической модели. Предполагаем, что основным фактором, влияющим на физические процессы смешения газовых примесей с воздухом и дальнейшее распространение смеси при испарении СГ в открытое пространство, является конвективный перенос массы, импульса и энергии. Поэтому достаточно использовать упрощенные уравнения Навье-Стокса, которые получены отбрасыванием вязких членов в уравнениях движения газовой смеси (эйлеров подход с источниковыми членами) [5].

Полная система уравнений, описывающая нестационарное трехмерное течение двухкомпонентной смеси газов в данной постановке имеет вид [5]:

, (1)

где – вектор-столбцы вида:

, (2)

(3)

(4)

(5)

, (6)

где – время; – составляющие вектора скорости ; – давление и плотность; – полная энергия единицы объема газовоздушной смеси:

, (7)

где – внутренняя энергия единицы массы газа; компоненты вектора; – суть проекции распределенных объемных источников; – ускорение свободного падения.

Закон переноса компоненты смеси с учетом скорости диффузии имеет вид [12]:

, (8)

где – относительная массовая плотность примеси (отношение плотности газообразного вещества примеси к плотности смеси),
– интенсивность изменения плотности примеси вследствие диффузии в соответствии с законом Фика (коэффициент диффузии определялся по методике, предложенной М. Е. Берляндом [18]).

Система уравнений (1–8) является незамкнутой. Дополним ее уравнениями, определяющими теплофизические свойства компоненты смеси. Для идеального политропного газа величина связана с и смеси зависимостью: .

Граничные условия. На входе будем задавать граничные условия на поверхностях тех граней, которые примыкают к границам расчетной области и через которые в расчетную область поступает атмосферный воздух. Набегающий поток на входе определяется величинами:

– полной энтальпии:

; (9)

– функции энтропии:

; (10)

– направлением вектора скорости потока (углами );

– относительной массовой плотностью примеси .

Параметры потока на входе определяются из равенств (9, 10) с учетом заданных с привлечением соотношения для «левого» инварианта Римана [19]. На непроницаемых участках, ограничивающих расчетную область поверхностей, выполняются условия «непротекания»: , где – вектор нормали к рассматриваемой границе. На поверхности испарения выставляется условие протекания примесного газа с заданной интенсивностью. Граничные условия на выходе будем задавать на поверхностях тех граней конечно-разностных ячеек, которые примыкают к границам расчетной области и через которые предполагается вытекание или втекание смеси. В выходных областях, кроме атмосферного давления , задаваемого либо взятого из эксперимента, использовались соотношения для «правого» инварианта Римана [20].

Начальные условия. В начальный момент времени во всех «газообразных» ячейках расчетной области принимаются параметры окружающей среды. В ячейках с испарением или истечением газа задается закон изменения расхода примеси.

Алгоритм численного решения. Векторное уравнение (1) является следствием законов сохранения массы, импульса и энергии, которые могут быть представлены в интегральной форме для каждой расчетной ячейки:

, (11)

где – объем элементарной расчетной ячейки; – ограничивающая поверхность данной ячейки, которая имеет внешнюю нормаль ; – тензор плотности потока консервативных переменных , столбцами которого являются векторы соответственно.

Закон переноса компоненты смеси (8) может быть также представлен в интегральной форме для каждой расчетной ячейки:

. (12)

Метод расчета. Компьютерное решение системы фундаментальных уравнений газовой динамики для смеси, дополненной законами сохранения массы примесей в интегральной форме, получено явным методом С. К. Годунова [19]. Для аппроксимации уравнений Эйлера применяется конечноразностная схема первого порядка. Центральные разности второго порядка используются для диффузионных источниковых членов в уравнениях сохранения примесей. Простая интерполяция давления применяется в вертикальном направлении. Метод Годунова характеризуется робастным алгоритмом, устойчивым к большим возмущениям параметров потока (например, давления).

В основе метода расчета лежит идея использования для построения разностной схемы точных решений уравнений с кусочно-постоянными начальными данными. Для гиперболической системы такие решения распадаются на совокупности независимых и сравнительно просто рассчитываемых деталей – «распадов разрывов».

Уравнения (11, 12) допускают возникновение и существование поверхностей разрыва двух видов: ударных волн и тангенциальных разрывов. Функции, удовлетворяющие уравнениям (11, 12), можно рассматривать в качестве обобщенных решений уравнений газовой динамики. Использование интегральных законов сохранения массы, импульса, энергии и концентрации газообразной примеси в качестве исходных для построения разностных уравнений обеспечивает построение разрывных решений без выделения разрывов.

Совокупность газодинамических параметров во всех ячейках в момент времени представляет собой известное решение на временном слое с индексом . Параметры в момент времени (на слое ) рассчитывались посредством применения явных разностных аппроксимаций для соотношений в рамках интегро-интерполяционного метода С. К. Годунова [20]. На первом этапе непрерывное распределение параметров заменяется кусочно-постоянными среднеинтегральными значениями в каждой расчетной ячейке. При этом границы ячейки представляют собой неустойчивые поверхности произвольного разрыва, которые распадаются на устойчивые волновые элементы: ударную волну, контактную поверхность и волну разрежения. Для каждого такого разрыва определяются потоки массы, импульса и энергии через грани газовых ячеек. Устойчивость конечноразностной схемы обеспечивается выбором величины шага по времени .

Моделирование испарения. В результате дискретизации расчетной области поверхность пятна пролива СГ также разбивается на ряд конечноразностных ячеек у земли в плоскости XOZ (рис. 4).

Рис. 4. Дискретизация пятна пролива:
1 – земля; 2 – ячейки пятна пролива; 3 – ячейки
атмосферы; 4 – суммарный расход газа через пятно пролива; 5 – расход газа через одну ячейку пятна пролива

Fig. 4. Discretization of spillage spot:
1 – earth; 2 – spillage spot cells; 3 – atmosphere cells;
4 – total gas flow through spillage stain;
5 – gas flow through one cell of spillage spot

При равномерном разбиении в направлении осей OX и OZ площади граней «испаряющих» ячеек одинаковы. Сделав допущение о равномерности потока с пятна пролива, можно определить индивидуальный заданный расход газа для каждой из ячеек «испарения» , где k – количество ячеек, примыкающих к пятну пролива сжиженного газа.

Допустим, имеется газовая среда, для которой термодинамические величины – давление P, плотность , внутренняя энергия единицы массы  – подчиняются уравнению состояния. Предположим, что в начальный момент времени t для левого полупространства X<0 среда характеризуется значениями параметров P₁, ₁, u₁, а для правого полупространства X>0 – значениями P₂, ₂, u₂ (здесь u – компонента вектора скорости в направлении координаты X, а другие ее компоненты равны нулю) (рис. 5).

Рис. 5. Расчетная схема для определения
расхода испарения газа:
1 – фиктивная вычислительная ячейка со стороны пятна пролива; 2 – граница «пролив–воздух»;
3 – воздушная расчетная ячейка

Fig. 5. Calculation model for determining
the gas evaporation rate:
1 – fictitious computational cell from the side
of spillage stain; 2 – spillage-air boundary;
3 – air calculated cell

Если привести в соприкосновение две массы газа, сжатые до различных давлений (P₁ – давление со стороны пятна пролива, P₂ – давление со стороны атмосферы), и убрать перегородку между ними, то поверхность их соприкосновения будет поверхностью разрыва в начальном распределении давления. Начальный разрыв распадается на несколько разрывов, которые с течением времени будут отходить друг от друга. На контактном разрыве испытывает скачок плотность, а значит, и внутренняя энергия (R₁, E₁ – для левой и R₂, E₂ – для правой областей), а давление P и поперечная компонента скорости U непрерывны. В свою очередь, эти области отделены от невозмущенных областей с параметрами (P₁, ₁, u₁) снизу («слева») и (P₂, ₂, u₂) сверху («справа») или ударной волной УВ, или волной разрежения ВР.

Решая задачу распада разрыва на грани конечноразностной ячейки, примыкающей к вентиляционному проему, можно определить плотность R и скорость u, а значит, и индивидуальный расход газа G_i через рассматриваемую грань. Используя метод итераций, можно подобрать давление P₂ таким образом, чтобы расчетный расход газа G_i отличался от заданного G_з на наперед заданную малую величину  (рис. 6). Тестирование такого итерационного алгоритма показало быструю сходимость процесса подбора давления «испарения» и незначительное увеличение общего времени нестационарного расчета движения газовой смеси в расчетной области. Так как информация о предыдущем шаге итерации по времени запоминается в специальной структуре данных, то итерационный процесс подбора противодавления в процессе общего расчета ускоряется.

Интерполяция функции интенсивности испарения. При моделировании испарения с пятна пролива интенсивность «выброса» примеси в газовой фазе в атмосферу обычно принимают постоянной (рис. 7).

Если имеется суммарная масса m пролитого СГ и время t₁ начала и t₂ конца процесса испарения, тогда текущая интенсивность испарения может быть найдена из соотношения:

(13)

Рис. 6. Итерационный алгоритм подбора
противодавления в текущий момент времени
процесса испарения

Fig. 6. Iterative algorithm of selection of counterpressure at the current time of the evaporation process

Рис. 7. Постоянный закон
интенсивности испарения:
t₁, t₂ – время начала и конца процесса испарения;
G_з, – заданная интенсивность испарения

Fig. 7. Constant law of evaporation intensity:
t₁, t₂ – time of beginning and end of evaporation process;
G_з, –specified evaporation intensity

Иногда закон интенсивности испарения (рис. 8) задан табличной функцией, которая получена либо из эксперимента, либо с помощью аналитической модели. В этом случае возникает проблема интерполирования таблично заданной функции. При выборе подходящего метода интерполирования следует учитывать возможный сложный характер функции и неравномерность расположения узлов интерполяции [14].

Рис. 8. Интерполирование интенсивности
испарения кусочно-кубическим сплайном:
1 – начало испарения; 2 – конец испарения

Fig. 8. Interpolation of the evaporation intensity by piecewise-cubic spline:
1 – evaporation beginning; 2 – evaporation end

Ю. К. Чернышев в работе [7] провел аналитический обзор наиболее известных методов интерполирования (метод конформных отображений, эрмитовы сплайны, сплайн Безье), рассмотрел их преимущества и недостатки и пришел к выводу, что наиболее подходящими являются кусочно-кубические эрмитовы сплайны, в основе которых лежит методика Х. Акимы построения нелинейных приближений первой производной в узлах интерполирования и ее обобщения [12]. Были введены дополнительные промежуточные узлы интерполяции и доказано, что данный прием устраняет колебания знака второй производной.

Площадь пролива (м²) можно определить по следующей формуле:

(14)

где – плотность жидкого опасного вещества, кг/м³; – суммарная масса пролитого жидкого опасного вещества, кг/м³; – масса опасного вещества, переходящая в газовую фазу в первичное облако при мгновенном вскипании перегретого опасного вещества, кг;
– масса опасного вещества, переходящая в аэрозоль в первичное облако, кг.

Скорость испарения с поверхности пролива и расход аммиака во вторичном облаке, образующемся на стадии испарения из пролива, определяем по формуле [3]:

, кг/с, (15)

где – молярная масса опасного вещества, кг/моль; – начальная эффективная скорость вторичного облака, образующегося на стадии испарения из пролива, м/с; – давление насыщенного пара опасного вещества при температуре воздуха, мм рт. ст., которое можно определить по формуле:

(16)

где – температура воздуха, К; – температура кипения жидкого опасного вещества при давлении в окружающей среде P₀ (при нормальных условиях принимается равным 101 325 Па), К; – теплота испарения (кипения) жидкого опасного вещества, Дж/кг.

Вероятностная оценка безопасности. В результате моделирования рассеяния токсичной газовой примеси в атмосфере можно получить функции изменения во времени и пространстве массовой концентрации примеси . На основании этого определяется опасность воздействия ОХВ на обслуживающий персонал (рис. 2).

Рассмотрим вычислительную технологию определения вероятности поражения. Пусть интегральная степень воздействия является случайной, удовлетворяющей нормальному закону распределения с математическим ожиданием равным 5 и дисперсией равной 1.
В этом случае величина вероятности поражения (измеряется в долях единицы) может быть оценена по следующей формуле:

(17)

С другой стороны, интегральную степень воздействия можно оценить с помощью уравнения регрессии вида , где – количественные оценки поражающих факторов; , – коэффициенты уравнения регрессии.

Задавая верхний предел интеграла (18) (пробит-функцию ), можно определить вероятность поражения. При воздействии отравляющего химического вещества (ОХВ) на персонал техногенного объекта основным поражающим фактором является ингаляционная токсодоза – интеграл по времени концентрации ОХВ в воздухе:

(19)

где – время экспозиции (время, за которое набирается ингаляционная токсодоза), с;
– пространственно-временное значение массовой концентрации ОХВ, ppm; – табличный коэффициент (например, для цианистого водорода ).

Пробит-функция для поражения человека вследствие воздействия ингаляции смертельной токсодозы отравляющего химического вещества в общем случае определяется по формуле:

(20)

где и – табличные полуэмпирические коэффициенты (например, для цианистого водорода ; ).

Обычно в инженерной практике получение вероятности поражения происходит визуально с использованием таблицы зависимости вероятности от пробит-функции поражающего фактора. Это неудобно и не позволяет использовать данный аппарат интегрированным в компьютерной системе для оценки техногенной безопасности. Поэтому для автоматизации процесса связь пробит-функции и вероятности смертельного поражения человека проинтерполируем кусочно-кубическим эрмитовым сплайном, обобщенным Ю. К. Чернышевым (рис. 9).

Результаты

Апробация математической модели. Предложенный алгоритм и метод учета переменной интенсивности испарения сжиженного токсичного газа с поверхности пятна пролива в процессе выброса газовой примеси в приземный слой атмосферы был реализован в виде подсистемы исследовательского программного комплекса «Fire». Комплекс позволяет производить трехмерный анализ рассеяния токсичных газообразных примесей во времени и пространстве в практически приемлемое время и делать прогноз рисков летального исхода вследствие воздействия испарений токсичного газа на организм человека (рис. 10, 11).

Тестирование разработанной информационной технологии и анализ эффективности алгоритма проводились на примере моделирования испарения 6925 кг сжиженного цианистого водорода (токсичного взрывоопасного вещества плотностью 689 кг/м³, молярной массой 0,027 кг/моль, температурой кипения 298,6 К, теплотой испарения 933 кДж/кг) с пятна пролива в форме окружности радиусом R, образовавшегося в результате разрушения цистерны в районе железнодорожной станции в черте населенного пункта с большими зданиями (коэффициент степенной зависимости для аппроксимации скорости в атмосферном слое над поверхностью земли k = 0,4) (рис. 10).

Рис. 9. Интерполирование зависимости
вероятности P поражения от пробит-функции P_r
поражающего фактора

Fig. 9. Interpolating the dependence of P lesion probability on the probit-function P_r of the adverse factor

Центр окружности пятна пролива располагался на расстоянии Xc = 16 м, Zc = 16 м от начала координат в расчетной области с габаритами Lx×Ly×Lz = 85×10×85 м и вариантом по количеству ячеек вдоль координатных осей 85×10×85.

На расстоянии Xa = 30 м и Za = 28 м от начала координат располагалось здание станции с габаритами Dx×Dy×Dz = 15×5×25 м.

Ветер набегает со скоростью 3 м/с под углом 45º к оси OZ на высоте 0,5 м. В этом случае начальная эффективная скорость вторичного облака, образующегося на стадии испарения из пролива, составляла 1,19 м/с.

Если пренебречь массой газовой фракции первичного облака, образовавшегося в результате аварийного разрушения емкости, и учесть, что температура кипения вещества выше температуры окружающей среды (масса опасного вещества, переходящая в аэрозоль в первичное облако  = 0), по формуле (14) можно определить площадь пятна пролива 201 м² и соответствующий радиус пятна R = 8 м. Тогда в соответствии с формулой (15) интенсивность испарения цианистого водорода с пятна пролива будет составлять 0,00106 кг/с/м².

Рис. 10. Карта объектов у земли:
1 – вектор скорости ветра; 2 – пятно пролива;
3 – здание

Fig. 10. Map of objects near the ground:
1 – vector of wind speed; 2 – spillage stain;
3 – building

Считалось, что испарение начиналось с момента времени t₁ = 0 с и принудительно прекращалось по истечению t₂ = 5 с, например, с помощью заливки пятна пролива специальной пеной. Время окончания расчетов было принято таким, чтобы газовоздушное токсичное облако покинуло пределы расчетной области.

Было выполнено два варианта расчета: 1-й – без учета наличия здания железнодорожной станции и 2-й – с учетом загромождения области расчета зданием (рис. 11, 12).

Анализ изменения полей массовой концентрации цианистого водорода проводился в плоскости XOZ в слое вычислительных ячеек у земли. На рис. 11, а представлена динамика распространения токсичного вещества без учета присутствия здания. Видно, что образовавшееся после испарения облако беспрепятственно покидает расчетную область, постепенно рассеиваясь.

Влияние загромождения расчетной области зданием существенно изменяет картину рассеивания отравляющего химического вещества (см. рис. 11, б). Токсичное облако меняет форму и размеры в соответствии с течением газовой смеси, которое вырабатывается в процессе обтекания здания. Наличие препятствия увеличивает не только максимальные значения массовой концентрации токсичного вещества, но и время прохождения облака по расчетной области, а значит и время экспозиции обслуживающего персонала вредному воздействию.

Нестационарные поля массовой концентрации цианистого водорода можно рассматривать как пространственно-временное распределение опасного параметра исследуемого физического процесса. Эти данные используются для расчета токсодозы как поражающего фактора, пробит-функции для цианистого водорода и, соответственно, условной вероятности летального исхода для человека при ингаляции данного токсичного вещества (рис. 12).

Используем для сравнения рассмотренных вариантов вычислительного эксперимента характерную площадь Sх расчетной области у земли, которую занимает зона с условной вероятностью летального исхода 50 % и выше.

В случае учета загромождения (рис. 12, б) характерная опасная зона составляла 1290 м², что существенно больше, чем в первом варианте расчета (рис. 12, а) – 1001 м2.

            а – а

б – b

Рис. 11. Поля относительной массовой
концентрации токсичной примеси в момент
времени t = 20 c после начала испарения:
а – без здания; б – со зданием

Fig. 11. Fields of relative mass concentration of toxic content at the moment t = 20 sec after
the beginning of evaporation:
а – without building; b – with building

Качественный и количественный анализ представленных полей условной вероятности летального исхода вследствие поражения человека ингаляционной токсодозой позволяет сделать вывод о существенном влиянии загромождения пространства зданием на расположение и форму опасной для обслуживающего персонала зоны, что обосновывает необходимость учета

            а – а

б – b

Рис. 12. Поле условной вероятности летального исхода человека у земли, %:
а – без здания; б – со зданием

Fig. 12. The field of conditional probability of a person’s death at the ground, %:
a – without building; b – with building

данного фактора при анализе и прогнозе развития техногенной аварии рассмотренного типа на предприятиях с целью выработки рекомендаций по снижению рисков летального исхода.

Научная новизна и практическая
значимость

В математической модели учитываются: сжимаемость потока, сложный рельеф местности (загромождение пространства зданиями), трехмерный характер процесса рассеяния, наличие испарения с пятна пролива токсичного вещества с переменной интенсивностью. Модель позволяет получать пространственно-временные распределения опасного параметра – относительной массовой концентрации токсичного газа и поражающего фактора – ингаляционной токсодозы, необходимые для автоматизированного определения нестационарных трехмерных полей вероятности поражения обслуживающего персонала на основе пробит-анализа.

Выводы

В работе представлен метод вероятностной оценки летального поражения обслуживающего персонала техногенного объекта, находящегося в зоне аварийного рассеяния облака токсичного химического вещества в условиях застройки. В основу метода положена трехмерная математическая модель процесса испарения с заданной интенсивностью сжиженного токсического газа с пятна пролива, образовавшегося в результате разрушения емкости транспортировки или хранения, и его рассеяния в приземном слое атмосферы. Модель позволяет получать поля ингаляционной токсодозы и с помощью процедуры пробит-анализа численно оценивать условную вероятность поражения персонала, который подвергается воздействию токсичного газа как показателя безопасности техногенного объекта. Дальнейшее совершенствование данного метода лежит в области рассмотрения комбинации аварийных сценариев с различными поражающими факторами.

СПИСОК ИСПОЛЬЗОВАННЫХ ИСТОЧНИКОВ

1. Беляев, Н. Н. Компьютерное моделирование загрязнения окружающей среды при разливе аммиака / Н. Н. Беляев, О. В. Коптилая // Екологія і природокористування : зб. наук. пр. – Дніпропетровськ, 2002. – Вип. 2. – С. 158–162.

2. Беляев, Н. Н. Расчет территориального риска при теракте: экспресс модель / Н. Н. Беляев, И. В. Калашников, В. А. Козачина // Наука та прогрес транспорту. – 2018. – № 1 (73). – С. 7–14.
   doi: 10.15802/ stp2018/123474

3. Мацак, В. Г. Гигиеническое значение скорости испарения и давления пара токсических веществ, применяемых в производстве / В. Г. Мацак, Л. К. Хоцянов. – Москва : Медгиз, 1959. – 231 с.

4. РД-03-26-2007. Методические указания по оценке последствий аварийных выбросов опасных веществ. – Введ. 2008–01–25. – Москва : НТЦ Промышленная безопасность, 2008. – 122 с.

5. Скоб, Ю. А. Математическое моделирование струйного истечения газовоздушной смеси с различной концентрацией примеси в атмосферу / Ю. А. Скоб // Авиационно-космическая техника и технология. – 2017. – № 4. – С. 83–92.

6. Стоецкий, В. Ф. Оценка риска при авариях техногенного характера / В. Ф. Стоецкий, В. И. Голинько, Л. В. Дранишников // Наук. вісн. НГУ. – 2014. – № 3. – С. 117–124.

7. Чернышев, Ю. К. Выпуклые векторные сплайны в применении к профилированию лопаток ГТД / Ю. К. Чернышев // Авиационно-космическая техника и технология. – 2000. – № 21. – С. 16–18.

8. Assael, M. J. Fires, Explosions, and Toxic Gas Dispersions: Effects Calculation and Risk Analysis / Marc J. Assael, Konstantinos E. Kakosimos. – New York : CRC Press, 2010. – 349 p. doi: 10.1201/9781439826768

9. Biliaiev, M. M. Numerical simulation of toxic chemical dispersion after accident at railway / M. M. Biliaiev, L. Ya. Muntian // Наука та прогрес транспорту. – 2016. – № 2 (62). – С. 7–15.
   doi: 10.15802/stp2016/67279

10. Brauer, R. L. Safety and Health for Engineers / R. L. Brauer. – Hoboken : Wiley, 2016. – 608 p.

11. Computational Fluid Dynamics for Engineers / B. Andersson, R. Andersson, L. Hakansson [et al.]. – Cambridge : Cambridge University Press, 2012. – 212 p.

12. Engeln-Müllges, G. Numerik-Algorithmen: Verfahren, Beispiele, Anwendungen / G. Engeln-Müllges, K. Niederdrenk, R. Wodicka. – Berlin : Xpert.press, 2010. – 756 p.

13. Hughes, P. Introduction to Health and Safety at Work: The Handbook for the NEBOSH National General Certificate / P. Hughes, E. Ferrett. – Kidlington, Oxford : Butterworth-Heinemann, 2011. – 608 p.

14. Knott, G. D. Interpolating Cubic Splines / G. D. Knott. – Boston : Birkhäuser, 2012. – 254 p.

15. Nolan, D. P. Handbook of Fire and Explosion Protection Engineering Principles: for Oil, Gas, Chemical and Related Facilities / Dennis P. Nolan. – Burlington : Gulf Professional Publishing, Elsevier, 2011. – 351 p.

16. Risk assessment methodology for onboard hydrogen storage / M. Dadashzadeh, S. Kashkarov, D. Makarov , V. Molkov // Intern. Journal of Hydrogen Energy. – 2018. – Vol. 43. – Іss. 12 – Р. 6462–6475.
    doi: 10.1016/j.ijhydene.2018.01.195

17. Simulation-based safety investigation of a hydrogen fueling station with an on-site hydrogen production system involving methylcyclohexane / Jo Nakayama, Hitoshi Misono, Junji Sakamoto, Naoya Kasai, Tadahiro Shibutani, Atsumi Miyake // International Journal of Hydrogen Energy. – 2017. – Vol. 42. – Іss. 15. – Р. 10636–10644. doi: 10.1016/j.ijhydene.2016.11.072

18. Skob, Y. A. Mathematical modeling of hydrogen explosion consequences at fueling station [Электронный ресурс] / Y. A. Skob, E. A. Granovskiy, M. L. Ugryumov // 7th Intern. Conference on Hydrogen Safety. – Hamburg (Germany), 2017. – P. 1–12. – Режим доступа: https://www.hysafe.info/wp-content/uploads/2017_papers/159.pdf – Загл. с экрана. – Проверено : 28.05.2018.

19. Tashvigh, A. A. Soft computing method for modeling and optimization of air and water gap membranedistillation – a genetic programming approach / Akbar Asadi Tashvigh, Bahram Nasernejad // Desalinationand Water Treatment. – 2017. – Vol. 76. – Р. 30–39. doi: 10.5004/dwt.2017.20696

20. Toro, E. F. Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction / E. F. Toro. – Berlin : Springer, 2009. – 724 p.

Ю. О. СКОБ^1*, М. Л. УГРЮМОВ^2*

^1*Каф. «Інформатика», Харківський національний аерокосмічний університет ім. М. Є. Жуковського «ХАІ», вул. Чкалова, 17, Харків, Україна, 61070, тел. +38 (057) 315 11 31,
ел. пошта yuriy.skob@gmail.com, ORCID 0000-0003-3224-1709
^2*Каф. «Інформатика», Харківський національний аерокосмічний університет ім. М. Є. Жуковського «ХАІ», вул. Чкалова, 17, Харків, Україна, 61070, тел. +38 (057) 315 11 31,
ел. пошта m.ugryumov@khai.edu, ORCID 0000-0003-0902-2735

МАТЕМАТИЧНЕ МОДЕЛЮВАННЯ НАСЛІДКІВ ВИПАРОВУВАННЯ

АВАРІЙНОГО ПРОЛИТТЯ ТОКСИЧНИХ РЕЧОВИН

НА ЗАЛІЗНИЧНОМУ ТРАНСПОРТІ

Мета. Основною метою роботи є розрахунок просторових полів розподілу умовної ймовірності летального ураження обслуговуючого персоналу залізничної станції, викликаного інгаляцією токсичного газу, який розсіюється в приземному шарі атмосфери в умовах заданої вітрової обстановки, для чисельної оцінки рівня безпеки техногенного об’єкту. Методика. Розроблено тривимірну математичну модель процесів випаровування токсичної хімічної речовини з поверхні плями пролиття в результаті аварійного руйнування ємності зберігання або транспортування зрідженого газу і подальшого розсіювання газової домішки в приземному шарі атмосфери з урахуванням захаращення простору будівлями. Також розроблено обчислювальну технологію визначення умовної ймовірності ураження людини токсичним газом на основі пробіт-аналізу ступеня впливу вражаючого фактора (інгаляційної токсодози) на організм людини. Для автоматизації обчислювального процесу табличну залежність «пробіт-функція – ймовірність ураження» замінено узагальненим кусочно-кубічним сплайном. Результати. На основі розроблених моделей отримані результати розрахунків просторово-часових полів умовної ймовірності смертельного ураження персоналу залізничної станції, який піддався інгаляційному впливу хмари газоподібного ціаністого водню. Визначено, що наявність будівель на шляху розсіювання токсичної хмари збільшує площу концентрації та час проходження хмари по розрахунковій області, що, відповідно, збільшує час експозиції обслуговуючого персоналу шкідливому впливу.
Наукова новизна. У розробленій математичній моделі враховуються: стисливість потоку, складний рельєф місцевості (захаращення розрахункового простору будівлями техногенного об’єкту), тривимірний характер процесу розсіювання газоподібної домішки, наявність випаровування з плями пролиття токсичної речовини зі змінною інтенсивністю в залежності від вітрової обстановки, фізичних характеристик домішки і класу шорсткості приземного шару атмосфери. Математична модель дозволяє отримувати просторово-часові розподіли небезпечного параметра – відносної масової концентрації токсичного газу і вражаючого фактора – інгаляційної токсодози, які необхідні для визначення нестаціонарних тривимірних полів умовної ймовірності ураження обслуговуючого персоналу техногенного об'єкту на основі апарату пробіт-аналізу.
Практична значимість. Розроблена обчислювальна технологія дозволяє експерту на етапі прийняття рішення здійснювати автоматизований чисельний аналіз і прогноз в часі та просторі умовної ймовірності смертельного ураження обслуговуючого персоналу, який піддається інгаляційному впливу токсичного газу як складової частини показника безпеки техногенного об'єкту – індивідуального ризику.

Ключові слова: газові суміші; чисельні методи; диференціальні рівняння зі частинними похідними; вплив шкідливих речовин; забруднення

Y. O. SKOB^1*, M. L. UGRYUMOV^2*

^1*Dep. «Informatics», Kharkiv National Aerospace University named after M. Zhukovsky «Kharkiv Aviation Institute», Chkalov St., 17, Kharkiv, Ukraine, 61070, tel. +38 (057) 315-11-31,
e-mail yuriy.skob@gmail.com, ORCID 0000-0003-3224-1709

^2*Dep. «Informatics», Kharkiv National Aerospace University named after M. Zhukovsky «Kharkiv Aviation Institute», Chkalov St., 17, Kharkiv, Ukraine, 61070, tel. +38 (057) 315-11-31,
e-mail m.ugryumov@khai.edu, ORCID 0000-0003-0902-2735

MATHEMATICAL MODELING OF EVAPORATION

CONSEQUENCES OF TOXIC

SUBSTANCE EMERGENCY SPILLAGE

AT RAILWAY TRANSPORT

Purpose. The main purpose of the article is calculation of spatial distribution fields of the conditional probability of lethal damage to the railway station personnel, caused by the inhalation of toxic gas, which is dissipated in the surface layer of the atmosphere under the conditions of a given wind situation, for a numerical assessment of the safety level of the technogenic object. Methodology. The authors developed a three-dimensional mathematical model of the evaporation processes of toxic chemical substance from the surface of the spillage stain as a result of emergency destruction of the storage or transportation capacity of liquefied gas and further dispersion of the gaseous admixture in the ground layer of the atmosphere, taking into account the cluttering of space by buildings. Also it was developed a calculation technology for determining the conditional probability of human injury by toxic gas on the basis of probit-analysis of the impact degree of damaging factor (inhalation toxodose) on human body. To automate the calculation process, the tabular dependence «probit-function-injury probability» is replaced by a generalized piecewise cubic spline. Findings. Based on the developed model we obtained the results of calculations of the space-time fields of the conditional probability of lethal injury to the railway station personnel who underwent inhalation exposure of a cloud of hydrogen cyanide gas. We also determined that the presence of buildings on the way of the toxic cloud dispersion increases the concentration area and the time of cloud passage along the calculated area, which, accordingly, increases the exposure time of station personnel to harmful impact. Originality. The developed mathematical model takes into account: the flow compressibility, the complex terrain (cluttering of the calculation space by the buildings of technogenic object), the three-dimensional nature of dispersion process of the gaseous admixture, the evaporation of a toxic substance with a variable intensity depending on the wind conditions, physical characteristics of admixture and the roughness grade of atmosphere surface layer. The mathematical model makes it possible to obtain spatio-time distributions of a dangerous parameter – the relative mass concentration of toxic gas and the damaging factor – inhalation toxodose, which are necessary to determine the nonstationary three-dimensional fields of the conditional probability of injury to the technogenic object personnel on the basis of the probit-analysis apparatus. Practical value. The developed calculation technology allows the expert at the decision-making stage to perform automated numerical analysis and forecast in time and space of the conditional probability of lethal injury to service personnel exposed to the inhalation effect of toxic gas as an integral part of the safety index of a technogenic object - individual risk.

Key words: gas mixtures; numerical methods; partial differential equations; exposure to harmful substances; pollution

REFERENCES

1. Belyaev, N. N., & Koptilaya, O. V. (2002). Kompyuternoe modelirovanie zagryazneniya okruzhayushchey sredy pri razlive ammiaka. Ekolohiia i pryrodokorystuvannia, 2, 158-162. (in Russian)

2. Biliaiev, M. M., Kalashnikov, I. V., & Kozachyna, V. A. (2018). Territorial Risk Assessment after Terrorist Act: Express Model. Science and Transport Progress, 1(73), 7-14. doi: 10.15802/ stp2018/123474 (in Russian)

3. Matsak, V. G., & Khotsyanov, L. K. (1956). Gigienicheskoe znachenie skorosti ispareniya i davleniya para toksicheskikh veshchestv, primenyaemykh v proizvodstve. Moscow: Medgiz. (in Russian)

4. RD-03-26-2007 Metodicheskie ukazaniya po otsenke posledstviy avariynykh vybrosov opasnykh veshchestv. (2008). Moscow: NTTs Promyshlennaya bezopasnost. (in Russian)

5. Skob, Y. A. (2017). Matematicheskoe modelirovanie struynogo istecheniya gazovozdushnoy smesi s razlichnoy kontsentratsiey primesi v atmosferu. Aerospace Technic and Technology, 4, 83-92. (in Russian)

6. Stoetsky, V. F., Golinko, V. I., & Dranishnikov, L. V. (2014). Risk assessment in man-caused accidents. Scientific Bulletin of National Mining University, 3, 117-124. (in Russian)

7. Chernyshev, Y. K. (2000). Vypuklye vektornye splayny v primenenii k profilirovaniyu lopatok GTD. Aerospace Technic and Technology, 21, 16-18. (in Russian)

8. Assael, M. J., & Kakosimos, K. E. (2010). Fires, Explosions, and Toxic Gas Dispersions: Effects Calculation and Risk Analysis. New York: CRC Press Publisher. doi: 10.1201/9781439826768 (in English)

9. Biliaiev, M. M., & Muntian, L. Y. (2016). Numerical simulation of toxic chemical dispersion after accident at railway. Science and Transport Progress, 2(62), 7-15. doi: 10.15802/stp2016/67279 (in English)

10. Brauer, R. L. (2016). Safety and Health for Engineers. Hoboken: Wiley Publisher. (in English)

11. Andersson, B., Andersson, R., Hakansson, L., Mortensen, M., Sudiyo, R., & Frontmatter, B. W. (2012). Computational Fluid Dynamics for Engineers. Cambridge: Cambridge University Press Publisher. (in English)

12. Engeln-Müllges, G., Niederdrenk, K., & Wodicka, R. (2010). Numerik-Algorithmen: Verfahren, Beispiele, Anwendungen. Berlin: Xpert.press Publisher. (in Germany)

13. Hughes, P., & Ferrett, E. (2011). Introduction to Health and Safety at Work: The Handbook for the NEBOSH National General Certificate. Kidlington: Oxford, Butterworth-Heinemann. (in English)

14. Knott, G. D. (2012). Interpolating Cubic Splines. Boston: Birkhäuser Publisher. (in English)

15. Nolan, Dennis P. (2011). Handbook of Fire and Explosion Protection Engineering Principles: for Oil, Gas, Chemical and Related Facilities. Burlington: Gulf Professional Publishing, Elsevier. (in English)

16. Dadashzadeh, M., Kashkarov, S., Makarov, D., & Molkov, V. (2018). Risk assessment methodology for onboard hydrogen storage. International Journal of Hydrogen Energy, 43(12), 6462-6475.
    doi: 10.1016/j.ijhydene.2018.01.195 (in English)

17. Nakayama, J., Misono, H., Sakamoto, J., Kasai, N., Shibutani, T., & Miyake, A. (2017). Simulation-based safety investigation of a hydrogen fueling station with an on-site hydrogen production system involving methylcyclohexane. International Journal of Hydrogen Energy, 42(15), 10636-10644. doi: 10.1016/j.ijhydene.2016.11.072 (in English)

18. Skob, Y. A., Granovskiy, E. A., & Ugryumov, M. L. (2017). Mathematical modeling of hydrogen explosion consequences at fueling station. Proceedings of 7-th International Conference on Hydrogen Safety, 1-12. Retrieved from https://www.hysafe.info/wp-content/uploads/2017_papers/159.pdf (in English)

19. Tashvigh, A. A., & Nasernejad, B. (2017) Soft computing method for modeling and optimization of air andwater gap membrane distillation – a genetic programming approach. Desalination and Water Treatment, 76, 30-39. doi: 10.5004/dwt.2017.20696 (in English)

20. Toro, E. F. (2009). Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction. Berlin: Springer Publisher. (in English)

Статья рекомендована к публикации д.т.н., проф. Н. Н. Беляевым (Украина)

Поступила в редколлегию: 06.03.2018

Принята к печати: 07.06.2018.

doi: 10.15802/stp2018/133637 © Ю. А. Скоб, М. Л. Угрюмов, 2018