Сопряженные задачи переноса и диффузии в проблеме оценки и прогноза состояния окружающей среды

МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ

Федеральное государственное бюджетное образовательное учреждение

высшего профессионального образования

“КУБАНСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ“

(ФГБОУ ВПО “КубГУ“)

Кафедра математического моделирования

ВЫПУСКНАЯ КВАЛИФИКАЦИОННАЯ (ДИПЛОМНАЯ)

РАБОТА

Сопряженные задачи переноса и диффузии в проблеме оценки и прогноза состояния окружающей среды

Работу выполнил М.А. Кучумов

Факультет ФКТиПМ

Специальность/направление 010501 – Прикладная математика и информатика

Научный руководитель, д-р физ. – мат наук, профессор А.В. Павлова

Нормоконтролер

канд. физ. – мат наук, доцент М.С. Капустин

Краснодар 2013

Реферат

Дипломная работа __ с., 10 рис., 13 источников, 2 приложения.

ЗАГРЯЗНЯЮЩЕЕ ВЕЩЕСТВО, ДИФФУЗИЯ, ПЕРЕНОС, ПРЯМАЯ и СОПРЯЖЕННАЯ ЗАДАЧИ, задачи регулирования загрязнения

Объектом исследования является процесс распространения загрязнителей.

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

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

Исследования проводились с помощью аналитических и разностных методов.

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

Содержание

• Введение
• 1. Модели оценки загрязнения атмосферы на основе уравнения турбулентной диффузии
• 1.1 Постановка прямой задачи диффузии и переноса загрязнений
• 2. Методы решения прямой и сопряженной задач
• 3. Описание программной реализации решения прямой и сопряженной задач
• 3.1 Реализация алгоритма метода покомпонентного расщепления для прямой задачи
• 3.2 Реализация алгоритма решения сопряженной задачи
• 3.3 Программный интерфейс
• 4. Задачи оптимизации и регулировки мощности источников
• 4.1 Определение функционала чувствительности
• 4.2 Задача оптимального размещения источников промышленных выбросов
• 4.3 Задачи регулировки мощности источников
• 4.3.1 Оптимизация с помощью основных уравнений
• 4.3.2 Оптимизация с помощью сопряженной задачи
• Заключение
• Список использованных источников
• Приложения

Введение

Атмосфера является составной частью биосферы и представляет собой оболочку Земли, вращающуюся вместе с ней как единое целое. Эта оболочка условно разделена на несколько слое, каждый из которых имеет своё название и характерные физико-химические особенности. Условно принято атмосферу делить на две большие составные части верхнюю и нижнюю. Наибольший интерес с точки зрения экологии представляет нижняя часть атмосферы главным образом тропосфера, поскольку в ней происходят основные метеорологические явления, влияющие на загрязнение атмосферного воздуха.

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

В связи с интенсивным функционированием всех отраслей промышленности, энергетики, транспорта увеличение массы загрязняющих выбросов оказывает всё более серьёзное влияние на окружающую среду. Биосфера обладает значительной устойчивостью по отношению к загрязняющим примесям. Эта устойчивость основана на естественной способности различных компонентов природной среды к самоочищению. Вместе с тем способность “усваивать” без заметного ущерба различного рода примеси в биосфере не безгранична. Перед обществом возникает сложная проблема оптимальных соотношений между объёмами производства и необходимостью обеспечивать должный уровень чистоты окружающей среды. Сложность этой проблемы определяется не только трудностями, возникающими при решении частных задач, связанных с охраной окружающей среды. Изучение этой проблемы в целом требует, прежде всего, интеграции исследований целого ряда научных дисциплин: биологии и географии, экономики и медицины, химии и юриспруденции, физики атмосферы и математики, космических исследований и технологий.

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

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

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

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

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

Этапы эволюции методов исследования отображены, например, в работах М.Е. Берлянда [1], Л.Н. Бызовой [2] др. Применение моделей атмосферной диффузии к проблемам охраны окружающей среды получило значительное развитие в работах Г.И. Марчука [3,4].

Целью дипломной работы является построение математической модели, описывающей процесс распространения пассивных загрязняющих веществ от сосредоточенных источников, использование аппарата сопряженных задач для определения наиболее безопасных зон размещения объектов, загрязняющих атмосферу, расчет распределения концентрации примеси и заданных функционалов. В качестве исходного уравнения, было взято уравнение турбулентной диффузии, применяемое Марчуком Г.И. и его методика построения и решения сопряженных задач, применяемые в работах по математическому моделированию проблем, связанных с охраной окружающей среды [3-5].

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

математическая модель загрязняющее вещество

1. Модели оценки загрязнения атмосферы на основе уравнения турбулентной диффузии

Численная модель распространения примесей в пограничном слое атмосферы основана на полуэмпирическом уравнении турбулентной диффузии [3,4].

1.1 Постановка прямой задачи диффузии и переноса загрязнений

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

Рассматриваемое в работе уравнение, описывающее перенос загрязняющей примеси имеет вид

, (1.1)

где q – величина, характеризующая концентрацию ЗВ, – коэффициент концентрации, – вертикальный коэффициент диффузии, – горизонтальный коэффициент диффузии, дельта-функция Дирака моделирует сосредоточенный источник, постоянной мощности Q.

Рассматривается конечная область рассеяния примеси с границей S при следующих начальных и граничных условиях:

q=q₀ при t=0 – начальное распределение концентрации загрязнителя;

при y=0 – частичное поглощение/отражение на подстилающей поверхности ( – коэффициент, характеризующий взаимодействие примеси с подстилающей поверхностью);

при y=H – отсутствие вертикальных потоков на высоте Н;

q=0 при x=0, x=a – выход на фоновые значения концентрации на боковых границах.

Поскольку для нижних слоев атмосферы выполняется закон сохранения масс, к уравнениям, описывающим перенос и диффузию примеси, необходимо добавить уравнение неразрывности

.

Дискретные аппроксимации для рассматриваемых задач строятся на базе вариационного принципа в сочетании с методом расщепления [3,4,6,7]. При таком подходе осуществление технологически простой реализации дискретной модели на ЭВМ.

1.2 Сопряженная задача диффузии и переноса загрязнений

Проблема планирования размещения промышленных предприятий с учетом санитарно-допустимых норм загрязнения тесно связана с проблемами управления качеством атмосферы и охраны окружающей среды. В работах [3-5, 7-9] сформулирован ряд математических моделей для решения такого рода задач и изложены подходы к их решению.

При решении некоторых задач, в частности, связанных с возможным размещением предприятий с соблюдением санитарных норм загрязнения для экологически значимых зон, когда объектом изучения является не само поле концентрации, а некоторые функционалы от этого поля, удобно применять сопряженные уравнения [3,5,9].

Пусть рассматривается функционал

, (1.2)

где

. Если, например,

, , , (1.3)

то функционал J (q) представляет собой концентрацию q (x,y,t) в выделенной подобласти области D,

где , взвешенную с весом ₀. Область соответствует зоне, где производится оценка загрязнения атмосферы.

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

Для оценки функционала J (q) целесообразно реализовать подход, основанный на использовании сопряженных уравнений.

Следуя [5], сформулируем сопряженную задачу.

С этой целью уравнение (1) умножим на функцию q^*, которую будем называть сопряженной, и результат проинтегрируем по .

Причем в силу того, что у нас выполняется условие неразрывности

, где ,

внесем под знак соответствующей производной. Тогда

Рассмотрим первое слагаемое в левой части. Меняя порядок и интегрируя по t по частям, получим

Рассмотрим второе слагаемое в левой части. Аналогично проделанному выше меняем порядок интегрирования, но по частям проинтегрируем по .

Так же поступим с третьим и четвертым слагаемым левой части рассматриваемого равенства.

.

Сопряженное уравнение для (1.1) будет иметь вид

. (1.4)

Следует отметить, что сопряженная задача корректна при решении от к .

В качестве “начальных” условий рассмотрим при t=T q* (T,M) =q* (0,M) =0. Граничные условия примут вид

при 0 < x < a,

при 0 < y < H,

.

В случае если на боковых гранях параллелепипеда были заданы граничные условия вида =0, для сопряженной задачи на боковых гранях будут справедливы условия ^* = 0. Тогда

.

В рассматриваемом случае

при y=0, при y=H, q*=0 при x=0, x=a.

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

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

2. Методы решения прямой и сопряженной задач

Численное моделирование проблем распространения загрязняющих веществ проводится чаще всего на основе использования разностных схем.

Запишем поставленную в 1.1 задачу операторном виде [6], т.е.

, в области ,

где . Здесь,

, , в D при t=0.

Для решения этой задачи вводится сеточная область . Чтобы построить приближенное решение, необходимо заменить исходную дифференциальную задачу некоторой конечномерной дискретной задачей, обычно представляющей собой систему линейных алгебраических уравнений (СЛАУ) для компонентов вектора . Одной и той же дифференциальной задаче можно поставить в соответствие множество различных дискретных моделей, однако, далеко не все из них пригодны для практической реализации. Вопросы построения разностных схем подробно рассмотрены в [6,11].

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

. (2.1)

Соотношение (2.1) является системой обыкновенных дифференциальных уравнений для компонентов вектора . В дальнейшем индекс ставить не будем, полагая что (2.1) есть разностный аналог по пространственным переменным исходной задачи.

С учетом изложенного рассмотрим задачу Коши

, при . (2.2)

Рассмотрим простейшие методы аппроксимации задачи (2.2) по времени, полагая, что не зависит от времени. Одна из них – явная схема первого порядка аппроксимации на сетке

(2.3)

Неявная схема первого порядка аппроксимации имеет вид

(2.4)

Это схемы первого порядка аппроксимации (в предположении, что существуют вторые производные по времени функции ).

Схема Кранка – Николсон имеет вид [6]:

(2.5)

Эта схема аппроксимирует исходную задачу со вторым порядком по времени.

При этом схема (2.3) будет устойчива при выполнении определенного условия, например, если – симметричная положительно определенная матрица с собственными числами из интервала [с,b], а удовлетворяет соотношениям .

Рассмотрим метод расщепления, применяемых при решении двумерных задач, описывающих процессы различной природы. Рассматриваемое эволюционное уравнение имеет вид (2.1) в . Оператор не зависит от времени и представим в виде при условии, что . Будем полагать, что эта задача уже редуцирована к разностному виду. Рассмотрим так называемый метод стабилизации. Для этого рассмотрим разностную схему решения (2.1) в предположении f=0:

. (2.6)

Эта схема аппроксимирует исходную задачу со вторым порядком аппроксимации по .

С помощью ряда преобразований (2.6) приводится к виду

(2.7)

Отсюда видно, что (2.7) при достаточной гладкости решения совпадает со схемой Кранка – Николсона, т.е. имеет второй порядок аппроксимации по .

Эта разностная схема допускает удобную компьютерную реализацию.

В случае неоднородной задачи

при t=0. (2.8)

В этом случае задача запишется следующим образом:

(2.9) где

. (2.10)

При условии (2.10) схема (2.8) будет обладать вторым порядком аппроксимации по .

Рассмотрим первый этап, который соответствует переносу. Разобьем весь промежуток [0,T] на элементарные. Переносу соответствует следующий оператор:

.

Соответствующий разностный аналог имеют вид

,

где – шаги сетки вдоль соответствующих осей.

Рассмотрим следующий этап, соответствующий диффузии. А именно будем рассматривать оператор

.

Разностный аналог оператора будет иметь вид

.

За при этом принимается решение задачи на предыдущем этапе (когда рассматривался только перенос).

На третьем этапе рассматривается трансформация примеси, при этом рассматриваются лишь трансформация компонентов примеси, связанная с разложением их на атмосферные составляющие.

,

где – значение для вектора концентраций полученное на предыдущем этапе.

В результате получается схема второго порядка аппроксимации по и по .

Структура алгоритмов для решения сопряженных задач полностью определяется структурой алгоритмов решения основной задачи.

3. Описание программной реализации решения прямой и сопряженной задач

Для решения поставленных задач была написана программа, в которой был реализован методы: покомпонентного расщепления и итерационный. Для написания программы была использована среда Borland Delphi 7. Этот выбор обусловлен тем, что язык программирования Delphi представляет широкий инструментарий для создания удобного пользовательского интерфейса, визуализации и численного решения задач.

3.1 Реализация алгоритма метода покомпонентного расщепления для прямой задачи

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

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

Все конечно-разностные аналоги исходных операторов являются положительно полуопределёнными и аппроксимирующими со вторым порядком точности. Производные входящие в граничные условия аппроксимируются с первым порядком точности.

Система алгебраических уравнений решается при помощи метода прогонки [6,11] (так как матрица этой системы является блочнодиагональной с трехдиагональными блоками, первое и последнее уравнение в блоке аппроксимируют соответствующие граничные условия). Потом происходит удовлетворение граничных условий, неучтенных в этой системе. Промежуточный результат записывается в трехмерный массив fipr, затем, используя получившиеся данные, переходим к следующей системам и проделываем аналогичные действия. После того как завершился шаг покомпонентного расщепления, результаты записываются в массив fi и управление передаётся на начало цикла. После завершения процесса вычислений получившиеся данные записываются в массив fi. Для построения графика необходимо ввести интересующую высоту y и нажать кнопку “Прямая” в программном окне.

3.2 Реализация алгоритма решения сопряженной задачи

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

Для построения графика необходимо ввести интересующую высоту y и нажать кнопку “Сопряженная” в программном окне.

3.3 Программный интерфейс

Интерфейс программы представлен на рисунке 1. Форма приложения содержит поля для ввода данных, две вкладки с компонентами вывода графиков функций. На первой вкладке выводится график решения прямой задачи, на второй – график решения сопряженной задачи.

Рисунок 1 – Вид окна приложения

В программе можно изменить следующие параметры: L – правая граница по оси X; T – задает время; const f, const x, const y – задает мощность и координаты источника; n – количество точек разбиения по x; U,V – скорости перемещения примеси; my,ny – задают коэффициенты диффузии; k1, k2 – коэффициенты третьего граничного условия при y=0; tau – шаг по времени. После нажатия кнопки “Расчет” происходит расчет решения и появляются: компоненты графиков; кнопка “Прямая”, строящая график решения прямой задачи на первой вкладке; поле изменения высоты y и кнопка “Сопряженная”, строящая график решение сопряженной задачи на второй вкладке.

4. Задачи оптимизации и регулировки мощности источников

4.1 Определение функционала чувствительности

Состояние атмосферы находится в непосредственной зависимости от уровня развития промышленности в регионе, действие любого источника выброса загрязняющего вещества представляет некоторую опасность для окружающей среды. При решении ряда оптимизационных задач, связанных с возможным размещением новых предприятий или регламентации работы существующих с целью соблюдения санитарных норм, объектом изучения являются некоторые функционалы от поля концентрации загрязняющего вещества. К числу изучаемых функционалов относятся, например, полное количество примесей, осевших в конкретном регионе, среднегодовая концентрация различных субстанций, экономический ущерб, наносимый окружающей среде и др. В таких задачах использование сопряженных функций является весьма эффективным [3,5,9].

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

Пусть рассматривается функционал (1.2)

, где .

Если, например, определяется соотношением (1.3), то функционал J представляет собой суммарную концентрацию примеси в выделенной подобласти области D, взвешенную с весом ₀ или коллективную экспозицию, которой будет подвержено население региона в результате эмиссии загрязнителя источником. Область G соответствует охраняемой зоне. Задача сводится к оценке функционала J, определяемого на множестве функций состояния, удовлетворяющих исходной задаче.

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

. (4.1)

Будем рассматривать модельный пример в области . Пусть в рассматриваемой области функционируют n предприятий. Они характеризуются координатами r_i (x_i,y_i), временем и мощностью выброса, соответственно _i и Q_i. Пусть экологически значимые зоны характеризуются координатами в этой области r_T (x_T,y_T). Рассматривается их состояние в момент Т, 0 _i < T, i=1,2,…,n [10]. Уравнения переноса примесей в атмосфере (1.1) позволяют получать распределение осредненных значений концентраций примесей в заданной области. Распространение примеси описывается уравнением (1.1), которое в случае рассмотрения плоской задачи при наличии n импульсных сосредоточенных источников примет вид

(4.2)

при , . Здесь принято

Пусть в качестве искомого функционала рассматриваем уровень загрязнения в точке r_T в момент времени Т. Возникает необходимость рассмотрения сопряженных уравнений диффузии.

Сопряженное уравнение для (4.2) будет иметь вид

,

при ,

Решения этих задач может быть представлено аналитически в следующем виде [10]:

(4.3)

(4.4)

(4.5)

Для решения модельных задач несложные вычисления могут быть проведены с помощью Microsoft Excel.

Введем сеточную область G , с шагом 1.

Для узлов сетки нетрудно получить следующие значения: – величина, характеризующая концентрацию ЗВ, q* – чувствительность точки r к загрязнению, J – уровень загрязнения в точке r_T в момент Т.

Предположим, существует два источника выброса со следующими характеристиками:

1) мощность выброса Q_i=1, время выброса _i=30,6, r_i (2,66, 2,66);

2) мощность выброса Q_i=2, время выброса _i=33,3, r_i (10,13, 10,5).

Горизонтальный и вертикальный коэффициенты диффузии примем равными 0,1 (==0,1), ₁=0,02, ₂=0,5; экологически значимую зону, располагающуюся в точке r_T= (15, 15), будем рассматривать в момент времени Т=50. Наблюдения будем проводить в t=40.

Рисунки 2, 3 иллюстрируют распределение q и линии уровня функции чувствительности q* для области G.

Рисунок 2 – Диаграмма концентраций ЗВ при двух источниках выбросов

Из рисунка 1 видно, что наибольшая концентрация загрязняющих веществ достигается в точках, практически соответствующих координатам выбросов (скорость ветра мала). Причем, чем больше мощность источника, тем выше концентрация, что соответствует реальной картине.

Рисунок 2 – Диаграмма чувствительности точек области к загрязнению

Рисунок 2 показывает, что наиболее предрасположены к загрязнению точки охраняемой области. При этом уровень загрязнения в точке r_Т (х_Т,y_Т) в момент времени Т функционал J=8,96.

4.2 Задача оптимального размещения источников промышленных выбросов

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

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

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

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

Пусть в рассматриваемой области G имеются охраняемые зоны G_k. При этом необходимо разместить источник загрязнения заданной мощности Q так, чтобы среднегодовое значение концентрации вредной примеси в охраняемых зонах не превысило предельно допустимое значение концентрации (ПДК) С=const.

В качестве оцениваемых функционалов выберем осредненное по времени и зонам G_k значение концентрации, где, согласно введенным обозначениям, q – решение прямой задачи для уравнения (1.1)

.

Но рассматривать проблему размещения промышленных предприятий посредством решения прямой задачи для уравнения (1.1) малоэффективно, т.к. приходится решать прямую задачу раз (n_i – число шагов сетки вдоль соответствующие оси, i=1,2). В то же время использование сопряженного уравнения (1.4) значительно уменьшает трудности в нахождении решения, поскольку для каждой из зон G_k требует только однократного решения соответствующей сопряженной задачи.

Двойственный функционал будет иметь вид [5]

.

Пусть , при , где M – точка размещения источника загрязнения. Остается определить такие точки M, в которых значение функционала меньше ПДК. Эти точки и образуют область возможного размещения источников выбросов.

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

Пусть источник вредной примеси нужно разместить в точке r₀= (x₀,y₀). В случае если рассматриваются точечные источники, функционал (1.3) примет вид

, (4.6)

где – решение сопряженной задачи, которое находится из уравнения (1.4).

Использование двойственного представления функционала позволяет определить область для каждого k=1,2,…m.

Пересечение охраняемых областей непусто. Точка r₀= (x₀,y₀,z₀) является искомой, если удовлетворяет условию, где .

4.3 Задачи регулировки мощности источников

Некоторые способы борьбы с загрязнением промышленными предприятиями окружающей среды приводят к закрытию производства, которое в свою очередь, ведет к сокращению рабочих мест. Таким образом, с одной стороны экономика и экология могут понести значительный ущерб из-за наличия вредных примесей в окружающей среде, с другой – зачастую невозможно отказаться от предприятий, по вине которых происходит загрязнение окружающей среды. Возникает проблема нахождения оптимума между полным закрытием “грязных” предприятий и бесконтрольным загрязнением окружающей среды. Для решения такого класса задач могут использоваться модели оптимизации регулирования мощности источников с использованием аппарата линейного и нелинейного программирования [12,13].

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

Пусть в заданном регионе G с границей S в точках r_i (i=1,2,…,n) расположены n промышленных объектов A_i, c мощностями выбросов Q_i соответственно (i=1,2,…,n). Состав выбросов для простоты считается одинаковым.

В области G выделим m экологических зон G_k (k=l, 2,., m), для каждой из которых заданы предельно допустимые концентрации выпавшего за интервал времени [0, Т] аэрозоля. В результате приходим к следующей математической постановке задачи.

Уравнение рассеяния загрязняющей примеси от n индустриальных объектов (сосредоточенных источников постоянной мощности) имеет вид

. (4.7)

Пусть заданы следующие условия:

при y=0; (4.8) при y=H.

Считая задачу (4.7), (4.8) климатически периодической (с периодом по времени, равным году), получим начальные условия

q (T,x,y) = q (0,x,y). (4.9)

Рассмотрим функционал вида (1.2), характеризующий санитарную дозу аэрозоля, выпавшего на подстилающую поверхность (y=0) в области экологической зоны G_k.

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

, (4.10)

где – исходная мощность выбросов, Q_i – планируемая мощность выбросов, – коэффициент, определяющий капитальные вложения технологию, обеспечивающую выпуск того же объема продукции при уменьшении выбросов (в расчете на единицу мощности выбросов). Тогда функционал I представляет полные затраты, необходимые для улучшения технологии всех предприятий А_i при переходе от выбросов планируемым выбросам Q_i. В результате приходим к задаче (4.7) – (4.9) о нахождении таких выбросов Q_i, чтобы выполнялись условия

, (4.11) J_k?c_k, k=1,2,…,m.

Задачу (4.7) – (4.9), (4.11) можно свести к задаче линейного программирования [12]. При этом возможны два различных подхода, один из которых реализуется с помощью основных уравнений, другой – с помощью сопряженных.

4.3.1 Оптимизация с помощью основных уравнений

Решение задачи (4.7) – (4.9) представим в виде суперпозиции решений элементарных задач. Пусть

, (4.12)

где – решение прямой задачи.

Подставляя (4.12) в (2.1), получаем

, (4.13)

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

,, i=1,2,…,n; k = l, 2,., m.

Здесь и – известные константы. С учетом введенных обозначений переходя к , получим задачу линейного программирования по отысканию оптимального набора q_i на основе решения задачи

, k = l, 2,., m, q_i?0, i = 1,2,…,n,

где .

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

4.3.2 Оптимизация с помощью сопряженной задачи

С учетом (4.1) функционал сопряженной задачи будет иметь вид

.

Таким образом, аналогично случаю основной задачи, приходим к оптимизационной задаче для сопряженных уравнений:

,, k = l, 2,., m,q_i?0, i = 1,2,…,n,

где

Таким образом, снова приходим к задаче линейного программирования.

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

Заключение

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

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

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

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

Список использованных источников

Берлянд, М.Е. Современные проблемы атмосферной диффузии. – Ленинград: Гидрометеоиздат, 1975. – 448 с.

Бызова, Н.Л. Экспериментальные исследования атмосферной диффузии и расчет распространения примеси / Н.Л. Бызова, Е.К. Гаргер, В.Н. Иванов – Ленинград: Гидрометеоиздат, 1991. – 279 с.

Марчук, Г.И. Математическое моделирование в проблеме охраны окружающей среды. – М.: Наука, 1982. – 320 с.

Марчук, Г.И. Численное решение задач динамики атмосферы и океана. – Ленинград: Гидрометеоиздат, 1974. – 303 с.

Марчук Г.И. Сопряженные уравнения, Москва: Институт вычислительной математики РАН, 2001. – 241 с.

Марчук Г.И. Методы вычислительной математики. – М.: Наука, 1989. – 608 с.

Агошков В.И. Методы решения задач математической физики / В.И. Агошков, П.Б. Дубовский, В.П. Шутяев. – М.: ФИЗМАТЛИТ, 2002. – 303 с.

Пененко В.В., Модели и методы для задач охраны окружающей среды / В.В. Пененко, А.Е. Алоян. – Новосибирск: Наука, 1985. – 254 с.

Алоян, А.Е. Моделирование динамики и кинетики газовых примесей и аэрозолей в атмосфере. – М.: Наука, 2008. – 415 с.

Петросян, Л.А. Математические модели в экологии / Л.А. Петросян, В.В. Захаров. – СПб: Изд. СпбГУ, 1997. – 256 с.

Самарский, А.А. Численные методы решения задач конвекции-диффузии / А.А. Самарский, П.Н. Вабищевич. – М.: URSS: Либроком, 2009. – 246 с.

Винтер Г. Модель экологического менеджмента. – Минск: Технопринт, 2001. – 320 с.

Пахомова Н. Экологический менеджмент / Н. Пахомова, А. Эндрес, К. Рихтер. – СПб.: Питер, 2003. – 544 с.

Приложения

Приложение А

Примеры решения прямых и сопряженных задач

Было решено несколько прямых и сопряженных задач с постоянным точечным источником в области D={0<x<2, 0<y<2}. В первой задаче источник располагался в точке М (0.5; 0.5). Так же полагалось: U=0.2, V=0, k1=1, k2=0.5, my=ny=0.5. Шаг по пространственной переменной выбирался равный 0,1, а по временной 0,01.

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

Рисунок А.1 ? Распределение концентрации при y=0

При у=0.6, т.е. выше источника распределение концентрации следующее:

Рисунок А.2 ? Распределение концентрации при y=0.6

Второй пример. В области D={0<x<4, 0<y<2}. Источник располагается в точке М (0,1; 0,8). Так же полагалось: U=V=0.5, k1=0.7, k2=0.4, my=0.6, ny=0.3. Шаг по пространственной переменной выбирался равный 0,2, а по временной 0,01.

При решении этой задачи методом покомпонентного расщепления были получены следующие результаты при Т=3. При y=0.4 (ниже источника) распределение концентрации выглядит следующим образом:

Рисунок А.3 ? Распределение концентрации при y=0.4

При у=0.6, т.е. выше источника распределение концентрации следующее:

Рисунок А.4 ? Распределение концентрации при y=0.6

Третий пример. В области D={0<x<2, 0<y<2}. Источник располагается в точке М (1; 0,5). Так же полагалось: U=-0.25, V=0.25, k1=1, k2=0.5, my=1, ny=0.5. Шаг по пространственной переменной выбирался равный 0.1, а по временной 0.01.

При решении спряженной задачи были получены следующие результаты при Т=2. При y=0,1 и при y=1.

Рисунок А.5 ? Функция чувствительности при y=0.1

Рисунок А.6 ? Функция чувствительности при y=1

Четвертый пример. В области D={0<x<2, 0<y<2}. Источник располагается в точке М (1; 0,5). Так же полагалось: U=0, V=0 (штиль), k1=10, k2=5, my=2, ny=1. Шаг по пространственной переменной выбирался равный 0.1, а по временной 0.01.

При решении этой задачи были получены следующие результаты при Т=1. При y=0,1 и при y=0,5.

Рисунок А.7 ? Распределение концентрации при y=0.1

Рисунок А.8 ? Распределение концентрации на уровне источника

Приложение Б

ЛИСТИНГ ПРОГРАММЫ

const

n1=400;

delta=0.001;

Nr=201;

type

vect=array [0. n1] of Extended;

matr=array [0. n1,0. n1] of Extended;

vectr=array [0. nr] of Extended;

var

L,M,T,k1,k2,U,V,my,ny,h,tao,cons,tc: Extended;

i,j,k,n,z: integer;

fi,fipr,fi1: matr;

a1,a2,a3,b,pr: vect;

yr,xr,ydr,fr,farr,ftr,ar,br,cr,dr,er,ffr: vectr;

function func (o,s: integer): real;

begin

if (o=strtoint (frmmain. LabeledEdit14. Text)) and

(s=strtoint (frmmain. LabeledEdit15. Text)) then

func: =cons*cons

else func: =0;

end;

procedure Progonka (var a1,a2,a3,b,pr: vect);

var

pr1,pr2: vect;

i1: integer;

begin

pr1 [n-1]: =-a1 [n] /a2 [n];

pr2 [n-1]: =b [n] /a2 [n];

for i1: =n-1 downto 1 do

begin

pr1 [i1-1]: =-a1 [i1] / (pr1 [i1] *a3 [i1] +a2 [i1]);

pr2 [i1-1]: = (b [i1] – pr2 [i1] *a3 [i1]) / (pr1 [i1] *a3 [i1] +a2 [i1]);

end;

pr [0]: = (b [0] – a3 [0] *pr2 [0]) / (a2 [0] +a3 [0] *pr1 [0]);

for i1: =1 to n do

pr [i1]: =pr1 [i1-1] *pr [i1-1] +pr2 [i1-1];

end;

procedure Prog5 (var nr: integer; ar,br,cr,dr,er,fr,yr: vectr; itask: integer);

var

i: integer;

begin

if itask=1 then

begin

dr [1]: =dr [1] /cr [1];

er [1]: =er [1] /cr [1];

cr [2]: =cr [2] – dr [1] *br [2];

dr [2]: = (dr [2] – er [1] *br [2]) /cr [2];

er [2]: =er [2] /cr [2];

for i: =3 to nr do

begin

cr [i]: =er [i-2] *ar [i] +dr [i-1] * (dr [i-2] *ar [i] – br [i]);

dr [i]: = (dr [i] +er [i-1] * (dr [i-2] *ar [i] – br [i])) /cr [i];

er [i]: =er [i] /cr [i];

end;

itask: =2;

end;

fr [1]: =fr [1] /cr [1];

fr [2]: = (fr [2] +fr [1] *br [2]) /cr [2];

for i: =3 to n do

fr [i]: = (fr [i] – fr [i-2] *ar [i] – fr [i-1] * (dr [i-2] *ar [i] – br [i])) /cr [i];

yr [nr]: =fr [nr];

yr [nr-1]: =dr [nr-1] *yr [nr] +fr [nr-1];

for i: =nr-2 downto 1 do

yr [i]: =dr [i] *yr [i+1] – er [i] *yr [i+2] +fr [i];

end;

procedure Prog3 (var nr: integer; ar,cr,br,fr,yr: vectr; itask: integer);

var

i: integer;

begin

if itask=1 then

begin

br [1]: =br [1] /cr [1];

for i: =2 to nr do

begin

cr [i]: =cr [i] – br [i-1] *ar [i];

br [i]: =br [i] /cr [i];

end;

itask: =2;

end;

fr [1]: =fr [1] /cr [1];

for i: =2 to n do

fr [i]: = (fr [i] +fr [i-1] *ar [i]) /cr [i];

yr [nr]: =fr [nr];

for i: =nr-1 downto 1 do

yr [i]: =br [i] *yr [i+1] +fr [i];

end;

procedure TfrmMain. cmdRunClick (Sender: TObject);

begin

L: = StrToFloat (LabeledEdit1. Text);

M: = StrToFloat (LabeledEdit2. Text);

T: = StrToFloat (LabeledEdit3. Text);

k1: = StrToFloat (LabeledEdit4. Text);

k2: = StrToFloat (LabeledEdit5. Text);

U: = StrToFloat (LabeledEdit6. Text);

V: = StrToFloat (LabeledEdit7. Text);

my: = StrToFloat (LabeledEdit8. Text);

ny: = StrToFloat (LabeledEdit9. Text);

n: = StrToInt (LabeledEdit10. Text);

cons: = StrToFloat (LabeledEdit11. Text);

tao: = StrToFloat (LabeledEdit16. Text);

h: =L/n cmdRun. Enabled: =false;

for i: =1 to n-1 for j: =1 to n-1 do

fi [i,j]: =0;

for j: =1 to n-1 do begin

fi [0,j]: =0;

fi [n,j]: =0;

end;

for i: =1 to n-1 do

begin

fi [i,n]: =0;

fi [i,0]: = ( (k1*fi [i,1]) / (k1+k2*h));

end;

tc: =0.0;

fi1: =fi;

Z: =0;

while (abs (tc-T) >=tao) and (z<=10000000) do

begin

z: =z+1;

tc: =tc+tao;

fi: =fi1;

for j: =1 to n-1 do

begin

a1 [0]: =0.0;

a2 [0]: =1.0;

a3 [0]: =-1.0;

b [0]: =0.0;

a1 [n]: =1.0;

a2 [n]: =-1.0;

a3 [n]: =0.0;

b [n]: =0.0;

for k: =1 to n-1 do

begin

a1 [k]: = ( ( (-u*tao) / (4*h)) – (my*tao) / (2*h*h));

a2 [k]: = (1.0+ (my*tao) / (h*h));

a3 [k]: = ( ( (u*tao) / (4*h)) – (my*tao) / (2*h*h));

b [k]: = (fi [k,j] – ( (tao*u* (fi [k+1,j] – fi [k-1,j])) / (4*h) +

(tao*my* (fi [k+1,j] – 2*fi [k,j] +fi [k-1,j])) / (2*h*h)));

Application. ProcessMessages;

end;

Progonka (a1,a2,a3,b,pr);

for k: =0 to n do

fipr [k,j]: =pr [k];

end;

{for k: =1 to n-1 do

begin

fipr [0,k]: =0;

fipr [n,k]: =0;

end;

}

for k: =1 to n-1 do

begin

fipr [k,n]: =0;

fipr [k,0]: = ( (k1*fipr [k,1]) / (k1+k2*h));

end;

fi: =fipr;

for i: =1 to n-1 do

begin

a1 [0]: =0.0;

a2 [0]: =-k1/h-k2;

a3 [0]: =k1/h;

b [0]: =0.0;

a1 [n]: =0;

a2 [n]: =1.0;

a3 [n]: =0.0;

b [n]: =0.0;

for k: =1 to n-1 do

begin

a1 [k]: = ( ( (-v*tao) / (4*h)) – (ny*tao) / (2*h*h));

a2 [k]: = (1.0+ (ny*tao) / (h*h));

a3 [k]: = ( ( (v*tao) / (4*h)) – (ny*tao) / (2*h*h));

b [k]: = (fi [i,k] – ( (tao*v* (fi [i,k+1] – fi [i,k-1])) / (4*h) +

(tao*ny* (fi [i,k+1] – 2*fi [i,k] +fi [i,k-1])) / (2*h*h))) +tao*func (i,k);

end;

Progonka (a1,a2,a3,b,pr);

for k: =0 to n do

fipr [i,k]: =pr [k];

end;

fi: =fipr;

fi1: =fipr;

end;

chart1. Visible: =true;

cmdGraf. Visible: =true;

cmdRun. Enabled: =true;

labeledEdit12. Visible: =true;

label1. Caption: =IntToStr (z);

z: =0;

end;

end.