Расчет аварийного рассеивания токсичного газа в атмосфере с учетом рельефа местности на базе кода terrain-2i Текст научной статьи по специальности «Математика»

УДК 519.6 : 504.3.054

Н. Н. БЕЛЯЕВ, П. Б. МАШИХИНА (ДИИТ)

РАСЧЕТ АВАРИЙНОГО РАССЕИВАНИЯ ТОКСИЧНОГО ГАЗА В АТМОСФЕРЕ С УЧЕТОМ РЕЛЬЕФА МЕСТНОСТИ НА БАЗЕ КОДА TERRAIN-2I

Розглянуто 2D чисельна модель для прогнозування забруднення атмосфери з урахуванням рельефу. Наведено результати розрахуншв на баз1 розроблено! модель

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

Рассмотрена 2D численная модель для прогноза уровня загрязнения атмосферы с учетом рельефа местности. Представлены результаты численного моделирования.

Ключевые слова: рассеивание токсичного газа, прогнозирование загрязнения атмосферы, учет рельефа местности, двумерная численная модель

The 2D numerical model to simulate the pollutant dispersion over complex terrain was proposed. The results of numerical experiment are presented.

Keywords: dispersion of toxic gas, forecasting of air pollution, taking into account the terrain, two-dimensional numerical model

Введение

Особенностью аварий на железной дороге является то, что во многих случаях эти аварии происходят в условиях сложного рельефа - насыпи, выемки и т.д. (рис. 1). Поэтому при прогнозе последствий аварий с химически опасными веществами весьма важно учесть влияние рельефа на процесс переноса загрязнителя, т.к. рельеф местности способен кардинально изменить процесс рассеивания. Как известно, влияние рельефа местности приводит, в первую очередь, к формированию неравномерного воздушного потока. В этом случае нельзя найти аналитическое решение, как гидродинамической задачи, так и задачи транспорта загрязнителя [1, 7, 8]. Поэтому применение на практике аналитических моделей для прогноза загрязнения атмосферы с учетом реального рельефа -исключено. Нормативная (эмпирическая) методика прогноза последствий аварий с химически опасными веществами не учитывает не только влияние рельефа на процесс переноса загрязнителя, но также и изменение скорости ветра с высотой, возможность выброса переменной интенсивности и т.д. [3, 6]. Единственным теоретическим инструментом решения прогнозных задач в данном случае является применение численных моделей (СББ модели) [2, 4, 5, 10 - 12]. Необходимо отметить, что при проведении расчетов по оценке возможных последствий аварий характерной особенностью является требование минимальных затрат ком-

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

Математическая модель

Процесс расчета переноса загрязняющих веществ в атмосфере разбивается на два этапа. На первом этапе решается задача гидродинамики и определяется поле скорости воздушного потока с учетом его деформации при обтекании рельефа. Для учета отрыва потока будем использовать модель отрывных вихревых течений идеальной несжимаемой жидкости [2]. В этом случае для расчета поля скорости воздушного пока необходимо найти решение нелинейной системы уравнений гидродинамики [2]:

© Беляев Н. Н., Машихина П. Б., 2010

уравнения переноса завихренности

д v д и

ю =---в воздушном потоке:

д х д у

дю ди ю дv ю

д t д x д y

= 0;

(1)

уравнения Пуассона для расчета функции тока у:

д2у д2 у

дХ2 ду2

= -ю .

(2)

Значение компонент скорости определяется по следующим зависимостям:

и =

ду

д у'

v = -

ду

д x

При использовании рассматриваемой модели будем считать, что отрыв потока происходит в угловых точках рельефа. При отрыве потока образуются вихри, интенсивность которых определяется из суммарной циркуляции скорости по контрольному контуру и физических допущений о «нулевой» относительной толщине пограничного слоя при реальном обтекании. Алгоритм расчета интенсивности вихрей рассмотрен в работе [2].

Решение данных уравнений гидродинамики позволяет рассчитать поле скорости воздушного потока при обтекании рельефа. На втором этапе осуществляется решение уравнения переноса загрязнителя в атмосфере [1, 4]:

дС диС DvC — +-+-

дt

N

сС = div (ц gradC) -

дх ду Xa (t)5(Х - х )(у - у ),

г=1

где С - концентрация загрязнителя; и, V -компоненты вектора скорости воздушной среды; |о, = (ц,х, цу ) - коэффициент турбулентной диффузии; Q - интенсивность выброса загрязнителя; 5(х - )5(у - у1) - дельта-функция Дирака; xi, yi - координаты источника выброса; с - коэффициент, учитывающий химический распад загрязнителя, вымывание осадками; ^ - время.

Постановка краевых условий для данного уравнения рассмотрена в работах [2, 4].

Для рассмотренных выше уравнений гидродинамики осуществляется постановка следую-

щих граничных условий. Так, для уравнения Пуассона на поверхности препятствий (рельеф) ставится граничное условие «непротекания» [2]:

у = 0.

На границе входа воздушного потока в расчетную область задается профиль скорости воздушного потока и = f(у) и соответствующее значения функции тока у и завихренности ю:

у1 вход = у(у); ю I вход = ю(у).

На верхней границе расчетной области также ставится условие «непротекания»:

у = const.

На границе выхода воздушного потока из расчетной области ставятся «мягкие» граничные условия, необходимые для замыкания разностных уравнений [2].

Решение стационарной гидродинамической задачи находится методом установления решения по времени t. Для этого в начальный момент t = 0 для завихренности ставится начальное условие типа ю 11=0 = 0 или ю 11=0 = ю0(х, у).

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

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

Рис. 1. Авария на железной дороге (сход состава с насыпи)

Применение метода маркирования (porosity technique) дает возможность пользователю формировать любую форму рельефа, не внося

изменении в расчетный код, делая изменения только в файле исходных данных.

Метод решения

Для численного интегрирования уравнения переноса завихренности применяется попеременно-треугольная неявная разностная схема [2]. Разностные уравнения для расчета завихренности, на каждом шаге расщепления, имеют вид:

- на первом шаге расщепления

i

ю

i

ю

1

1 1 n+— n+—

i Ю,i2 - <1 Ю-w

At

"h+1^ г, i

Ч 1Ш1, i-i

Ay

Ax

= 0;

- на втором шаге расщепления

и+1 п+1

Ю j j

At

u+1,jюП+1 j - u-j

Ax

vl,j = 0,

u + u

где u =■

2

Ay

ит. д.

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

.п + 1

ном слое ^ .

Для расчета функции тока применяется разностная схема суммарной аппроксимации [5]. Разностные уравнения на каждом шаге расщепления имеют вид: - на первом шаге

1

n+—

Vij 4 -V;

У

ю

У .

At

- на втором шаге

1

n +—

Vi,j2 -V

i

1 1 n +— n+—

Vi,j2 -V*-W

1 1

n+— n +—

Vi,j2 -Vi,j-1 .

.2

At Ax

- на третьем шаге

3 13 3

n+— n+— n+— n+—

Vij4 -Vu2 %+14, -V,4

Ay

2

ri. i

i

n+—

Vi.j +41

Vi,

At Ax

- на четвертом шаге

2

Ay2

n+1

Vij -V At

3

n+—

4

Ю-

где ю1=-4 (ю

i

"Юг-1, 1+1Юг-1,1 -1 + Юг, j

В угловых точках, где допускается отрыв потока, располагается источник завихренности. Алгоритм расчета интенсивности завихренности в угловой точке рассмотрен в работе [2]. В дискретной модели это реализуется с помощью маркеров.

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

Практическая реализация модели

На основе рассмотренной численной модели разработан код TERRAIN-2I на алгоритмическом языке FORTRAN. Данный код был применен для решения следующей задачи - затекание загрязненного воздуха в каверну, расположенную на подветренной стороне насыпи (рис. 2) в случае аварийного выброса токсичного газа. Полагалось, что на месте аварии возле насыпи произошел выброс аммиака, и образовалось облако сложной формы: максимальная высота «гриба» - 12 м; максимальная ширина «шляпки» - 80 м; ширина «ножки» - 20 м. При моделировании учитывался неравномерный профиль ветра на входе в расчетную область и изменение вертикального коэффициента диффузии с высотой:

( Z >

u = u

Z1

ц = 0,11 г; где 21 = 10 м, п = 0,15.

Моделирование рассеивания загрязнителя выполнено при таких параметрах: размеры расчетной области 250 м х 21 м; У10 = 6,1 м/с; ц = = 4 м2/с. Максимальная высота «гриба» - 12 м; максимальная ширина «шляпки» - 80 м; ширина «ножки» - 20 м; средняя концентрация токсичного газа в облаке 13 г/м3. Внутри насыпи располагается каверна (выемка), имеющая

2

2

сложную «Г-образную» форму. Верхняя часть этой выемки отмечена буквой «В» (рис. 2).

Рассмотрим результаты вычислительного эксперимента. На рис. 2-4 показана динамика загрязнения атмосферы для различных моментов времени.

Рис. 2. Зона загрязнения атмосферы, г = 0,01 с

Рис. 3. Зона загрязнения атмосферы, г = 35 с

Рис. 4. Зона загрязнения атмосферы, г = 65 с

Видно, что к моменту времени г = 35 с уже большая часть каверны заполнена токсичным газом. К моменту времени г = 65 с загрязненный воздух заполнил всю каверну. К этому моменту времени дифракция облака на насыпи практически закончилась и наветренная сторона насыпи находится вне зоны загрязнения.

Для оценки риска поражения людей, которые спрятались в выемке, был выполнен расчет величины токсодозы

где г - время экспозиции; С - концентрация токсичного вещества в точке расположения рецептора.

В табл. 1 представлены данные о динамике изменения величины токсодозы в выемке (точка х = 145 м; у = 3,5 м). В табл. 2 представлены данные о динамике изменения концентрации токсичного газа у подножия насыпи, с её подветренной стороны (точка х = 205 м; у = 3,5 м). Это место также можно рассматривать как естественное укрытие.

Таблица 1 Динамика изменения величины токсодозы

г, с 44 55 64

ТО, 0,043 0,139 0,275

г/м3 • мин

Таблица 2

Динамика изменения величины концентрации токсичного газа за холмом

г, с 33 44 55

С, г/м3 3,45 2,76 1,30

ТВ = | С(х, у, г, г) Л,

Анализ данных, представленных в таблицах, показывает, что значение токсодозы, при которой могут появиться начальные симптомы поражения (ТВ = 0,25 г/м3 • мин), будет достигнуто к моменту времени г = 64 с. За холмом будет формироваться опасная зона загрязнения, с риском тяжелого токсичного поражения людей (ПДК для аммиака, в рабочей зоне, составляет 20 мг/м3).

Результаты данной задачи иллюстрируют возможность применения построенной численной модели для решения задач с крайне сложным рельефом местности, который включает в себя комбинацию таких элементов как «холм», «впадина», «выемка».

На расчет задачи потребовалось 10 с компьютерного времени.

Выводы

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

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Бруяцкий, Е. В. Теория атмосферной диффузии радиоактивных выбросов [Текст] / Е. В. Бруяцкий. - К. Ин-т гидромеханики НАН Украины, 2000. - 443 с.

2. Численное моделирование распространения загрязнения в окружающей среде [Текст] / М. З. Згуровский [и др.]. - К.: Наук. думка, 1997. - 368 с.

3. Методика прогнозування насладив виливу (ви-киду) небезпечних х1м1чних речовин при аварь ях на промислових об'ектах i транспорт! [Текст]. - К., 2001. - 33 с.

4. Марчук, Г. И. Математическое моделирование в проблеме окружающей среды [Текст] / Г. И. Марчук. - М.: Наука, 1982. - 316 с.

5. Самарский, А. А. Теория разностных схем [Текст] / А. А. Самарский. - М.: Наука, 1983. -616 с.

6. Хмшь, Г. А. Концептуально-методичний апарат анал1зу й оцшки техногенних та природних ри-зишв [Текст] / Г. А. Хмшь // Еколопя довкшля та безпека життедгяльностг - 2007. - № 5. -C. 47-55.

7. Шаталов, А. А. Методика расчета распространения аварийных выбросов, основанная на модели рассеивания тяжелого газа [Текст] / А. А. Шаталов, М. В. Лисанов // Безопасность труда в промышленности. - 2004. - № 9. -С. 46-52.

8. Hanna, S. Air Quality Modeling Over Short Distances [Text] / S. Hanna // College on Atmospheric

Boundary Layer and Air Pollution Modeling: 16 May - 3 June 1994. № SMR/760-2. - P. 712-743.

9. Murakami, S. Comparison of "k-e" model, ASM and LES with wind tunnel test for flow field around cubic model [Text] / S. Murakami, A. Mo-chida, H. Yoshihiko // 8th Int'l. Conf. on Wind Engineering (Western Ontario, July 8-11, 1991). -№ 12 - 3.

10. Belyaev, M. Air Pollution Modelling of Tech-nogenic Catastrophes Caused by Terrorist Attacks on Rail Transport of Chemical Substances [Text] / M. Belyaev // 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.

11. Biliaiev, M. M. Numerical simulation of indoor air pollution and atmosphere pollution for regions having complex topography [Text] / M. M. Biliaiev, M. M. Kharytonov // Conf. Abstracts of 31st NATO / SPS Int'l Technical Meeting on Air Pollution Modelling and Its Application (27 September - 01 October, Torino, Italy, 2010). - P. 1.7.

12. Tedeschi, G. Study of vertical transport of marine aerosol using an unsteady 2D model [Text] / G. Tedeschi // Conf. Abstracts of 31st NATO / SPS Int'l Technical Meeting on Air Pollution Modelling and Its Application (27 September - 01 October, Torino, Italy, 2010). - P. 4.9.

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

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