Введение Величина прямого ущерба от загрязнения атмосферы выбросами промышленных предприятий, тепловых электростанций, выхлопными газами автомобилей и других атмосферных загрязнителей, по оценкам ряда авторов, ежегодно составляет сотни миллиардов долларов. От последствий загрязнения атмосферы ежегодно умирает примерно 2,1 миллиона человек, еще около 470 тысяч погибает из-за разрушения озонового слоя. Один из перспективных путей частичного решения указанных проблем непосредственно связан с применением аппарата системного анализа, теории управления и современных средств обработки информации, используемых для мониторинга процессов формирования и распространения атмосферных поллютантов, прогнозирования степени загрязненности контролируемых объектов и территорий, а также для управления технологическим оборудованием промышленных предприятий с целью снижения ущерба, причиняемого их атмосферными выбросами. В результате практического применения результатов этих исследований созданы и хорошо зарекомендовали себя на практике многие современные математические, аппаратные и программные средства экологического мониторинга и управления, например, единая государственная система мониторинга окружающей среды Российской Федерации, системы экологического мониторинга США, Европейского Союза и пр. Прогноз динамики распространения химически опасного вещества, выброшенного в атмосферу в результате аварии, является важной задачей в области экологической безопасности. Целью такого прогноза является получение информации для оценки последствий аварийных ситуаций, а именно: определение размеров зон токсичного поражения, близости этих зон к жилой застройке, скорости распространения, а также возможности взрыва или возгорания. Данная информация необходима для оперативного устранения последствий аварий с химически опасными веществами, а именно: выбора способа их нейтрализации, маршрутов эвакуации, расчета критичных сроков эвакуации. Кроме того, возможность моделирования аварии при различных природных условиях служит основой для разработки адекватных мер защиты населения и окружающей среды [1-6]. Задачей нашего исследования является разработка математической модели для анализа динамики распространения химически опасного вещества, выброшенного в атмосферу в результате аварии, и прогнозирования его концентрации в заданных точках контролируемой территории. Математическая модель Для решения поставленной задачи рассмотрим прогностическое уравнение турбулентной диффузии, представленное в [2]: (1) где с - искомая концентрация токсичного вещества; x, y, z - координаты точки расчета по осям абсцисс, ординат и аппликат соответственно; t - время; u, w, v - проекции вектора средней скорости перемещения вещества на оси абсцисс, ординат и аппликат соответственно; α - коэффициент изменения концентрации вещества из-за химических превращений; kx, ky, kz - составляющие коэффициента обмена по осям x, y, z соответственно; Q - интенсивность выброса токсичного вещества; δ(х - xs), δ(y - ys), δ(z - zs) - дельта-функция Дирака; xs, ys, zs - координаты источника эмиссии токсичного вещества. Уравнение (1) должно быть дополнено начальным и тремя граничными условиями по каждой из пространственных координат: Данная математическая модель позволяет вычислить концентрацию химически опасного вещества в любой точке исследуемого пространства с течением времени, учитывая начальные данные, направление ветра и интенсивность выброса вещества в атмосферу. Построение разностной схемы Решение уравнения (1) предлагается искать численным методом. Построим разностную схему для данного уравнения. Пусть для независимых переменных заданы следующие интервалы их изменения: Получим разностную сетку, разбивая каждый из этих интервалов на некоторое количество равных частей. Введем следующие обозначения: - n, i, j, k - порядковые номера точек деления по t, x, y, z соответственно; - - интервалы между точками по t, x, y, z соответственно; - - значение функции c, соответствующее точкам tn, xi, yj, zk; - - значение функции c, соответствующее точкам tn, xi±1, yj, zk ; - - значение функции c, соответствующее точкам tn+1, xi, yj, zk. Используя введенные обозначения, а также аппроксимацию дифференциальных операторов [7], составляющих уравнение (1) в точке c(tn, xi, yj, zk), запишем явную разностную схему, аппроксимирующую уравнение (1): (2) Очевидно, что из полученного выражения можно явно выразить искомую величину концентрации в каждый следующий момент времени . Однако для возможности использования данной разностной схемы необходимо исследовать ее устойчивость. Исследование устойчивости разностной схемы Исследовать устойчивость явной разностной схемы (2), аппроксимирующей дифференциальное уравнение (1), предлагается с помощью спектрального метода. Возможность применения этого метода в данном случае обоснована в [7]. Представим решение в виде гармоники где λn - собственные числа оператора перехода; i - мнимая единица. Тогда разностная схема (2) примет вид (3) Произведя необходимые преобразования, выразим значение λ из (3): (4) В [8] показано, что для того, чтобы разностная схема (2) была устойчива, достаточно, чтобы все собственные числа оператора перехода удовлетворяли условию |λ| ≤ 1+cΔt, где c > 0 - независимая константа. Так как по условию теоремы константа c > 0, а шаг по времени Δt > 0, то очевидно, что cΔt > 0. Тогда неравенство |λ| ≤ 1 будет более строгим, чем |λ| ≤ 1+cΔt. Для упрощения практических расчетов в дальнейшем воспользуемся неравенством |λ| ≤ 1. (5) Воспользовавшись равенством (4), неравенством (5) и проведя необходимые преобразования, можем записать: (6) Рассмотрим некоторые частные случаи. Случай 1. Предположим самый упрощенный вариант, что ветра нет. Это означает, что проекции вектора средней скорости перемещения примеси на оси абсцисс, ординат и аппликат равны 0 (u = 0; v = 0; w = 0). Тогда неравенство (6) примет вид (7) Неравенство (7) является условием устойчивости конечно-разностной схемы (2) при условии отсутствия ветра. Случай 2. Предположим, что проекция вектора средней скорости перемещения примеси на оси абсцисс не равна 0, а проекции на оси ординат и аппликат равны 0 (u ≠ 0; v = 0; w = 0). Тогда можем полагать, что равенство (4) примет вид Таким образом, для устойчивости разностной схемы (2) требуется, чтобы собственные числа оператора перехода были расположены внутри или на границе круга радиусом 1, центр которого находится в начале координат комплексной плоскости. Проведя необходимые преобразования, получаем: (8) Неравенство (8) и является условием устойчивости конечно-разностной схемы (2) для случая, когда одна из проекций вектора скорости отлична от 0, а остальные равны 0. Случай 3. Предположим, что проекции вектора средней скорости перемещения примеси на оси абсцисс, ординат и аппликат не равны 0 (u ≠ 0; v ≠ 0; w ≠ 0). Тогда равенство (4) преобразуется к виду Для устойчивости разностной схемы (2), аналогично случаю 2, требуется, чтобы собственные числа оператора перехода были расположены внутри или на границе круга радиусом 1, центр которого находится в начале координат комплексной плоскости. Проведя необходимые преобразования, приходим к неравенству (9) Неравенство (9) и является условием устойчивости конечно-разностной схемы (2) для случая, когда проекции вектора средней скорости перемещения примеси на оси абсцисс, ординат и аппликат не равны 0. Таким образом, разностная схема (2) является условно устойчивой, в общем виде условие ее устойчивости получено в результате рассмотрения случая 3. Программная реализация Разработанный программный продукт был применен при проведении серии вычислительных экспериментов по прогнозированию и оценке негативных последствий, наступивших в результате аварии на железной дороге при транспортировке химически опасного груза. В качестве примера рассматривается возможная авария на железнодорожной станции с выбросом винилхлорида в атмосферу. Был произведен расчет изменения концентрации поллютанта во времени в заданных точках контролируемой территории, времени подхода шлейфа атмосферного загрязнителя к зонам жилой застройки и интенсивности их поражения для различных метеоусловий и размеров эмиссии. В качестве иллюстрации представлены результаты эксперимента в сечении z = 3 м при средних скоростях ветра u = 5 м/с, v = 2 м/c и интенсивности выброса Q = 1 кг/с через временные интервалы в 10, 60 и 90 с от начала аварии соответственно (рис.). Зона загрязнения атмосферы в случае гипотетической аварии: t = 10 с (а); t = 60 с (б) Зона загрязнения атмосферы в случае гипотетической аварии (окончание): t = 90 с (в) Следует отметить, что при подобной аварии скорость покрытия поллютантом зон жилой застройки настолько велика, что за такой короткий промежуток времени не удастся не только эвакуировать население, но даже предупредить о возникновении опасной ситуации. Предлагаемый метод позволяет не только оперативно прогнозировать динамику распространения опасного вещества в случае возникновения реальной аварии, но и моделировать последствия различных опасных ситуаций для разработки необходимых мер безопасности населения и окружающей среды. Выводы Разработана математическая модель, позволяющая прогнозировать интенсивность и направление распространения химически опасного вещества, выброшенного в атмосферу в результате аварийной ситуации. Построенная разностная схема позволяет быстро моделировать динамику распространения поллютанта в условиях определенной ограниченности начальных данных, характерных для данного класса задач. Разработанный программный продукт позволяет не только производить необходимые расчеты, но и предоставляет возможность визуально оценить распространение шлейфа атмосферного поллютанта по контролируемой территории. Согласно представленным результатам вычислений, при аварии с выбросом опасного химического вещества на некоторых участках железной дороги может произойти загрязнение зон жилой застройки, расположенных в непосредственной близости от железнодорожных путей.