УДК 519.6 : 504.3.054
Н. Н. БЕЛЯЕВ (ДИИТ)
РАСЧЕТ ЗАГРЯЗНЕНИЯ АТМОСФЕРЫ ПРИ АВАРИЯХ С ХИМИЧЕСКИ ОПАСНЫМИ ВЕЩЕСТВАМИ НА БАЗЕ КОДА AIR-3
Розглянуто принципи побудови пакета прикладних програм AIR-3 чисельного моделювання xiMi4Horo забруднення атмосфери у випадку аварш. Моделювання базусться на чисельному штегруванш рiвняння конвективно-дифузiйного переносу домiшки. Наводяться результати обчислювального експерименту.
Рассмотрены принципы построения пакета прикладных программ AIR-3 численного моделирования химического загрязнения атмосферы в случае аварий. Моделирование основывается на численном интегрировании уравнения конвективно-диффузионного переноса примеси. Приводятся результаты вычислительного эксперимента.
The basic principles of development of the software AIR-3 aimed at simulation of the atmosphere chemical pollution after accidents are under consideration. The modeling is based on the numerical integration of the ^-gradient transport model. The results of numerical experiment are presented.
Введение
Прогноз химического загрязнения атмосферы при авариях на железной дороге - важная задача в области экологической безопасности. Целью такого прогноза является определение интенсивности и размеров зон химического загрязнения, выявление подзон токсичного поражения различной тяжести, зон, где может быть взрыв, возгорание и близость этих зон к жилой застройке, и другим объектам. Данная информация служит основой для разработки адекватных мер защиты окружающей среды, выбора маршрутов эвакуации, количества транспортных средств, необходимых для эвакуации, обоснования сроков эвакуации и т.д. Основой получения прогнозных данных является расчет уровня загрязнения окружающей среды, и в первую очередь атмосферы, при аварии. Нормативная методика [7], используемая в настоящее время для прогноза последствий аварий с химически опасными веществами на транспорте, обладает рядом существенных недостатков:
1. В случае аварии с перегретыми жидкостями, загрязнение атмосферы происходит в два этапа - это загрязнение от первичного облака, а потом от шлейфа испаряющейся от зоны разлива жидкости. Учет такого развития событий нормативная методика не предусматривает.
2. При оперативном прогнозе данная методика не учитывает направление ветра, и прогноз выполняется только скорости ветра V = 1 м/с.
3. Согласно нормативной методике принимается, что в случае аварии население понесет такие потери: до 25 % населения -легкое поражение; до 40 % - поражение
кое поражение; до 40 % - поражение средней тяжести; до 35 % - смертельное поражение (см. [7], с. 7). Т.е., величина потерь просто назначается, а не обосновывается расчетом.
4. В данной методике принято, что при свободном разливе на грунт испарение жидкости осуществляется со слоя жидкости на грунте высотой до 0,05 м. Таким образом, данная методики не учитывается процесс испарения жидкости из грунта, куда жидкость фильтруется с поверхности земли.
5. Для определения глубины загрязнения территории используют данные таблиц, однако искомые величины в таблицах приводится только для ограниченного диапазона количества жидкости (от 0,5 до 300 т), хранящегося в резервуарах. Поэтому, например, в случае разливов менее 0,5 т методика не позволяет выполнить прогноз.
6. В применяемой методике прогноз относится к зоне до 10 м от поверхности земли. Таким образом, данная методика не дает возможности оценить уровень загрязнения, в частности на более высоких уровнях от поверхности земли с целью оценки возможности затекания загрязненного воздуха в систему вентиляций в здании.
7. Нормативная методика не определяет изменение величины концентрации опасного вещества в атмосфере с течением времени, а значит, не дает картину о динамике распространения поражающего фактора, что в свою очередь, не дает возможность научно-обоснованно прогнозировать риск поражения людей.
© Беляев Н. Н., 2010
159
8. Как известно, скорость и направление ветра могут меняться. Если это произошло, то нормативная методика не может «уточнить» прогноз для новой метеорологической ситуации, потому, что уже произошло частичное испарение разлитой жидкости, но какой части -неизвестно, и поэтому выполнить прогноз для нового значения скорости и направления ветра с помощью данной методики невозможно.
Поэтому значительный интерес представляет создание математических моделей, позволяющих оперативно получать прогнозные данные и не имеющих перечисленные недостатки. В этой связи важным вопросом является разработка численных моделей (СББ модели) [1, 2, 4, 5, 8 - 10], которые позволяют существенно расширить возможности исследователей по получению прогнозной информации, необходимой для оценки последствий аварийных ситуаций. Для экспресс-прогноза, на этапе экспертной оценки, можно ограничится кинематической моделью переноса примеси (модель градиентного типа) для расчета динамики загрязнения атмосферы без привлечения уравнений гидродинамики, применяемых для моделирования переноса примеси с учетом рельефа, зданий, лесных массивов и т.д.
В работе рассматривается описание трехмерной численной модели и применение созданного на ее основе кода ЛШ-3 для решения задачи прогноза загрязнения атмосферы при аварийной утечке химически опасных веществ на различных участках Приднепровской железнодорожной дороги вблизи г. Днепропетровска.
Математическая модель
Моделирование процесса переноса химически опасного вещества в атмосфере в случае аварийного выброса (разлива) осуществляется на базе уравнения переноса примеси (модель градиентного типа) [5]
дС диС ЗуС д(ш - )С _
дt дх
ду
дг
д . дС. д . дС. д . дС. _ — (цх-) +—(К —) +—(цъ-) +
^г х ^ ' ^г у л / л \г ъ '
дх дх ду ду дг дг
(Щг - Г), (1)
где С - концентрация примеси, попадающей в атмосферу при аварии; и, V, ш - компоненты вектора скорости воздушной среды; - скорость оседания примеси (в данной работе принята равной нулю); ц = (ц , цу , - коэффициент турбулентной диффузии; Q - интенсив-
ность выброса токсичного вещества от зоны разлива; 5 (г - г,) - дельта-функция Дирака; г, =
(х, , у,, *) - координаты источника выброса.
Постановка краевых условий для данного уравнения рассмотрена в работах [5, 4]. Если рассматривается задача о миграции облака, образовавшегося на месте аварии, то в начальный момент времени задается форма облака и концентрация загрязнителя в нем.
В построенной численной модели используется неравномерный профиль скорости ветра и коэффициентов диффузии:
( г >
и _ и
\ г
; _ 0,11 *;
ц, = ко • и ; Цу = ц ,
где и\ - скорость ветра на высоте * п = 0,15; к0 - параметр [3].
= 10 м;
Метод решения
Рассмотрим построение разностной схемы для численного интегрирования уравнения переноса примеси [4, 1]. Выполним следующую аппроксимацию производных, входящих в уравнение переноса примеси. Заменим производную по времени разделённой разностью «назад»:
дС сп+ - с,
ф
ы
ы
Конвективные производные представим в виде:
диС ди+С ди - С
дх дх ь-: дх
дvС дv+С дv-С
ду ду ду
дшС дш+С дш~С ь-
дг
и + и
дг
и - и
дг
v + M
где и _■
•; и =■
V _ -
V _
v-|v| + ш +
w - ш
;ш _-
2 2 2 Для аппроксимации конвективных производных используем выражения:
ди + С ~ Ф• к Сф - и^к ^^фк _ Сп+1 .
дх
Ых
2
2
2
du С ui+i, j, k C;'+i! j ,k uijk Ci~jk _ L_ С n+i .
д x A x
dv+C ^ Vi, j+1,kCijk _ VijkCi, j_1,k
ду Ay
dv_C ^ Vi, j+1,kCi, j+1,k _ VijkCijk
ду Ay
dw+C ^ w+j,k+1Cjk _ wijkCi, j,k_1 _
dz Az
dw~C , wi,j,k+1Ci, j,k+1 _ w_kCi, j,k _
x '
dz
Az
dx l"x dx]""* Ax2
_ MxCn+1 + M+xCn+1;
/ rn+1 _ rn+1 rn+1 _ rn+1
d|.. dC 4, j+1,k 4jk _" 4, jkk 4,j-u _
" "Ay2
1"+ ф 1*"+ A_"y
Aj2
_ M'_yCn+x + M+C1
„ / „ ~ \ ,^n+1 ,^n+1 ,^n+1 fin+1
д I dC I _ Ci, i,k+1 Cijk ~ Ci, j,k Cij,k_1
dz Гz dz
i, j,k+1 ijk i, j,k ij,k_1 _ "z ~z Az2 "z Az2
_ M-Cn+1 + M+Cn+1.
At
. T_ f-1 n+1 . T+ n+1 . T_ n+1 . _n+1
+LyС + Lz С + Lz С + H Сг]к _
_ (M+ Сn+1 + M~ Сn+1 + M+ Сn+1 +
\ xx xx yy
С n+k _Сп i
С + 1 Т+Ск +L+ Ск +L+Ск) + н Скк _ At 2 x * z ' 4 ij
_ -4 (M+x Сk + М_Сn + Mjy Сk +
+ MJ+ Сn+M+ Сk + M_ Сn ) ; (3)
7 1 1
- на втором шаге k_n + 2 c_n + 4:
Сk _ Сc 1
Сijk 4jk 1
At
+ 2 ( Сk + L-уСk + L_ Сk ) Ск}
Компоненты вектора скорости будут рассчитываться, при решении гидродинамической задачи, на гранях разностной ячейки (контрольного объема). Вторые производные аппроксимируем следующими выражениями:
д ( дС | - /+1,к ук - - г-1,_
~ Ах2
В данных выражениях Ьх, Ьх, Ьу, Ь- , 1+, Ь-, М.+,, Мух и т.д. - обозначения разностных операторов. С учетом этих обозначений разностный аналог трехмерного уравнения переноса примеси будет иметь вид:
С п +1 С п
ик - ик у Ь Спу1 + П Спу1 + Ь Спу1 +
+M+0n+1 + M+ССn+1 + M_Сn+1) .
Расщепим решение данного разностного уравнения при интегрировании на временном интервале & так:
к 1
- на первом шаге к _ 4:
_ 1 (M-JZk + M+Сc + MJy Сk +
+ MJ+Сc + M_z Сk +M+ Сс); (4)
J , 3 , 1
- на третьем шаге k _ n+4; c _ n + ^ используется формула (4);
- на четвертом шаге k _ n + 1; c _ n + 4 используется формула (3)
На пятом шаге расщепления учитывается влияние источников на изменение концентрации примеси, и расчетная зависимость на этом шаге имеет вид:
5 n+1 5 n
Сj,j,k_Сi,j,k _fql (tn+12) 5 .
A t f_1 Ax Ay Az 1
В дискретном виде дельта-функция Дирака «размазывается» по объему разностной ячейки с учётом сбережения суммарного количества выбрасываемого загрязнения. Функции 51 тождественно равняются нулю, кроме ячеек, где расположен 1-й источник загрязнения. На каждом шаге расщепления расчет неизвестной концентрации осуществляется по явной формуле бегущего счета.
На базе рассмотренной численной модели создан код AIR-3, реализованный на алгоритмическом языке FORTRAN.
Дополнительные расчетные процедуры
Одним из наиболее частых следствий аварий является разлив опасных веществ. В разработанном коде для расчета интенсивности испарения от зоны аварийного разлива используется несколько подходов. При экспресс прогнозе применяется зависимость [6]
Q _ (5,83 + 4,1 • V)Pn4M , (5)
где Q - интенсивность испарения;
V - скорость ветра;
Рн - давление насыщенных паров;
М - молекулярная масса вещества.
Для более детального расчета процесса испарения при разливе на грунт и фильтрации загрязнителя внутрь грунта (с учетом свойств грунта) рассматривается решение сопряженной задачи испарения на границе «поверхность грунта - атмосфера» [1].
В разработанном коде есть возможность моделировать загрязнение от движущегося источника эмиссии токсичного вещества, например, движущейся дымящей цистерны. В этом случае параметр интенсивности эмиссии Q зависит как от времени, так и от координат х, у. В качестве исходной информации также вводятся данные о скорости движения источника эмиссии и его маршрут.
Для оценки риска токсичного поражения людей в коде используется сопоставление величины расчетного значения концентрации со смертельной концентрацией в точке расположения рецептора или расчет токсодозы для людей, попавших в зону химического загрязнения и выходящих из этой зоны согласно маршруту эвакуации [1]:
TD = j C(х, y, z, t) dt,
где ^ - время экспозиции; С - концентрация токсичного вещества в точке расположения рецептора.
В разработанном коде для задания маршрута эвакуации применяются маркеры. В коде предусмотрено задание не только маршрута эвакуации, но и скорости передвижения на нем людей.
Кроме данной расчетной зависимости в коде А1Я-3 для прогноза риска токсичного поражения используется зависимость для расчета «пробит»-функции. Эта функция представляет собой комбинацию ряда эмпирических констант и величины токсодозы и имеет вид:
P = а + р + Tdd .
где
Тёё = |СV (х, у, г, ¿) dt,
о
а, в, V - константы, зависящие от вида токсичного вещества и человека (например, ребенок, пожилой, мужчина).
После расчета поля загрязнения на момент
времени t после аварии определяется величина токсодозы (пробит-функции) в расчетной области, отдельно рассчитывается зона, где значение токсодозы равно или превышает ее смертельное значение для данного токсичного вещества и строится поле этой зоны с целью визуальной оценки ее размеров относительно жилой застройки.
В построенной численной модели используется метод маркирования [7] («porosity technique») для выделения разностных ячеек, которые относятся к области аварийного выброса (разлива) химически опасного вещества и области жилой застройки.
Практическая реализация модели
Разработанный код был применен для проведения серии вычислительных экспериментов по оценке интенсивности химического загрязнения атмосферы и подстилающей поверхности в случае возможных аварий при транспортировке опасных грузов на участке Приднепровской железной дороги: станция Днепропетровск-Южный - станция Днепропетровск-Грузовой.
В результате расчетов определялось время подхода шлейфа токсичного вещества (трифто-рид бора, винилхлорид, аммиак) к жилой застройке для различных метеоситуаций и размеров эмиссии, время формирования зон с поражающими концентрациями, возможные потери среди населения. В качестве иллюстрации результатов расчетов на рис.1, 2 представлена зона загрязнения атмосферы в случае гипотетической аварии на ст. Днепропетровск-Грузовой (выброс винилхлорида, скорость ветра V10 = = 7 м/c), а на рис. 3-5 - динамика развития зоны загрязнения в случае гипотетической аварии на ст. Встречная (выброс трифторида бора, скорость ветра V10 = 4,5 м/c). Из представленных рисунков хорошо видно, что в случае аварии в течение короткого промежутка времени шлейф токсичного вещества достаточно быстро достигнет жилой застройки, и создастся угроза токсичного поражения людей. Проведенные вычислительные эксперименты показали, что за столь короткий промежуток времени не удастся не только осуществить эвакуацию людей из зоны поражения, но и оповестить население об угрозе поражения.
Отметим, что для решения каждой задачи потребовалось минимальное время - порядка 5 с, что является немаловажным фактором при решении прогнозных задач данного класса.
Рис. 1. Зона загрязнения атмосферы в случае гипотетической аварии на ст. Днепропетровск-Грузовой, выброс винилхлорида (сечение г = 2 м), г = 25 c
Рис. 2. Зона загрязнения атмосферы в случае гипотетической аварии на ст. Днепропетровск-Грузовой, выброс винилхлорида (сечение г = 2 м), г = 34 с
Рис. 3. Зона загрязнения атмосферы в случае гипотетической аварии на ст. Встречная, выброс трифторида бора, г = 18 с (сечение г = 2 м)
Рис. 4. Зона загрязнения атмосферы в случае гипотетической аварии на ст. Встречная (выброс трифторида бора), г = 60 с (сечение г = 2 м)
Рис. 5. Зона загрязнения атмосферы в случае гипотетической аварии на ст. Встречная, (выброс трифторида бора, г = 160 с (сечение г = 2 м)
Выводы
Рассмотрены основные принципиальные особенности разработанного кода AIR-3, предназначенного для прогноза последствий аварий с химически опасными веществами. Прогноз загрязнения атмосферы осуществляется на базе численного интегрирования трехмерного уравнения переноса примеси. Рассмотренная численная модель может быть отнесена к группе «screening models», позволяющих быстро получать прогнозную информацию при определенной ограниченности исходных данных, характерной для задач рассматриваемого класса.
Выполненные с помощью кода AIR-3 расчеты загрязнения атмосферы при гипотетическом аварийном выбросе токсичных веществ на ряде конкретных участков железной дороги показали, что в случае аварийной ситуации за короткий промежуток времени произойдет опасное
загрязнение воздушной среды в селитебной зоне, расположенной вблизи железой дороги, по которой осуществляется транспортировка этих грузов.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Антошкина, Л. И. Оценка экологического риска при авариях с химически опасными веществами [Текст] / Л. И. Антошкина, Н. Н. Беляев, Е. Ю. Гунько. - Д.: Наука и образование, 2008. - 132 с.
2. Беляев, Н. Н. Численное моделирования загрязнения воздушной среды на промплощад-ках [Текст] / Н. Н. Беляев, П. Б. Машихина, Е. Ю. Гунько // Вюник Дншропетр. нац. ун-ту залiзн. трансп. 1м. акад. В. Лазаряна. - 2007. -Вип. 16. - Д.: Вид-во ДНУЗТ, 2007. - С. 18-20.
3. Берлянд, М. Е. Прогноз и регулирование загрязнения атмосферы [Текст] / М. Е. Берлянд. - Л.: Гидрометеоиздат, 1985. - 273 с.
4. Численное моделирование распространения загрязнения в окружающей среде [Текст] / М. З. Згуровский [и др.]. - К.: Наук. думка, 1997. - 368 с.
5. Марчук, Г. И. Математическое моделирование в проблеме окружающей среды [Текст] / Г. И. Марчук. - М.: Наука, 1982. - 320 с.
6. Мацак, В. Г. Гигиеническое значение скорости испарения и давления пара токсических ве-
ществ, применяемых в производстве [Текст] / В. Г. Мацак, Л. К. Хоцянов. - М.: Медгиз, 1959. - 231 с.
7. Методика прогнозування насладив впливу (ви-киду) небезпечних х1м1чних речовин при аварь ях на промислових об'ектах i транспорт! [Текст]. - К., 2001. - 33 с.
8. Mahfound, K. Computation of Wind Flow and Air Pollution for Regions Having a Complex Topography [Text] / Kadja Mahfound, Anagnostopoulos [et al.] // Proc. of 3rd European & African Conf. on Wind Engineering (Eindhoven Univ. of Technology, Netherlands, July 2-6, 2001). - P. 355-358.
9. Maurizi, A. Numerical Simulation of the Atmospheric Flow in a Mountainous Region of North Portugal [Text] / A. Maurizi, J. M. L. M. Palma, F. A. Castro // Proc. of 2nd EACWE (Genova, 1997). - Vol. 1. - P. 293-300.
10. Belyaev, M. Air Pollution Modelling of Tech-nogenic Catastrophes Caused by Terrorist Attacks on Rail Transport of Chemical Substances [Text] // Simulation and Assessment of Chemical Processes in a Multiphase Environment. NATO Science for Peace and Security Series. - C.: Environmental Security, Springer, 2007. - P. 327-336.
Поступила в редколлегию 02.06.2010.
Принята к печати 14.06.2010.