CC BY
94
© А.А. Белослудцев, Д.В. Гусаров, М.А. Еремин, Н.М. Кузьмин, А.В. Хоперсков, С.С. Храпов, 2009
КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ
УДК 532.5; 519
ИНФОРМАЦИОННО-КОМПЬЮТЕРНЫЙ КОМПЛЕКС ДЛЯ МОДЕЛИРОВАНИЯ ДИНАМИКИ ПРИМЕСЕЙ ОТ ПРЕДПРИЯТИЙ ХИМИЧЕСКОЙ ПРОМЫШЛЕННОСТИ *
А.А. Белослудцев, Д.В. Гусаров, М.А. Еремин,
Н.М. Кузьмин, А.В. Хоперсков, С.С. Храпов
Работа посвящена вопросам определения характера загрязнения приземного слоя атмосферы от промышленных источников. В работе построена математическая модель для описания нестационарной трехмерной динамики загрязнений, в том числе от нестационарных источников для конкретного задаваемого физического состояния атмосферы. В предлагаемом подходе используется прямое численное интегрирование точных уравнений переноса примеси в атмосфере с учетом основных физических факторов, что приближает данный метод к проведению вычислительного эксперимента. На основе разработанной математической модели создана информационная система для компьютерного моделирования процесса распространения примесей от промышленных источников, расположенных на территории предприятия.
Введение
Расчеты динамики переноса примесей в атмосфере позволяют ответить на многие вопросы экологии промышленных зон. Получение карт загрязнений посредством расчета на ЭВМ особенно актуально применительно к веществам, для которых лабораторные измерения концентрации сопряжены с большими затратами, а также для примесей, не выбрасываемых непосредственно в атмосферу, но возникающих в ней в результате химических превращений обычных загрязнителей. С помощью компьютерных программ, основанных на прямом решении на ЭВМ математических уравнений, описывающих исследуемый процесс, строятся информационно-математические модели промышленных районов. Такие модели позволяют осуществлять экологический мониторинг за состоянием воздушно-водного бассейна, включая последствия различного рода аварий, газовок и т. п., примеры которых многочисленны.
Отличительной чертой математического моделирования, основанного на численном интегрировании нестационарных уравнений переноса, в отличие от официальных методик типа ОНД-86 (ОНД), является принципиальная возможность учета любых влияющих внешних факторов (метеоусловия, ветер, рельеф местности, солнечная радиация, химические и фотохимические превращения и т. п.) на «математическом уровне», то есть путем введения
* Работа выполнена при поддержке гранта РФФИ N° 10-07-97017-р_поволжье_а, ГК N° 02.740.11.5198.
требуемых факторов в систему математических уравнений. Заметим, что для построения карт загрязнения в рамках рассматриваемой модели не требуются многолетние наблюдения параметров состояния атмосферы и распределения массы загрязняющих веществ по источникам выбросов. Достаточно данных на текущий момент. Это позволяет учитывать последствия нестационарных процессов - аварии, газовки, изменение метеоусловий со временем и т. п. Проблема оценки уровня загрязнения нижних слоев атмосферы и подстилающей поверхности является одной из центральных в экологии промышленных районов. Теоретические исследования, основанные, как правило, на численном решении модельных уравнений и статистическом анализе, позволяют во многих случаях дать приемлемые прогнозы распределения концентраций примесей (см.: [3; 8; 17; 19].
Компьютерная модель предназначена для получения оперативной информации, необходимой для улучшения экологической обстановки в регионе. Ее можно использовать в качестве одного из инструментов, позволяющих проводить контроль предприятию за загрязнением окружающей среды. Система может служить основой при проведении природопользовательских и экологических экспертиз при размещении различных производств. Возможно определение метеорологических условий, при которых загрязнение приземного слоя атмосферы максимально и минимально.
Математическая модель динамики примесей
Математическая модель основана на численном интегрировании системы нестационарных уравнений в частных производных с источниками, которая описывает процесс переноса загрязняющих примесей с учетом основных физических факторов. Основными компонентами компьютерной системы являются:
1) Модель приземного слоя атмосферы.
Модель приземного слоя атмосферы учитывает метеорологические и климатические условия (направление и силу ветра, температуру воздуха, конвективное состояние атмосферы). Скорость ветра и амплитуда турбулентных пульсаций задается в м/с. Направление ветра задается в градусах от западного направления, либо выбирая из 16 определенных направлений. Конвективное состояние приземного слоя определяется в первую очередь вертикальным градиентом температуры.
2) Модель динамики загрязняющих примесей.
Расчет динамики распространения примесей учитывает основные физические факторы, определяющие этот процесс: ветровой снос, диффузионный перенос, осаждение на поверхность земли и повторное загрязнение. Значения коэффициентов турбулентной диффузии задаются на основе стандартной модели ОНД-86 (ОНД).
На процесс распространения примесей в атмосфере оказывают влияние множество факторов, но основными являются два: ветровой снос и диффузионное расплывание. Эти два физических эффекта можно учесть в рамках математической модели, в основе которой лежит дифференциальное уравнение в частных производных:
= -V[C(r,t) рк]+ V(D(r, Г) Vpk)+Q(r, t) -a(r, t)pk, (1)
ot
где t - время;
r - радиус-вектор;
pk - плотность k-й примеси;
с - вектор скорости ветра;
D(r, t) - тензор турбулентной диффузии;
Q(r, t) - совокупность источников и стоков рассматриваемой примеси;
а - параметр, описывающий «вымывание» примеси из атмосферы за счет различных факторов и вторичное загрязнение от подстилающей поверхности;
V = (5 / Ох, 5 / Оу, 0 / Oz) - дифференциальный векторный оператор набла.
Источники загрязнения устанавливаются функцией Q, которая определяет временную динамику выбросов, их параметры (координаты и мощность выбросов, высоту труб, температуру; для плоскостных источников задается пространственное распределение загрязнителей). Уравнение (1) записано в векторной форме. Для решения задачи распространения примеси в приземном слое атмосферы удобно работать в декартовой системе координат: х и у - координаты в плоскости поверхности земли, а г - вертикальная координата. Слагаемое др/дt описывает изменение концентрации примесей со временем в точке с координатами (х, у, г).
Рис. 1. Схема, поясняющая принцип и структуру работы созданного программного комплекса. Стрелками указаны направления потока данных
Запишем уравнение (1) в декартовой системе координат:
др
д
Я, ЛСх (Г-Ґ) Р]-^[СУ (Г -Ґ) р]
дґ дх ду
+
д
+ — ІРх(г - 0'^р)+ ду (ру -^)^р)+ дТ" Р (г-ґ )^р)+ 0$, ґ)- а(г,ґ )р.
(1а)
Величины В , В , В , С , С являются эмпирическими параметрами, учитывающими как со-
X’ У 2’ X’ у г ^ ^ ^
стояние атмосферы, так и рельеф местности. Характер функциональной зависимости коэффициента диффузии и скорости ветра в работах различных авторов принимается приблизительно одинаковым. Можно считать, что коэффициент диффузии В на небольших высотах линейным образом зависит от величины скорости ветра, не обращаясь, однако, в нуль при | С |= 0. Поэтому в нашей модели принято
Р = ■
]ртт с / с
1 Р™
при | С > Ст при | С |< Ст
(2)
В зоне перемешивания для коэффициента диффузии имеем:
1 + Р0 • г • ехр< - 4,0
\С(г) +1)
(3)
г
г
И
где значения параметров В0, zp сказываются только на высотах г < (2 - 3) zp (Zp - высота пограничного слоя). В зависимости от стратификации атмосферы высота эффективного перемешивания в вертикальном направлении может сильно меняться в пределах от 50 до 1 000 м.
Первое слагаемое справа в формуле (1) обусловлено ветром в плоскости земли С = (Сх,С ):
С = С 0 + С1, (4)
где необходимо учитывать наряду с постоянной составляющей С0 и флуктуирующую часть С1. Среднее направление ветра в каждой точке плоскости (х, у) характеризуется углом ф0 между осью х и направлением ветра, а дисперсия направлений задается параметром сф. Тогда при нормальном распределении направлений ветра
Р(ф) = ~г== ехр(-(ф- Ф0)2 /2^ф)
“~ф> (5)
Зависимость скорости ветра от вертикальной координаты для базовой модели выбираем в виде:
ад = с .(1+0,434• Ц1+2/20}), (6)
где величина С - характеризует скорость ветра у поверхности Земли; параметр 20 - определяет шероховатость поверхности. В общем случае можно задавать произвольное распределение С(2).
Модель источников и стоков
К исходным данным для проведения расчетов по этой методике необходимо знание источников загрязнения - координаты, высота, химический состав и объем выбросов, включая пространственно распределенные источники - пруды-накопители, дороги, карьеры. Информация по источникам соответствует принятой в общероссийской методике ОНД-86.
Характер источников может быть произвольный, определяется целями расчета, в частности, различают стационарные и нестационарные, точечные, плоскостные и объемные источники. Основными источниками для модели горения прудов-накопителей являются плоскостные источники, а при моделировании динамики выбросов от промышленных предприятий - точечные источники.
Высота выброса примеси в атмосферу, дальнейшее ее рассеяние и осаждение зависят не только от высоты источника (трубы), но и от таких параметров примеси, как температура и скорость ее вытекания из трубы. Пусть труба имеет высоту Н0, ее радиус на выходе примеси -R0, скорость выноса из жерла трубы - ю0, температура атмосферы - Та и температура примеси -Та + АТ. Тогда эффективная добавка АН к высоте трубы Н0 (Н = Н0 + АН) удовлетворительно аппроксимируется соотношением:
АН = 3,75 +1,6 МАГ *о. (7)
Здесь первый член описывает вертикальную инерцию примеси, а второй - ее дополнительную «архимедову» плавучесть.
В зависимости от величины и знака вертикального градиента температуры, исходящая из локализованных и распределенных источников примесь может либо «всплыть» в верхние слои атмосферы (турбулентная конвекция), либо прижиматься к земле и накапливаться в приземном слое толщиной порядка высоты деревьев и зданий. И в первом, и во втором случае соответствующие эффекты можно описать с помощью модельного параметра а. Однако второй из указанных факторов, как наиболее опасный в экологическом отношении, необходимо учесть с помощью локального уравнения «накопления»:
д# „
-Г7 = а зР-р#, (8) дt
где # - локальная плотность накапливающейся в приземном слое атмосферы примеси;
Р - параметр, описывающий полное «вымывание» накопленной примеси, а3 > 0.
Как известно, абсолютное значение вертикального градиента температуры и его знак характеризует температурную стратификацию, которая является важным фактором, способствующим увеличению или уменьшению концентрации вредных примесей. Влияние температурной устойчивости атмосферы на турбулентный обмен оценивается по безразмерной величине Р, аналогичной по физическому смыслу числу Ричардсона К1. Последняя величина характеризует отно-
__ п Т (г = 0,2 м) - Т (г = 1,5 м)
шение разности температур на выбранных уровнях к скорости ветра Р = т ■ .
С (г = 2м)
Можно говорить, что атмосфера стратифицированно неустойчива при р > 0,06, устойчива при Р < 0,05 и нейтральна при р « 0.
Общая структура информационной системы показана на рисунке 1.
Важным компонентом модели является электронная карта. Вся моделируемая область покрывается сеткой размером N х М х К. В вертикальном г-направлении в простейшем случае достаточно ввести несколько слоев (К > 5), которые учитывают перенос примесей в вертикальном направлении и взаимообмен между приземным слоем атмосферы, верхними слоями и подстилающей поверхностью в соответствии с вышеизложенным. С построенной моделью проводились несколько серий тестирований. Во-первых, исследовалась роль внутренних параметров модели (выбор сетки, временной шаг интегрирования, граничные условия). Во-вторых, варьировались физические условия моделирования (ветровой и диффузионный режимы, стратификация атмосферы, режим вымывания примесей из атмосферы). Наконец, рассмотрены различные режимы работы загрязняющих источников (см. рис. 2 и 3).
I “О ЭноКаустин 1.2.8.28 - СЬІог ЕИ®|
Файл Карта Проект Модель Расчет Просмотр результатов Источники Справка
Д [ш] Щ \20%
л Каргга Объекты Проект Источники проекта
Источники проекта (2 шт)
Цех№30: Дефлектор.Столыхраненияпроб
Свойство Описание Значен^^ІЛ Колонна ней-
Количество СП источников выделения загрязнаяющих веществ (шт) 1
У'&У^/-'/'к Количество П источников выделения загрязнаяющих веществ (шт) 1
Высота СП источника выброса, м 15
Высота П источника выброса, м 15
Температура СП газовоздушной смеси на выходе из источника выброса. гр.С 33
Температура П газовоздушной смеси на выходе из источника выброса. гр.С 33
Скорость СП газовоздушной смеси на выходе из источника выброса. м/с Є.9315
Уф /у)К *У у Скорость П газовоздушной смеси на выходе из источника выброса, м/с Выбросы загрязняющих веществ СП, г/с 6,9315 0,0014
Выбросы загрязняющих веществ П, г/с 0,0014
Чж# УУ+У Выбросы загрязняющих веществ СП. мг/мЗ 1,0213
Выбросы загрязняющих веществ П, мг/мЗ 1.0213
Выбросы загрязняющих веществ СП. т/год 0.0402
і X % І Выбросы загрязняющих веществ П. т/год 0.0402
Диаметр СП устья трубы, м 0.5
Диаметр П устья трубы, м 0.5
Координата X 4101
Координата V 2327
Время работы в фоновом режиме •1
Время работы в режиме выброса -1
Режим работы источника - начало исследования (время ТО) і 1
Режим работы источника • начало исследования (мощность МО) 0,0014
Режим работы источника • начало выброса (время Т1) 2
Режим работы источника - начало выброса (мощность N1) 3.3347
Подключено X: 4*102 ¥: 1714,6
Рис. 2. Вид главного окна программы «ЭкоПром» с открытым окном «Источники проекта»
Рис. 3. Вид главного окна программы «ЭкоПром» во время проведения расчетов
Численный алгоритм решения
Опишем алгоритм построения численного решения уравнения адвекции - диффузии на примере одномерного уравнения вида:
д0 + п дР Л[ D(x) ^1 + Q(х, ґ) дt дх д{ дх ' '
(9)
Явная конечно-объемная схема для уравнения (1) имеет вид:
»+1 _ о гп _ гп
Ні Ні + ] і +1/2 Ji _1/
(
Дґ
Дх
Дх
_п _п _п _п Л
О 0+1 _ 0І _ О 0І _ 0_1
иі+\и А і_1/ 2 .
ч Дх Дх у
+ап
(10)
где индекс і обозначает номер ячейки; п - номер временного слоя; индексами і ± 1/2 обозначены величины, относящиеся к границам ячеек;/= п х р - поток.
Схема (10) устойчива при ограничении на шаг интегрирования согласно условию:
(
dt < тіп
Дх Дх
Для подавления осцилляций, возникающих при численной аппроксимации адвективных слагаемых (второе слагаемое в 1 или 2), потоки через границы ячеек аппроксимируются TVD-схемой
Ки 2 = 0,5а + /м) - 0,5 \им, 2|[1 -ф(г,+17 2)](Р^1 -Р),
где ф(г) - функция-ограничитель;
~ п _^п
- _Рі_Р,_1
Рі +1 _Р,
отношение прилежащих градиентов.
1
В тестовых расчетах были использованы различные варианты ограничителей, мы остано-
2| Г
вили свой выбор на ограничителе ван Лира ф(г) =----, представляющем хороший компромисс
г + 1
между численной диффузией решения и затратами на его вычисление. Расчеты, проведенные по схеме (10), показали полное отсутствие осцилляций численного решения даже при больших числах Пекле, когда адвекция доминирует над диффузией.
Основные компоненты информационно-компьютерного комплекса
Основными компонентами данной информационной системы являются:
1. Основная графическая оболочка управления информационной системы.
2. Расчетный модуль для численного моделирования нестационарной трехмерной динамики примеси.
3. Картографический блок территории предприятия, включая санитарно-защитную зону предприятия.
4. Электронную базу данных источников промышленных выбросов.
5. Электронную базу данных входных параметров для расчетного модуля, который включает подсистему метеорологических и физических параметров.
6. Модуль визуализации динамики примеси непосредственно во время проведения моделирования.
7. Графический блок для визуализации проведенных расчетов.
8. Блок для генерации отчетов.
СПИСОК ЛИТЕРАТУРЫ
1. Белов, И. В. Транспортная модель распространения газообразных примесей в атмосфере города / И. В. Белов, М. С. Беспалов, Л. В. Клочкова [и др.] // Математическое моделирование. - 2000. - Т. 12, № 11. - С. 38-46.
2. Белослудцев, А. А. Информационная система для создания электронных карт с рельефом ^еоМарВшШег) / А. А. Белослудцев, Н. В. Наводченко, А. В. Хоперсков // Свидетельство об официальной регистрации программ для ЭВМ № 2007612682 от 21.06 2007 г.
3. Берлянд, М. Е. Прогноз и регулирование загрязнения атмосферы / М. Е. Берлянд. - Л. : Гидрометео-издат, 1985. - 272 с.
4. Берлянд, М. Е. Современные проблемы атмосферной диффузии и загрязнения атмосферы / М. Е. Берлянд. - Л. : Гидрометеоиздат, 1975. - 418 с.
5. Гуди, Р. М. Атмосферы : пер. с англ. / Р. М. Гуди, Дж. Уолкер. - М., 1975.
6. Замай, С. С. Модели оценки и прогноза загрязнения атмосферы промышленными выбросами в информационно-аналитической системе природоохранных служб крупного города / С. С. Замай, О. Э. Яку-байлик. - Красноярск : Изд-во Краснояр. гос. ун-та, 1998. - 109 с.
7. Крукиер, Л. А. Математическое моделирование процессов переноса в несжимаемых средах с преобладающей конвекцией / Л. А. Крукиер // Математическое моделирование. - 1997. - Т. 9, № 2. - С. 4-12.
8. Марчук, Г. И. Математическое моделирование в проблеме окружающей среды / Г. И. Марчук. - М. : Наука, 1982. - 320 с.
9. Матвеев, Л. Т. Основы общей метеорологии. Физика атмосферы / Л. Т. Матвеев. - СПб. : Гидрометеоиздат, 2000. - 751 с.
10. Методика оценки последствий аварийных выбросов опасных веществ (Методика «Токси». Редакция 3.1). - М., 2005.
11. Методика расчета концентраций в атмосферном воздухе вредных веществ, содержащихся в выбросах предприятий. Общесоюзный нормативный документ (ОНД-86) / Госкомгидромет. - Л. : Гидрометеоиздат, 1987. - 68 с.
12. Методическое пособие по расчету, нормированию и контролю выбросов загрязняющих веществ в атмосферный воздух. - СПб. : НИИ Атмосфера, 2005. - 215 с.
13. Митько, Л. А. Компьютерное моделирование физико-химических процессов в приземном слое атмосферы - необходимый компонент информационного обеспечения в системе эколого-гигиенической безопасности / Л. А. Митько. - Киев : Ин-т проблем моделирования в энергетике НАН Украины, 2005. - 221 с.
14. Никифоров, А. Н. Моделирование полей загрязненности атмосферы в мезометеорологическом пограничном слое / А. Н. Никифоров, Н. С. Бузало // Изв. вузов. Северо-Кавказский регион. Естественные науки. - 2001 (спецвыпуск). - С. 126-128.
15. Предельно допустимые концентрации вредных веществ в воздухе и воде. - Л. : Химия, 1975.
16. Предельно допустимые концентрации (ПДК) загрязняющих веществ в атмосферном воздухе населенных мест. ГН 2.1.6.1338-03.
17. Хоперсков, А. В. Математическое модел^ование экосистем / А. В. Хоперсков // Социальная экология : сборник. - Волто^ал : Изд-во ВолГУ, 1993. - С. 88-91.
18. Хргиан, А. Х. Физика атмосферы / А. Х. Хргиан. - М. : Изд-во МГУ, 1986. - 328 с.
19. Экологический программный комплекс для персональных ЭВМ / под ред. А. С. Гаврилова. - СПб. : Гидрометеоиздат, 1992. - 165 с.
20. Kamenetsky, E. Model of the flow and air pollution concentration in urban canyons / E. Kamenetsky, N. Viern // Boundary Layer Meteorol. - 1995. - V. 73, N° 1-2. - P. 203.
Summary
THE PROGRAM COMPLEX FOR MODELING AIR POLLUTANT DYNAMICS
EJECTED BY CHEMICAL PLANTS
A.A. Belosludtsev, D. V. Gusarov, M.A. Eremin,
N.M. Kuzmin, A.V. Khoperskov, S.S. Khrapov
This work is devoted to study ground layer pollution dynamics from industrial sources. The mathematical model is constructed for description the nonsteady three dimensional pollutant dynamics at different regimes of nonsteady ejections for a given meteorological atmosphere state. The direct integration of exact advection equations for pollutants is used with taking into account main physical factors. On the basis of developed mathematical model the program complex is created for computational modeling of pollutants dynamics from industrial sources situated at the plants territories.
