К расчету химического загрязнения атмосферы с учетом сложного рельефа местности Текст научной статьи по специальности «Математика»

УДК 519.6 : 504.3.054

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

К РАСЧЕТУ ХИМИЧЕСКОГО ЗАГРЯЗНЕНИЯ АТМОСФЕРЫ С УЧЕТОМ СЛОЖНОГО РЕЛЬЕФА МЕСТНОСТИ

На базi тривимiрноl чисельно! моделi виконано розрахунок процесу поширення токсично! речовини в атмосферi з урахуванням рельефу при аварiйному розливi. Модель базусться на чисельному iнтегруваннi рiвняння конвективно-дифузiйного переносу домiшки та моделi потенцiйного руху. Наводяться результати обчислювального експерименту.

На базе трехмерной численной модели выполнен расчет распространения токсичного газа в атмосфере с учетом рельефа после аварийного разлива. Модель основывается на численном интегрировании уравнения конвективно-диффузионного переноса примеси и модели потенциального течения. Приводятся результаты вычислительного эксперимента.

The 3D numerical model was used to simulate the toxic gas dispersion over a complex terrain after an accident spillage. The model is based on the K-gradient transport model and the model of potential flow. The results of numerical experiment are presented.

Введение

Рост требований к качеству прогнозной информации по оценке последствий аварийных ситуаций, сопровождающихся выбросом (разливом) химически опасных веществ ставит задачу создания моделей, учитывающих при расчете такой важный фактор как рельеф местности [1]. Нормативная методика [5], используемая для прогноза последствий при авариях с опасными веществами, не учитывает влияние рельефа местности на процесс рассеивания примеси в атмосфере и поэтому не может быть применена для прогноза масштаба загрязнения для условий сложной орографии. Применение аналитических моделей рассеивания примеси, например, модели Гаусса (данная модель реализована во многих коммерческих программных продуктах, в частности таких как «Аммиак», ISC3, SCIPUFF, SDM, PUFF-PLUME, HYROAD и т.д.) - также исключено, т.к. аналитические модели не учитывают деформацию поля скорости ветрового потока (первоочередное влияние рельефа) при обтекании сложного рельефа, а значит применение их - принципиально невозможно в условиях сложной орографии. Математическое моделирование рассеивания загрязнителя в атмосфере с учетом рельефа местности может быть выполнено только с помощью численных моделей (CFD модели) [8 - 10].

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

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

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

Моделирование процесса переноса токсичного газа в атмосфере проводится на основе уравнения переноса примеси [2, 3]:

дС диС дуС д(w - wS)С _

dt дх

ду

dz

д . дС, д . дС, д . дС.

=—(цх —)+—(к —)+—(м> —)■

vr х ^ ' vr y л / \Г z s

дх дх ду ду дz дz

+Z Q (t Жr - r),

(1)

где С - концентрация примеси, попадающей в атмосферу при аварии; и, V, w - компоненты вектора скорости воздушной среды; ws - скорость оседания примеси (в данной работе принята равной нулю); ц = (ц , цу , - коэффициент турбулентной диффузии; Q - интенсивность выброса токсичного вещества от зоны

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

разлива; 5( г - г) - дельта-функция Дирака;

г^ = (х , у^ , zi) - координаты источника выброса.

Постановка краевых условий для данного уравнения рассмотрена в работах [2, 4]. Если рассматривается задача о миграции облака, образовавшегося на месте аварии, то в начальный момент времени задается форма облака и концентрация загрязнителя в нем.

В построенной численной модели используется неравномерный профиль скорости ветра и вертикального коэффициента диффузии:

u = u,

( г >

V г1 J

ц г = 0,11г,

dt

Ai

ôx =--1 ÔX h-• ôx

ôvC ÔV+C dv~C

0У ày dy

ôwC dw+С dw~ С t--

дг

u + u

дг

u — u

дг

v+v

где u =-— ; u =-L

2 2

v = -

v =

v —M + w + \w\ — w — \w\ w =-

■; w =■

2 2 2

Для аппроксимации конвективных производных используем выражения:

ди + С ~ 1 •к - иЬк 1 ,к = ¿+ сп+1 .

дх

Ax

дu С ^ ui+1, j,k Ci+1, j,k uijk Сi jk = L— С n+1 ;

д х

Ах

= L— С" ,

dv+С ^ vi, j+1,kCijk vijkCi, j—1,k = c"+1 • ду Ay y '

дv С ^ vi, j+1,kCi, j+1,k vijkCijk = j-c"+1 ; ду Ay y

dw+C wthk+1C,jk — w+C j.k—1

дг Аг

ôw~C ^ wi, j,k+1Ci, j,k+1 — wijkCi, j,k дг ~

= L+C"

Аг

= L—Cn+1.

где и1 - скорость ветра на высоте ^ = 10 м; п = 0,15.

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

Произведем следующую аппроксимацию производных, входящих в уравнение переноса примеси. Заменим производную по времени разделённой разностью «назад»:

дС СП+ - сппк

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

_ / _ ~ \ ^п+1 ^п+1 ^п+1 ^п+1

д ( дС \ _ Ч+1,1,к - СЦк _ Ч, 1,к - Сг-\,з,к

—I - РМ* х-—2—--1-1 х—-2— =

дх I дх I А х А х

=M—xCn+1 + M+C

n+1 .

дС_

ду V у ду

fi n+1 _г~т+1 г~т+1 _ s~m+1

- Ci, j+1,k~Cijk - Ci, j,k~Ci, j—1,k

y-r^--^ y"

A y2

= MyCn+1 + M+yCn+1;

A y

2

Конвективные производные представим в виде:

дuС du+С du—С

д ( дС \ _ Ci, jkk+1 — Cijk _ Ci, j,k — Cijkk-1

ъ V" г ^ r —ц г-^г^

n+1

= MZCn+l + M+C

гг гг

В данных выражениях LLx, LLX, LLy, LLy, L+, L—,

M+, MXcx и т.д. - обозначения разностных операторов. С учетом этих обозначений разностный аналог трехмерного уравнения переноса примеси будет иметь вид:

с n+1 с п

ijk — ijk + L+ с n+1 + LX с n+1 + JL с n+1 +

At x x y

+L— С n+1 + L+ С n+1 + L— С n+1 +a =

y г г ijk

= ( M+хс n+1 + M-хс n+1 + M+ycn+1 +

+M- сn+1 + M+cn+1 + M-сn+1 ) .

yy г zz J

Расщепим решение данного разностного уравнения при интегрировании на временном интервале dt так:

k 1

- на первом шаге k = — :

4

n+k n

С ij — С ij

A t 2

+1 ÍLt С k + L+, С k + L + С k ) +

л \ x У z '

2

4 1 4

4 (

+ - С * = -{ М+х Ск + М- Сп + М+У Ск +

уу

+Му Сп + М+2 Ск + М-22 Сп );

7 1 1

- на втором шаге к = п +—; с = п + —:

(3)

ленного интегрирования этого уравнения используется идея установления решения по времени, т.е. интегрируется уравнение эволюционного вида

дР д2 Р д2 Р д2 Р

дп дх2 ду2 д22

(6),

с к -с ']к 1

Д/

2 (

т V

+1 (( ск + гу ск + г ск) +

+ - Ск1. =1 (Мх Ск + М+х Сс + М~ Ск +

4 1 ] 4 \ ^ хх у

+ М+у Сс + М- Ск + М+ Сс);

уу 22 22 I '

(4)

1 3 1

- на третьем шаге к = п +—; с = п +— ис-

4 2

пользуется формула (4);

7 1 3

- на четвертом шаге к = п +1; с = п + — используется формула (3)

На пятом шаге расщепления учитывается влияние источников на изменение концентрации примеси, и расчетная зависимость на этом шаге имеет вид:

где п - фиктивное время.

При п решение уравнения (6) будет стремиться к «установлению», т.е. к решению уравнения (5) [7].

Численное интегрирование проводится с помощью схемы суммарной аппроксимации [7]. В этом случае разностные уравнения на каждом шаге расщепления имеют вид: - на первом шаге:

1

п +— Р 2 — Рп

1,1,к Ч ,],к Д/

1

п+—

— Р 2 1,1,к

+ Р 2 + 4—1,1,к

Дх2

1 1 п+— п+— -Р 2 + Р 2

4,1,^4,1—1,к

Д У

2

1

1

п+— п+--Р 2 + Р 2

Д г

2

5 п+1 5 п

ск—сиик_д, (/п+12) 5 1=1 Д х Д у д 2 1

д/

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

Для расчета скорости ветра при обтекании рельефа используется модель потенциального течения. Как известно [3], уравнение для определения потенциала скорости имеет вид:

д2 Р д2 Р д2 Р

д х2 д у2 д г2

= 0.

(5).

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

1 п+—

рп+1 _ р 2

1,1,к 1,1, к

Д/

рп+1 рп+1

Р+\,3 , к — '1,1,к

рп+1 _ Г>п+1

Р ,1+1, к — Р ,1,к

Д у2

Д х2

рп+1 _ рп+1

Р ,1,к+1 — Р ,1,к

Д 2 2

На каждом шаге расщепления неизвестная величина потенциала скорости определяется по схеме бегущего счета:

- на первом шаге расчетная зависимость имеет вид:

Компоненты вектора скорости определяются по зависимостям

оР =дР дх ду

Постановка граничных условий для данного уравнения рассмотрена в работе [3]. Для чис-

1

п+

= РЫ + Д/ •

1 1 п+— п+— -Р 2 + Р 2

1 1 п+— п+— — Р 2 + Р 2

4,1,к ^4—1,1,к

Дх

2

Д У2

1

1

п+— п+--Р 2 + Р 2

Д 2

2

- на втором шаге расчетная зависимость имеет вид:

1

n+— _ 2 .

i, j, k i, j,k

pn:i _ p

At •

P

n+1

- P

n+1

i+1, j, k i, j ,k

2

P

- P

n+1

n+1

i, j+1,k 1 i, j ,k

A У

2

A x

Т)П+1 т)П+1 Pijk+1 — Pi, j,k

A z2

В численной модели используется метод маркирования («porosity technique») для выделения разностных ячеек, которые относятся к области течения и к твердым границам (рельеф) [9]. Применение этого метода дает возможность быстро «построить» форму рельефа местности при вводе исходных данных для моделирования.

На основе рассмотренной численной модели создан код TERRAIN-3, реализованный на алгоритмическом языке FORTRAN. Разработанный код построен на модульном принципе, что дает возможность его применения для решения широкого круга задач, без внесения корректировки в структуру программ.

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

Рассматривается применение построенной трехмерной численной модели для прогноза загрязнения воздушной среды при гипотетическом аварийном разливе синильной кислоты на станции Днепропетровск-Южный (рис. 1, зона аварии условно показана «кружком», В - здание вокзала). Ветер направлен в сторону здания вокзала. Зона разлива находится на склоне насыпи (рис. 2) . Моделирование выполнено при следующих данных: размеры расчетной области 90 м х 90 м х 36 м; V10 = 3,8 м/с; цх= цу = = 1 м2/с; ц = 0,11-Z (Z - высота); площадь зоны разлива составляет порядка 112 м2. Интенсивность испарения от зоны разлива рассчитывалась в коде по зависимости [5]:

Q _ (5,83 + 4,1 • V) РП4М,

(7)

противоположной (северной) стороны. Расчет поля скорости и переноса примеси, как в атмосфере, так и внутри здания осуществляется сквозным счетом.

Рассмотрим результаты моделирования. На рис. 3 представлена зона загрязнения атмосферы возле вокзала для момента времени ^ = 80 с. Хорошо видно, как шлейф токсичного газа вытягивается вдоль склона насыпи в направлении вокзала и внутри здания вокзала формируется зона загрязнения. Также отчетливо видно, что над участком, где располагаются пути, образуется обширная область загрязнения. Наличие насыпи и здания приводит к деформации формирующейся зоны загрязнения.

где Q - интенсивность испарения,

V - скорость ветра;

Рн - давление насыщенных паров;

М - молекулярная масса вещества.

Данная величина составляла порядка 4,3 кг/с от зоны разлива.

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

Рис. 1. Место гипотетической аварии на ст. Днепропетровск-Южный (В - вокзал)

Рис. 2. Схема расчетной области: а - зона аварийного разлива на склоне насыпи; с - место расположения железнодорожных путей; в - здание вокзала

Рис. 3. Зона загрязнения атмосферы t = 80 с (сечение y =41,25 м)

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

Таблица 1 Изменение концентрации токсичного газа

в здании вокзала

t, с С, г/м

14,6 1,51

20 2,08

36 3,02

Отметим, что для решения данной задачи потребовалось 10 с.

Выводы

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

Дальнейшее совершенствование рассмотренной в работе модели необходимо проводить в направлении ее адаптации к моделированию процесса нейтрализации токсичного газа.

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

1. Аварии и катастрофы. Предупреждение и ликвидация последствий [Текст] : учеб. пособие в 5-ти кн. / под ред. В. А. Котляревского и

A. В. Забегаева. - М.: Изд-во АСВ, 2001. -200 с.

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

3. Лойцянский, Л. Г. Механика жидкости и газа [Текст] / Л. Г. Лойцянский. - М.: Наука, 1978. -735 с.

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

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

6. Мацак, В. Г. Гигиеническое значение скорости испарения и давления пара токсических веществ, применяемых в производстве [Текст] /

B. Г. Мацак, Л. К. Хоцянов. - М.: Медгиз, 1959. - 231 с.

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

8. Castino, F. Parameterization Of Convective And Stable Interval Boundary Layers Into Mass Consistent Models [Text] / F. Castino, M. Tombrou // 2 EACWE (Genova, Italy, 1997). - P. 309-314.

9. 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 University of Technology, Netherlands, July 2-6, 2001). - P. 355-358.

10. Maurizi, A. Numerical Simulation of the Atmospheric Flow in a Mountanious Region of North Portugal [Text] / A. Maurizi, J. M. L. M. Palma, F. A. Castro // Proc. of 2nd EACWE (Genova, Italy, 1997). - Vol. 1. - P. 293-300.

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

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

3