CC BY
41
Наука та прогрес транспорту. Вюник Дншропетровського нащонального унiверситету залiзничного транспорту, 2018, № 3 (75)
УДК 533.27:519.63:504.054
Ю. А. СКОБ1*, М. Л. УГРЮМОВ
2*
1 Каф. «Информатика», Харьковский национальный аэрокосмический университет им. Н. Е. Жуковского «ХАИ», ул. Чкалова, 17, Харьков, Украина, 61070, тел. +38 (057) 315 11 31, эл. почта yuriy.skob@gmail.com,
ШСГО 0000-0003-3224-1709
2 Каф. «Информатика», Харьковский национальный аэрокосмический университет им. Н. Е. Жуковского «ХАИ», ул. Чкалова, 17, Харьков, Украина, 61070, тел. 38 (057) 315 11 31, эл. почта m.ugryumov@khai.edu,
ОЯСГО 0000-0003-0902-2735
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПОСЛЕДСТВИЙ ИСПАРЕНИЯ АВАРИЙНОГО ПРОЛИВА ТОКСИЧНОГО ВЕЩЕСТВА НА ЖЕЛЕЗНОДОРОЖНОМ ТРАНСПОРТЕ
Цель. Основной целью работы является расчет пространственных полей распределения условной вероятности летального поражения обслуживающего персонала железнодорожной станции, вызванного ингаляцией токсичного газа, который рассеивается в приземном слое атмосферы в условиях заданной ветровой обстановки, для численной оценки уровня безопасности техногенного объекта. Методика. Разработана трехмерная математическая модель процессов испарения токсичного химического вещества с поверхности пятна пролива в результате аварийного разрушения емкости хранения или транспортировки сжиженного газа и дальнейшего рассеивания газовой примеси в приземном слое атмосферы с учетом загромождения пространства зданиями. Также разработана вычислительная технология определения условной вероятности поражения человека токсичным газом на основе пробит-анализа степени воздействия поражающего фактора (ингаляционной токсодозы) на организм человека. Для автоматизации вычислительного процесса табличная зависимость «пробит-функция - вероятность поражения» заменена обобщенным кусочно-кубическим сплайном. Результаты. На основании разработанных моделей получены результаты расчетов пространственно-временных полей условной вероятности смертельного поражения персонала железнодорожной станции, который подвергся ингаляционному воздействию облака газообразного цианистого водорода. Определено, что наличие зданий на пути рассеивания токсичного облака увеличивает площадь концентрации и время прохождения облака по расчетной области, что, соответственно, увеличивает время экспозиции обслуживающего персонала вредному воздействию. Научная новизна. В разработанной математической модели учитываются: сжимаемость потока, сложный рельеф местности (загромождение расчетного пространства зданиями техногенного объекта), трехмерный характер процесса рассеяния газообразной примеси, наличие испарения с пятна пролива токсичного вещества с переменной интенсивностью в зависимости от ветровой обстановки, физических характеристик примеси и класса шероховатости приземного слоя атмосферы. Математическая модель позволяет получать пространственно-временные распределения опасного параметра - относительной массовой концентрации токсичного газа и поражающего фактора - ингаляционной токсодозы, которые необходимы для определения нестационарных трехмерных полей условной вероятности поражения обслуживающего персонала техногенного объекта на основе аппарата пробит-анализа. Практическая значимость. Разработанная вычислительная технология позволяет эксперту на этапе принятия решения осуществлять автоматизированный численный анализ и прогноз во времени и пространстве условной вероятности смертельного поражения обслуживающего персонала, который подвергается ингаляционному воздействию токсичного газа как составной части показателя безопасности техногенного объекта -индивидуального риска.
Ключевые слова: газовые смеси; численные методы; дифференциальные уравнения с частными производными; воздействие вредных веществ; загрязнение
Технологические процессы предприятий авиационной промышленности включают использование, транспортировку и хранение
Введение
отравляющих химических веществ (ОХВ) в сжиженном состоянии [10]. Нарушение правил эксплуатации оборудования приводит к его отказам, которые сопровождаются выбросом в атмосферу ОХВ с образованием токсичных облаков [13]. Одним из наиболее опасных ви-
Наука та прогрес транспорту. Вюник Дншропетровського нащонального унiверситету залiзничного транспорту, 2018, № 3 (75)
дов техногенной аварии является разрушение емкости хранения или транспортировки сжиженного газа (СГ) с образованием пятна пролива [15] (рис. 1). Массовая концентрация токсичного газообразного вещества в газовоздушной смеси характеризует негативное отклонение от нормального химического состава воздуха и наряду с экспозицией является опасным параметром для обслуживающего персонала, оказавшегося в пределах зоны техногенной аварии.
Экспозиция обслуживающего персонала промышленного объекта определенным концентрациям ОХВ формирует поражающий фактор - ингаляционную токсодозу. Превышение пороговых значений токсодозы приводит к социальным последствиям - отравлению различной степени тяжести и человеческим жертвам. Поэтому определение риска предприятия для такого рода аварии является важной и актуальной инженерно-практической задачей.
Аварийное разрушение емкости хранения ОХВ
Выброс жидкой фазы из
Образование пролива ОХВ Испарение жидкости
Рис. 1. Развитие техногенной аварии Fig. 1. Development of technogenic accident
Оценка последствий техногенной аварии включает в себя определение вероятности поражения обслуживающего персонала, который может быть подвергнут воздействию ОХВ, на основе математического моделирования рассеяния токсичной примеси в атмосфере [8].
Математическое моделирование физических процессов выброса и рассеяния вредной примеси в приземном слое атмосферы позволяет
осуществить прогноз полей массовой концентрации ОХВ, определить ингаляционную токсодозу и вероятность поражения обслуживающего персонала [5, 17, 18].
Математическая модель физического процесса
Опасные параметры
поражающего фактора Вероятность поражения
Рис. 2. Схема вероятностной оценки последствий аварии
Fig. 2. Scheme of probabilistic assessment of accident consequences
Цель
Учитывая вышеизложенное, целью данной работы является разработка адекватной математической модели процесса испарения токсичной примеси с поверхности пятна пролива в результате техногенной аварии, формирования токсичного облака и его дальнейшего распространения в приземном слое атмосферы для получения пространственно-временных полей поражающего фактора ингаляционной токсодо-зы и определения полей вероятности поражения человека на основе пробит-анализа.
Анализ публикаций. Адекватное описание физических процессов смешения нереагирую-щих газов с воздухом и дальнейшего распространения смеси при истечении струи в открытое пространство или замкнутое помещение с принудительной (или естественной) вентиляцией возможно только с использованием системы нестационарных уравнений Навье-Стокса для сжимаемого газа. Ограниченные возможности современных компьютеров не позволяют эффективно осуществлять прямое
Наука та прогрес транспорту. Вюник Дншропетровського нащонального ушверситету залiзничного транспорту, 2018, N° 3 (75)
численное решение этих уравнений. В настоящее время численное моделирование турбулентных течений осуществляют путем решения осредненных по Рейнольдсу-Фавру уравнений Навье-Стокса, дополненных моделью турбулентности [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^. Свежий воздух со скоростью ветра поступает через входную грань расчетной
Наука та прогрес транспорту. Вюник Дншропетровського нацюнального ушверситету залiзничного транспорту, 2018, N° 3 (75)
области, перемешивается с примесью, образуя газовоздушное облако с массовой концентрацией Q.
Для упрощения математической модели и ускорения компьютерных вычислений принимаем предположение, что физические процессы перехода выбрасываемого в приземный слой атмосферы вещества (вскипание) из жидкого состояния в газообразную фазу проходят мгновенно и в бесконечно-тонком по высоте слое. Это допущение несколько огрубляет расчеты, но в общем достаточно адекватно позволяет описать процесс попадания газовой примеси в воздух окружающего пространства.
Течение газовой смеси в расчетной области будем определять параметрами окружающей атмосферы, площадью пятна пролива, параметрами газовой примеси, поступающей в результате испарения в атмосферу. В какой-то момент времени испарение может прекратиться, и поступления примеси в область не будет.
Основные уравнения математической модели. Предполагаем, что основным фактором, влияющим на физические процессы смешения газовых примесей с воздухом и дальнейшее распространение смеси при испарении СГ в открытое пространство, является конвективный перенос массы, импульса и энергии. Поэтому достаточно использовать упрощенные уравнения Навье-Стокса, которые получены отбрасыванием вязких членов в уравнениях движения газовой смеси (эйлеров подход с ис-точниковыми членами) [5].
Полная система уравнений, описывающая нестационарное трехмерное течение двухком-понентной смеси газов в данной постановке имеет вид [5]:
da db dc dd -— + — + — + — = pf ,
dt dx dy dz
(1)
где a, b, c, d, f - вектор-столбцы вида:
a = [p, pu, pv, pw, E]T, (2)
b = =ypu, P + pu2, puv, puw, (E + P)uJ , (3) c = |pv, pvu, P + pv2, pvw, (E + P)vJT , (4)
d = pwu, pwv, P + pw2, (Е + P)w] ,(5)
у = [0,0, - g ,0, - ^]Т, (6)
где Т - время; и, V, w - составляющие вектора скорости д ; Р, р - давление и плотность; Е - полная энергия единицы объема газовоздушной смеси:
Е = р(е + 2(и2 + V2 + w2)), (7)
где е - внутренняя энергия единицы массы газа; компоненты вектора; у - суть проекции распределенных объемных источников; g -ускорение свободного падения.
Закон переноса компоненты смеси с учетом скорости диффузии имеет вид [12]:
d(pQ) , d(puQ) , d(pvQ) , d(pwQ)
dt
dx
dy
dz
= Pq , (8)
где Q - относительная массовая плотность примеси (отношение плотности газообразного вещества примеси к плотности смеси), PQ - интенсивность изменения плотности примеси вследствие диффузии в соответствии с законом Фика ра = div(р§DgrаdQ) (коэффициент диффузии ^ определялся по методике, предложенной М. Е. Берляндом [18]).
Система уравнений (1-8) является незамкнутой. Дополним ее уравнениями, определяющими теплофизические свойства компоненты смеси. Для идеального политропного газа величина е связана с Р и р смеси зави-Р
симостью: е = ---— .
(к - 1)р
Граничные условия. На входе будем задавать граничные условия на поверхностях тех граней, которые примыкают к границам расчетной области и через которые в расчетную область поступает атмосферный воздух. Набегающий поток на входе определяется величинами: - полной энтальпии:
^00 = '
к P u2 + v2 + w
к -1 p 2 - функции энтропии:
(9)
Наука та прогрес транспорту. Вюник Дншропетровського нацюнального унiверситету залiзничного транспорту, 2018, № 3 (75)
5 - Р •
(10)
д_ д
Ш айУ + ф А й с — Щр/йУ, (11)
где У - объем элементарной расчетной ячейки; с - ограничивающая поверхность данной ячейки, которая имеет внешнюю нормаль п (с - сп); А - тензор плотности потока консервативных переменных а , столбцами которого являются векторы Ь, с, й соответственно.
Закон переноса компоненты смеси (8) может быть также представлен в интегральной форме для каждой расчетной ячейки:
- направлением вектора скорости потока (углами ах, ау, аг);
- относительной массовой плотностью примеси Q .
Параметры потока на входе определяются из равенств (9, 10) с учетом заданных ах, ау, аг с привлечением соотношения для
«левого» инварианта Римана [19]. На непроницаемых участках, ограничивающих расчетную область поверхностей, выполняются условия «непротекания»: qli — 0 , где п - вектор нормали к рассматриваемой границе. На поверхности испарения выставляется условие протекания примесного газа с заданной интенсивностью. Граничные условия на выходе будем задавать на поверхностях тех граней конечно-разностных ячеек, которые примыкают к границам расчетной области и через которые предполагается вытекание или втекание смеси. В выходных областях, кроме атмосферного давления РА , задаваемого либо взятого из эксперимента, использовались соотношения для «правого» инварианта Римана [20].
Начальные условия. В начальный момент времени во всех «газообразных» ячейках расчетной области принимаются параметры окружающей среды. В ячейках с испарением или истечением газа задается закон изменения расхода примеси.
Алгоритм численного решения. Векторное уравнение (1) является следствием законов сохранения массы, импульса и энергии, которые могут быть представлены в интегральной форме для каждой расчетной ячейки:
д_ а
\WpQdy + $рQqd с — \\\рQdУ . (12)
Метод расчета. Компьютерное решение системы фундаментальных уравнений газовой динамики для смеси, дополненной законами сохранения массы примесей в интегральной форме, получено явным методом С. К. Годунова [19]. Для аппроксимации уравнений Эйлера применяется конечноразностная схема первого порядка. Центральные разности второго порядка используются для диффузионных источни-ковых членов в уравнениях сохранения примесей. Простая интерполяция давления применяется в вертикальном направлении. Метод Поду-нова характеризуется робастным алгоритмом, устойчивым к большим возмущениям параметров потока (например, давления).
В основе метода расчета лежит идея использования для построения разностной схемы точных решений уравнений с кусочно-постоянными начальными данными. Для гиперболической системы такие решения распадаются на совокупности независимых и сравнительно просто рассчитываемых деталей -«распадов разрывов».
Уравнения (11, 12) допускают возникновение и существование поверхностей разрыва двух видов: ударных волн и тангенциальных разрывов. Функции, удовлетворяющие уравнениям (11, 12), можно рассматривать в качестве обобщенных решений уравнений газовой динамики. Использование интегральных законов сохранения массы, импульса, энергии и концентрации газообразной примеси в качестве исходных для построения разностных уравнений обеспечивает построение разрывных решений без выделения разрывов.
Совокупность газодинамических параметров во всех ячейках в момент времени ^ представляет собой известное решение на временном слое с индексом п. Параметры в момент
времени ¿"+1 — ^ + т (на слое п +1) рассчитывались посредством применения явных разностных аппроксимаций для соотношений в рамках интегро-интерполяционного метода
Наука та прогрес транспорту. Вюник Дншропетровського нащонального унiверситету залiзничного транспорту, 2018, № 3 (75)
С. К. Годунова [20]. На первом этапе непрерывное распределение параметров заменяется кусочно-постоянными среднеинтегральными значениями в каждой расчетной ячейке. При этом границы ячейки представляют собой неустойчивые поверхности произвольного разрыва, которые распадаются на устойчивые волновые элементы: ударную волну, контактную поверхность и волну разрежения. Для каждого такого разрыва определяются потоки массы, импульса и энергии через грани газовых ячеек. Устойчивость конечноразностной схемы обеспечивается выбором величины шага по времени т .
Моделирование испарения. В результате дискретизации расчетной области поверхность пятна пролива СГ также разбивается на ряд ко-нечноразностных ячеек у земли в плоскости Х02 (рис. 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
При равномерном разбиении в направлении осей ОХ и 02 площади граней «испаряющих» ячеек одинаковы. Сделав допущение о равномерности потока с пятна пролива, можно определить индивидуальный заданный расход газа для каждой из ячеек «испарения» = GE /к, где к - количество ячеек, примыкающих к пятну пролива сжиженного газа.
Допустим, имеется газовая среда, для которой термодинамические величины - давление Р, плотность р, внутренняя энергия единицы массы в - подчиняются уравнению состояния. Предположим, что в начальный момент времени t для левого полупространства Х<0 среда характеризуется значениями параметров Рь рь иь а для правого полупространства Х>0 - значениями Р2, р2, и2 (здесь и - компонента вектора скорости в направлении координаты 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
Наука та прогрес транспорту. Вюник Дншропетровського нацюнального ушверситету залiзничного транспорту, 2018, № 3 (75)
Если привести в соприкосновение две массы газа, сжатые до различных давлений (P1 - давление со стороны пятна пролива, P2 - давление со стороны атмосферы), и убрать перегородку между ними, то поверхность их соприкосновения будет поверхностью разрыва в начальном распределении давления. Начальный разрыв распадается на несколько разрывов, которые с течением времени будут отходить друг от друга. На контактном разрыве испытывает скачок плотность, а значит, и внутренняя энергия (R1, E1 - для левой и R2, E2 - для правой областей), а давление P и поперечная компонента скорости U непрерывны. В свою очередь, эти области отделены от невозмущенных областей с параметрами (Pi, рь и1) снизу («слева») и (P2, р2, и2) сверху («справа») или ударной волной УВ, или волной разрежения ВР.
Решая задачу распада разрыва на грани ко-нечноразностной ячейки, примыкающей к вентиляционному проему, можно определить плотность R и скорость и, а значит, и индивидуальный расход газа Gj через рассматриваемую грань. Используя метод итераций, можно подобрать давление P2 таким образом, чтобы расчетный расход газа Gj отличался от заданного Gs на наперед заданную малую величину в (рис. 6). Тестирование такого итерационного алгоритма показало быструю сходимость процесса подбора давления «испарения» и незначительное увеличение общего времени нестационарного расчета движения газовой смеси в расчетной области. Так как информация о предыдущем шаге итерации по времени запоминается в специальной структуре данных, то итерационный процесс подбора противодавления в процессе общего расчета ускоряется.
Интерполяция функции интенсивности испарения. При моделировании испарения с пятна пролива интенсивность «выброса» примеси в газовой фазе в атмосферу обычно принимают постоянной G = const (рис. 7).
Если имеется суммарная масса m пролитого СГ и время t1 начала и t2 конца процесса испарения, тогда текущая интенсивность испарения может быть найдена из соотношения:
G = m / (t2 -11) = const.
(13)
Исходные параметры потока (P1, р1, u1) и (P2, р2, u2); заданный расход Gs
Рис. 6. Итерационный алгоритм подбора противодавления в текущий момент времени процесса испарения
Fig. 6. Iterative algorithm of selection of counterpres-sure at the current time of the evaporation process
Рис. 7. Постоянный закон интенсивности испарения:
t1, t2 - время начала и конца процесса испарения;
G.J, - заданная интенсивность испарения
Fig. 7. Constant law of evaporation intensity:
t1, t2 - time of beginning and end of evaporation process;
G.j, -specified evaporation intensity
Иногда закон интенсивности испарения G3 = f (t) (рис. 8) задан табличной функцией, которая получена либо из эксперимента, либо с помощью аналитической модели. В этом случае возникает проблема интерполирования таблично заданной функции. При выборе подходяще-
Наука та прогрес транспорту. Вюник Дншропетровського нащонального ушверситету з&тзничного транспорту, 2018, N° 3 (75)
го метода интерполирования следует учитывать возможный сложный характер функции и неравномерность расположения узлов интерполяции [14].
О., г.'(с*ы2)
Г 1
ч N
А
/
О 1
6 %с
Рис. 8. Интерполирование интенсивности испарения кусочно-кубическим сплайном:
1 - начало испарения; 2 - конец испарения
Fig. 8. Interpolation of the evaporation intensity by piecewise-cubic spline: 1 - evaporation beginning; 2 - evaporation end
Ю. К. Чернышев в работе [7] провел аналитический обзор наиболее известных методов интерполирования (метод конформных отображений, эрмитовы сплайны, сплайн Безье), рассмотрел их преимущества и недостатки и пришел к выводу, что наиболее подходящими являются кусочно-кубические эрмитовы сплайны, в основе которых лежит методика Х. Акимы построения нелинейных приближений первой производной в узлах интерполирования и ее обобщения [12]. Были введены дополнительные промежуточные узлы интерполяции и доказано, что данный прием устраняет колебания знака второй производной.
Площадь пролива F (м2) можно определить по следующей формуле:
F =
M ж - Mг - ми
и
0,05 рж
(14)
Миж - масса опасного вещества, переходящая в аэрозоль в первичное облако, кг.
Скорость испарения с поверхности пролива и расход аммиака во вторичном облаке, образующемся на стадии испарения из пролива, определяем по формуле [3]:
ди = ^ 10-6 (5,38 + 4,1иоиэфф)Рн , кг/с,(15)
где ц - молярная масса опасного вещества, кг/моль; и0иэфф - начальная эффективная скорость вторичного облака, образующегося на стадии испарения из пролива, м/с; Рн - давление насыщенного пара опасного вещества при температуре воздуха, мм рт. ст., которое можно определить по формуле:
(
Рн = 760exp
(
AH„
iP
1
1
Л
T
\ кип
T
(16)
возд у
где рж - плотность жидкого опасного вещества, кг/м3; Мж - суммарная масса пролитого
жидкого опасного вещества, кг/м3; Мг - масса опасного вещества, переходящая в газовую фазу в первичное облако при мгновенном вскипании перегретого опасного вещества, кг;
где Твозд - температура воздуха, К; Ткип - температура кипения жидкого опасного вещества при давлении в окружающей среде Р0 (при нормальных условиях принимается равным 101 325 Па), К; АЯкип - теплота испарения (кипения) жидкого опасного вещества, Дж/кг.
Вероятностная оценка безопасности. В результате моделирования рассеяния токсичной газовой примеси в атмосфере можно получить функции изменения во времени и пространстве массовой концентрации примеси Q . На основании этого определяется опасность воздействия ОХВ на обслуживающий персонал (рис. 2).
Рассмотрим вычислительную технологию определения вероятности поражения. Пусть интегральная степень воздействия t является случайной, удовлетворяющей нормальному закону распределения с математическим ожиданием равным 5 и дисперсией равной 1. В этом случае величина вероятности поражения Р (измеряется в долях единицы) может быть оценена по следующей формуле:
P =
(17)
С другой стороны, интегральную степень воздействия t можно оценить с помощью
уравнения регрессии вида t — а + Ь /(х), где х — {хк} - количественные оценки поражающих факторов; а, Ь - коэффициенты уравнения регрессии.
Задавая верхний предел интеграла (18) (пробит-функцию Рг — ^х), можно определить
вероятность поражения. При воздействии отравляющего химического вещества (ОХВ) на персонал техногенного объекта основным поражающим фактором является ингаляционная токсодоза О - интеграл по времени концентрации ОХВ в воздухе:
тэ
О — | Qnd т, (19)
0
где тэ - время экспозиции (время, за которое набирается ингаляционная токсодоза), с; Q - пространственно-временное значение массовой концентрации ОХВ, ррт; п - табличный коэффициент (например, для цианистого водорода п — 1).
Пробит-функция для поражения человека вследствие воздействия ингаляции смертельной токсодозы отравляющего химического вещества в общем случае определяется по формуле:
Рг — А + В • 1п (О), (20)
где А и В - табличные полуэмпирические коэффициенты (например, для цианистого водорода А — -37,98; В — -3,7).
Обычно в инженерной практике получение вероятности поражения происходит визуально с использованием таблицы зависимости вероятности от пробит-функции поражающего фактора. Это неудобно и не позволяет использовать данный аппарат интегрированным в компьютерной системе для оценки техногенной безопасности. Поэтому для автоматизации процесса связь пробит-функции и вероятности смертельного поражения человека проинтерпо-лируем кусочно-кубическим эрмитовым сплайном, обобщенным Ю. К. Чернышевым (рис. 9).
Результаты
Апробация математической модели. Предложенный алгоритм и метод учета переменной интенсивности испарения сжиженного токсич-
Наука та прогрес транспорту. Вюник Дншропетровського нацюнального ушверситету затзничного транспорту, 2018, № 3 (75)
ного газа с поверхности пятна пролива в процессе выброса газовой примеси в приземный слой атмосферы был реализован в виде подсистемы исследовательского программного комплекса «Fire». Комплекс позволяет производить трехмерный анализ рассеяния токсичных газообразных примесей во времени и пространстве в практически приемлемое время и делать прогноз рисков летального исхода вследствие воздействия испарений токсичного газа на организм человека (рис. 10, 11).
Тестирование разработанной информационной технологии и анализ эффективности алгоритма проводились на примере моделирования испарения 6925 кг сжиженного цианистого водорода (токсичного взрывоопасного вещества плотностью 689 кг/м3, молярной массой 0,027 кг/моль, температурой кипения 298,6 К, теплотой испарения 933 кДж/кг) с пятна пролива в форме окружности радиусом R, образовавшегося в результате разрушения цистерны в районе железнодорожной станции в черте населенного пункта с большими зданиями (коэффициент степенной зависимости для аппроксимации скорости в атмосферном слое над поверхностью земли к = 0,4) (рис. 10).
Р
90
80
70
ВО
50
40
30
20
10
о
2 3 4 5 е 7 Рг
Рис. 9. Интерполирование зависимости вероятности P поражения от пробит-функции Pr поражающего фактора
Fig. 9. Interpolating the dependence of P lesion probability on the probit-function Pr of the adverse factor
Центр окружности пятна пролива располагался на расстоянии Xc = 16 м, Zc = 16 м от начала координат в расчетной области с габа-
Наука та прогрес транспорту. Вюник Дншропетровського нащонального ушверситету затзничного транспорту, 2018, № 3 (75)
ритами ЬххЬухЬг = 85^10^85 м и вариантом по количеству ячеек вдоль координатных осей 85x10x85.
На расстоянии Ха = 30 м и 2а = 28 м от начала координат располагалось здание станции с габаритами ВххВухВг = 15x5x25 м.
Ветер набегает со скоростью 3 м/с под углом 45° к оси 02 на высоте 0,5 м. В этом случае начальная эффективная скорость вторичного облака, образующегося на стадии испарения из пролива, составляла 1,19 м/с.
Если пренебречь массой газовой фракции
Мг первичного облака, образовавшегося в результате аварийного разрушения емкости, и учесть, что температура кипения вещества выше температуры окружающей среды (масса опасного вещества, переходящая в аэрозоль в первичное облако Мж = 0), по формуле (14) можно определить площадь пятна пролива 201 м2 и соответствующий радиус пятна Я = 8 м. Тогда в соответствии с формулой (15) интенсивность испарения цианистого водорода с пятна пролива будет составлять 0,00106 кг/с/м2.
щ W
Z
D
1 j 1
ч
¿а д
Щ *
1 г
2 л
\
Zc. /
/
1 j Л <
ч /
1
> А о
>
i Г1
1
Рис. 10. Карта объектов у земли: 1 - вектор скорости ветра; 2 - пятно пролива; 3 - здание
Fig. 10. Map of objects near the ground: 1 - vector of wind speed; 2 - spillage stain; 3 - building
Считалось, что испарение начиналось с момента времени t1 = 0 с и принудительно прекращалось по истечению t2 = 5 с, например, с помощью заливки пятна пролива специальной пеной. Время окончания расчетов было принято таким, чтобы газовоздушное токсичное облако покинуло пределы расчетной области.
Было выполнено два варианта расчета: 1-й -без учета наличия здания железнодорожной станции и 2-й - с учетом загромождения области расчета зданием (рис. 11, 12).
Анализ изменения полей массовой концентрации цианистого водорода проводился в плоскости Х02 в слое вычислительных ячеек у земли. На рис. 11, а представлена динамика распространения токсичного вещества без учета присутствия здания. Видно, что образовавшееся после испарения облако беспрепятственно покидает расчетную область, постепенно рассеиваясь.
Влияние загромождения расчетной области зданием существенно изменяет картину рассеивания отравляющего химического вещества (см. рис. 11, б). Токсичное облако меняет форму и размеры в соответствии с течением газовой смеси, которое вырабатывается в процессе обтекания здания. Наличие препятствия увеличивает не только максимальные значения массовой концентрации токсичного вещества, но и время прохождения облака по расчетной области, а значит и время экспозиции обслуживающего персонала вредному воздействию.
Нестационарные поля массовой концентрации цианистого водорода можно рассматривать как пространственно-временное распределение опасного параметра исследуемого физического процесса. Эти данные используются для расчета токсодозы как поражающего фактора, про-бит-функции для цианистого водорода и, соответственно, условной вероятности летального исхода для человека при ингаляции данного токсичного вещества (рис. 12).
Используем для сравнения рассмотренных вариантов вычислительного эксперимента характерную площадь Sх расчетной области у земли, которую занимает зона с условной вероятностью летального исхода 50 % и выше.
В случае учета загромождения (рис. 12, б) характерная опасная зона составляла 1290 м2, что существенно больше, чем в первом варианте расчета (рис. 12, а) - 1001 м2.
Наука та прогрес транспорту. Вюник Дншропетровського нацюнального ушверситету залiзничного транспорту, 2018, № 3 (75)
Рис. 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
Качественный и количественный анализ представленных полей условной вероятности летального исхода вследствие поражения человека ингаляционной токсодозой позволяет сделать вывод о существенном влиянии загромождения пространства зданием на расположение и форму опасной для обслуживающего персонала зоны, что обосновывает необходимость учета
Рис. 12. Поле условной вероятности летального исхода человека у земли, %: а - без здания; б - со зданием
Fig. 12. The field of conditional probability of a person's death at the ground, %: a - without building; b - with building
данного фактора при анализе и прогнозе развития техногенной аварии рассмотренного типа на предприятиях с целью выработки рекомендаций по снижению рисков летального исхода.
Научная новизна и практическая значимость
В математической модели учитываются: сжимаемость потока, сложный рельеф местности (загромождение пространства зданиями),
Наука та прогрес транспорту. Вюнпк Дншропетровського нацюнального ушверситету залiзничного транспорту, 2018, № 3 (75)
трехмерный характер процесса рассеяния, наличие испарения с пятна пролива токсичного вещества с переменной интенсивностью. Модель позволяет получать пространственно-временные распределения опасного параметра - относительной массовой концентрации токсичного газа и поражающего фактора - ингаляционной токсодозы, необходимые для автоматизированного определения нестационарных трехмерных полей вероятности поражения обслуживающего персонала на основе пробит-анализа.
Выводы
В работе представлен метод вероятностной оценки летального поражения обслуживающего персонала техногенного объекта, находящегося в зоне аварийного рассеяния облака ток-
сичного химического вещества в условиях застройки. В основу метода положена трехмерная математическая модель процесса испарения с заданной интенсивностью сжиженного токсического газа с пятна пролива, образовавшегося в результате разрушения емкости транспортировки или хранения, и его рассеяния в приземном слое атмосферы. Модель позволяет получать поля ингаляционной токсодозы и с помощью процедуры пробит-анализа численно оценивать условную вероятность поражения персонала, который подвергается воздействию токсичного газа как показателя безопасности техногенного объекта. Дальнейшее совершенствование данного метода лежит в области рассмотрения комбинации аварийных сценариев с различными поражающими факторами.
СПИСОК ИСПОЛЬЗОВАННЫХ ИСТОЧНИКОВ
1. Беляев, Н. Н. Компьютерное моделирование загрязнения окружающей среды при разливе аммиака / Н. Н. Беляев, О. В. Коптилая // Еколопя i природокористування : зб. наук. пр. - Дшпропетровськ, 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. Nie-derdrenk, 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.
Наука та прогрес транспорту. Вюник Дншропетровського нацюнального ушверситету залiзничного транспорту, 2018, № 3 (75)
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. - Iss. 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. - Iss. 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 membranedistilla-tion - 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 Каф. «1нформатика», Харювський нацюнальний аерокосм!чний ушверситет !м. М. £. Жуковського «ХА1», вул. Чкалова, 17, Харк!в, Украша, 61070, тел. +38 (057) 315 11 31, ел. пошта yuriy.skob@gmail.com, ORCID 0000-0003-3224-1709
2 Каф. «1нформатика», Харквський нацюнальний аерокосм!чний ушверситет !м. М. £. Жуковського «ХА1», вул. Чкалова, 17, Харюв, Украша, 61070, тел. +38 (057) 315 11 31, ел. пошта m.ugryumov@khai.edu, ORCID 0000-0003-0902-2735
МАТЕМАТИЧНЕ МОДЕЛЮВАННЯ НАСЛ1ДК1В ВИПАРОВУВАННЯ АВАР1ЙНОГО ПРОЛИТТЯ ТОКСИЧНИХ РЕЧОВИН НА ЗАЛ1ЗНИЧ-НОМУ ТРАНСПОРТ
Мета. Основною метою роботи е розрахунок просторових пол!в розпод!лу умовно! ймов!рност! летального ураження обслуговуючого персоналу зал1знично! станц!!, викликаного !нгаляц!ею токсичного газу, який розсшеться в приземному шар! атмосфери в умовах задано! в!трово! обстановки, для чисельно! оц!нки р!вня безпеки техногенного об'екту. Методика. Розроблено тривим!рну математичну модель процес!в випа-ровування токсично! х!м!чно! речовини з поверхн! плями пролиття в результат! авар!йного руйнування ем-ност! збер!гання або транспортування зр!дженого газу i подальшого розс!ювання газово! дом!шки в приземному шар! атмосфери з урахуванням захаращення простору буд!влями. Також розроблено обчислювальну технолог!ю визначення умовно! ймов!рност! ураження людини токсичним газом на основ! проб!т-анал!зу ступеня впливу вражаючого фактора (!нгаляц!йно! токсодози) на орган!зм людини. Для автоматизац!! обчи-слювального процесу табличну залежн!сть «проб!т-функц!я - ймов!рн!сть ураження» зам!нено узагальненим кусочно-куб!чним сплайном. Результати. На основ! розроблених моделей отриман! результати розрахуншв просторово-часових пол!в умовно! ймов!рност! смертельного ураження персоналу зал!знично! станц!!, який п!ддався !нгаляц!йному впливу хмари газопод!бного ц!ан!стого водню. Визначено, що наявн!сть буд!вель на шляху розс!ювання токсично! хмари зб!льшуе площу концентрац!! та час проходження хмари по розрахун-ков!й област!, що, в!дпов!дно, зб!льшуе час експозиц!! обслуговуючого персоналу шк!дливому впливу. Наукова новизна. У розроблен!й математичн!й модел! враховуються: стислив!сть потоку, складний рельеф м!сцевост! (захаращення розрахункового простору буд!влями техногенного об'екту), тривим!рний характер процесу розс!ювання газопод!бно! дом!шки, наявн!сть випаровування з плями пролиття токсично! речовини з! зм!нною !нтенсивн!стю в залежност! в!д в!трово! обстановки, ф!зичних характеристик дом!шки ! класу шорсткост! приземного шару атмосфери. Математична модель дозволяе отримувати просторово-часов! роз-под!ли небезпечного параметра - в!дносно! масово! концентрац!! токсичного газу ! вражаючого фактора -!нгаляц!йно! токсодози, як! необх!дн! для визначення нестац!онарних тривим!рних пол!в умовно! ймов!рнос-т! ураження обслуговуючого персоналу техногенного об'екту на основ! апарату проб!т-анал!зу.
Наука та прогрес транспорту. Вюник Дншропетровського нацюнального ушверситету залiзничного транспорту, 2018, № 3 (75)
Практична значимкть. Розроблена обчислювальна технологiя дозволяе експерту на етат прийняття рь шення здшснювати автоматизований чисельний аналiз i прогноз в чаа та просторi умовно! ймовiрностi смертельного ураження обслуговуючого персоналу, який пiддаеться iнгаляцiйному впливу токсичного газу як складово! частини показника безпеки техногенного об'екту - iндивiдуального ризику.
Ключовi слова: газовi сумшц чисельнi методи; диференцiальнi рiвняння зi частинними похвдними; вплив шшдливих речовин; забруднення
Y. O. SKOB1*, M. L. UGRYUMOV2*
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
Наука та прогрес транспорту. Вюник Дншропетровського нащонального ушверситету затзничного транспорту, 2018, № 3 (75)
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 razlich-noy 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-Mullges, 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: Birkhauser 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 and-water 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.
