Моделирование распространения примеси в атмосфере с учетом рельефа местности Текст научной статьи по специальности «Математика»

УДК 519.6

П. Б. МАШИХИНА (ДИИТ)

МОДЕЛИРОВАНИЕ РАСПРОСТРАНЕНИЯ ПРИМЕСИ В АТМОСФЕРЕ С УЧЕТОМ РЕЛЬЕФА МЕСТНОСТИ

На 6a3i моделi потенщально! Te4ii' та моделi переносу домiшки запропоновано 2D чисельну модель для прогнозування забруднення атмосфери з урахуванням рельефу. Наведеш результати розрахунк1в на 683i розроблено! моделi.

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

The 2D numerical model to simulate the pollutant dispersion over complex terrain was proposed. The model is based on the equation of potential flow and the equation of admixture transfer. Results of the numerical experiment are presented.

Введение

Применяемые на практике методики прогноза загрязнения атмосферы при авариях на химически опасных объектах, транспорте, такие как нормативная методика [3, 7], методика «ТОКСИ» [8], модель Гаусса [1, 9], методика «РЬа8^а1еу» [10] не учитывают влияние рельефа местности на процесс рассеивания загрязнителя. В этой связи важной задачей является разработка эффективных методов прогноза рассеивания токсичных веществ в атмосфере с учетом рельефа. Необходимо отметить, что особенностью задач данного класса является требование минимальных затрат компьютерного времени на получение прогнозных результатов, поэтому разрабатываемые методы прогноза должны обеспечивать получение результата за реальное время (желательно, в течении нескольких минут) и в то же время учитывать основные физические процессы, влияющие на перенос загрязнителя. Целью данной работы является создание эффективной численной модели прогноза рассеивания загрязнителей в атмосфере, учитывающей наиболее существенные физические факторы, влияющие на этот процесс и позволяющей учесть рельеф местности при проведении вычислительного эксперимента.

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

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

пользовать модель течения идеальной несжимаемой жидкости [2]. В этом случае для расчета поля скорости воздушного потока необходимо найти решение уравнения для потенциала скорости [1]

д2 Р д2Р

дх ду2

= 0.

(1)

Для уравнения (1) ставятся следующие граничные условия [2]:

дР 0

• на твердых стенках — = 0 , где п - еди-

дп

ничный вектор внешней нормали;

• на входной границе (границы втекания

дР

воздушного потока) — = Vn, где

дп

Vn = f (у) - известное значение скорости ветра;

• на выходной границе

Р = Р * (х = const, у) + const (условие Дирихле).

Верхняя граница - твердая «стенка».

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

u =-

дР

дх '

v = -

дР

ду

На втором этапе моделирования осуществляется расчет распространения загрязнителя в атмосфере на базе уравнения переноса приме-

си, осредненного по ширине расчетной области

[1, 4]:

дС диС ЗуС

сС = (ц gradC) -

дt дх ду

N

+1 й ()5(х - х,. )(у - у),

,=1

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

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

Формирование вида расчетной области

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

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

£.113Е*В1 о в.1№С*В1 г л.9е*с*ее л е.9еес->ее 1 л.вздЕ'ве п е.?бЭЕ*ве л л.ьзеЕ'ве 1 е.ыек'ве с

-В.^ЭВЕ'ВЧ у в.Э27Е'89

а. 1-. Е ■■' • 4».1в9Е*в9

а.ЭЬЭС-А£

ш

н: : : г I г Н

Рис. 1. Форма двух облаков, образовавшихся на месте аварии

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

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

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

дР_ дп

д2 Р д2 Р

дх2 ду2

(2)

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

При п решение уравнения (2) будет стремиться к «установлению», т.е. к решению уравнения (1).

Для численного интегрирования уравнения (2) используется попеременно-треугольный метод А. А. Самарского [5]. Разностные соотношения в этом случае имеют вид

Р

и+1/2

— Р"

рп _ рп

Г,+1,1 Г1,1

- Р

п+1/2

лп+1/2

-1,1

0,5 Ап

Ах

рп _ рп

г1,1+1 г1,1

Ау2

_ Р"+12

Ах2

Р"+12 Г1,1 -1

Ау2

pn+1 pn+1/2 p.i - p.i

0,5 An

pn+1 pn+1 p+1,i - p, i

Ax2

- pn+1/2 + pn+lj 2

', i

pn+1 pn+1 pn+1/2

p,i+1 - p,i . -p,i

Ax2

pn+12 "p, i-1

Ay2

Ay2

Aj

температура, равная Т0 = 100 °С. Начальная температура в области - 0 °С. Исследуется процесс изменения температуры в точке с координатами х = 1 м, у = 1 м с течением времени £ Расчетная сетка 23 х 20 ячеек. Аналитическое решение задачи имеет вид [5]

На первом шаге расщепления находится «промежуточное» значение потенциала Р"+12 на

временном слое ««+1/2», а на втором шаге -определяется «окончательное» значение потенциала Р/'+1 на временном слое ««+1». Неизвестное значение Рна каждом шаге осуществляется по методу бегущего счета [6, 4].

Компоненты вектора скорости рассчитываются по соотношениям:

p - p P - P

1i, i 1i-1,i 1i, i 1i, i-1

и =—--—; v =-

y Ax ' y

T = To

1 - erf

• erf

2*Jat) v 2*Jat

У

В табл. 1 представлены данные расчета температуры в указанной точке по разработанной модели и на основе приведенного аналитического решения.

Таблица 1

Сравнение численных и аналитических расчетных данных

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

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

Тестовые расчеты

Для тестирования разработанной численной модели рассматривалась задача, имеющая аналитическое решение. Как известно, сопоставление численных расчетов с аналитическим решением является традиционным подходом при тестировании численных моделей. Такой подход применяется для тестирования как конечно-разностных моделей, так и моделей, основанных на методе конечных элементов [11]. Как известно, уравнение вида (2) называют уравнением теплопроводности (с коэффициентом температуропроводности равным единице а = 1). Рассматривалась следующая задача, имеющая аналитическое решение [5]: имеется расчетная область в виде прямоугольника. Длина области - 9,2 м, ширина - 7,6 м. На границах области задано условие первого рода -

Время t, Расчет по Расчет по разра-

с аналитической ботанной чис-

зависимости ленной модели

9,76 96,79 97,79

12,76 97,56 98,99

17,30 98,20 99,68

Как видно из табл. 1, имеет место удовлетворительное согласование результатов, полученных на основе аналитического решения и численного.

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

На основе построенной численной модели и алгоритма расчета разработан код на алгоритмическом языке FORTRAN. Разработанный код был применен для решения следующей задачи. Произошел аварийный выброс в атмосферу токсичного газа (цианистый водород). В результате выброса на месте аварии образовалось два облака сложной (грибовидной) формы (рис. 1). На пути миграции облаков располагается гряда холмов (рис. 1, 2). Высота холмов порядка 14,2 м. Концентрация примеси в облаке при t = 0 принята равной 1 (в безразмерном виде). Размеры расчетной области: длина -350 м, высота - 60 м. Коэффициент атмосферной диффузии по направлению x (продольное направление) равен 3 м2/с, а по направлению у - 2.8 м2/ с. Высота первого облака - 31,42 м, ширина (в шляпке) - 42 м. Высота второго облака - 19,95 м, ширина (в шляпке) - 42 м. За облаками находится насыпь, высотой порядка 11,4 м. Коэффициент с равен нулю. В модели учитывается неравномерный профиль ветра на входе в расчетную область [1]:

и = и

(* >

V У

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

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

где и1 - скорость ветра (4,44 м/с) на высоте (10 м); п = 0.15.

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

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

Рис. 5. Зона загрязнения атмосферы, t =37 с

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

В табл. 2 представлена динамика изменения максимальной концентрации токсичного газа в облаке.

Таблица 2

Динамика изменения максимальной концентрации токсичного газа в облаке

Время, с Максимальная концентрация (в безразмерном виде)

16 0,324

32 0,209

48 0,136

60 0,104

65 0,089

Как видно из табл. 2, максимальная концентрация в облаке после аварии уменьшается достаточно быстро.

В заключение отметим, что на решение задачи потребовалось 5 с компьютерного времени.

Выводы

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

дить в направлении ее совершенствования для

расчета рассеивания токсичных газов над водной поверхностью.

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

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

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

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

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

5. Полянин, А. Д. Справочник по точным решениям уравнений тепло- и массопереноса [Текст] / А. Д. Полянин. - М.: Факториал, 1998. - 368 с.

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

7. Хмшь, Г. А. Концептуально-методичний апарат анал1зу й оцшки техногенних та природних ри-

зишв [Текст] / Г. А. Хмшь // Еколопя довшлля та безпека жптед1яльност1. - 2007. - № 5. -C. 47-55.

8. Шаталов, А. А. Методика расчета распространения аварийных выбросов, основанная на модели рассеивания тяжелого газа [Текст] /

A. А. Шаталов, М. В. Лисанов // Безопасность труда в промышленности. - 2004. - № 9. -С. 46-52.

9. Hanna, S. Air Quality Modeling Over Short Distances [Текст] / S. Hanna // College on Atmospheric Boundary Layer and Air Pollution Modeling: 16 May - 3 June 1994. - № SMR/760-2 -P. 712-743.

10. Model Evaluation Report on UDM Version 6.0 [Текст] // Cambridge Research Consultants Ltd., Cambridge CB2 1SJ UK, Ref. No. SMEGIS/00/9/E, 2002. - 51 р.

11. Thomas, B. C. Comparison of Numerical Modelling Techniques for Complex, Two-Dimensional Transient Heat-Conduction Problems [Текст] /

B. C. Thomas, I. V. Samarasekkera, J. K. Bri-macombe // Metallurgical Transactions. B. - June 1984. - Vol. 15B, No. 2. Process. Metallurgic. -12 р.

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