Размещено на http://www.allbest.ru/

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

Федеральное государственное автономное образовательное учреждение высшего образования

НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ

ТОМСКИЙ ПОЛИТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ

Институт неразрушающего контроля

Кафедра экологии и безопасности жизнедеятельности

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



СОДЕРЖАНИЕ

ЗАДАНИЕ

ВВЕДЕНИЕ

1. ТЕОРЕТИЧЕСКАЯ ЧАСТЬ

1.1 Состояние атмосферного воздуха

1.2 Энергетический источник

1.3 Преобразование вредных выбросов ТЭЦ в атмосферном воздухе

1.4 Математические модели распространения примеси

1.4.1 Гауссова модель

1.4.2. ANSYS/FLOTRAN

2. ПРАКТИЧЕСКАЯ ЧАСТЬ

2.1 Постановка физической задачи

2.2 Постановка математической задачи

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

2.4 Анализ полученных результатов

ЗАКЛЮЧЕНИЕ

СПИСОК ИСПОЛЬЗОВАННЫХ ИСТОЧНИКОВ

ПРИЛОЖЕНИЕ А



ЗАДАНИЕ

на выполнение курсовой работы

Студенту гр. 1Е31 Вилесову Максиму Вадимовичу

(ФИО)

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

2. Срок сдачи студентом готовой работы: 22.04.17г.

3. Исходные данные к работе: Источник - труба котельной, топливо - уголь, поллютант - диоксид серы (SO₂), (h_об)высота области 150м, (L_об) длина области 300м,h (высота трубы) = 100м, d (диаметр устья) = 6м, время расчета - 5мин, изменения скорости ветра от 1м/с до 5 м/с.

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

5. Дата выдачи задания на выполнение курсовой работы: 16.03.17г.

Руководитель __________________ Перминов В.А.

(подпись, дата)

Задание принял к исполнению _______________ Вилесов М.В.

(подпись, дата)



ВВЕДЕНИЕ

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

Тепловые электростанции работают на относительно дешевом органическом топливе - угле и мазуте, это невосполнимые природные ресурсы. Сегодня основными энергетическими ресурсами в мире являются уголь(40%), нефть (27%) и газ (21%).

При сжигании топлива на ТЭС образуются продукты сгорания, в которых содержатся: летучая зола, частички несгоревшего пылевидного топлива, серный и сернистый ангидрид, оксид азота, газообразные продукты неполного сгорания. При зажигании мазута образуются соединения ванадия, кокс, соли натрия, частицы сажи. В золе некоторых видов топлива присутствует мышьяк, свободный диоксид кальция, свободный диоксид кремния. [1]

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

Задачи данной курсовой работы:

1. Рассмотрение теоретических аспектов загрязнения атмосферы от энергетических источников

2. Составление программного кода для двухмерной задачи в программе математического моделирования MATLAB

3. Представление результатов в виде графиков, проверка и анализ полученных данных

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



1. ТЕОРЕТИЧЕСКАЯ ЧАСТЬ

1.1 Состояние атмосферного воздуха

Перенос примесей в верхние слои атмосферы определяется характером распределения температуры в атмосфере с высотой - устойчивостью атмосферы, степенью изменения температуры с высотой, то есть вертикальным градиентом температуры на единицу расстояния, обычно на сто метров. Состояние атмосферы бывает равновесным, устойчивым и неустойчивым. Степень устойчивости атмосферы определяет поведение воздушной частицы, выведенной из первоначального положения в выше или ниже лежащий слой атмосферы. Частицы теплого воздуха поднимаются вверх, холодного - опускаются вниз. Обычно в атмосфере происходит падение температуры с высотой. Если градиент температуры в сухой атмосфере равен 1 градус на сто метров, то воздух на любой высоте будет находиться в равновесии. Каждая воздушная частица, перемещающаяся вверх или вниз, принимает температуру окружающего воздуха и ее плотность становится равной плотности окружающих ее частиц, следовательно, нет причин к ее подъему или опусканию. В таком случае говорят, что атмосфера находится в состоянии безразличного равновесия, или что наблюдается равновесная стратификация. В естественных условиях влажного воздуха равновесное состояние наблюдается при меньшем градиенте температуры (примерно 0.6 градусов на сто метров).

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

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

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

Возрастание температуры с высотой называется инверсией температуры. Инверсия температуры может наблюдаться как у поверхности земли - приземная инверсия (нижняя граница совпадает с земной поверхностью) и приподнятая инверсия (нижняя граница расположена на некоторой высоте). Изотермия характеризуется ровным ходом температуры воздуха с высотой. Слои с изотермией и инверсией температуры являются задерживающим слоем для динамической турбулентности и конвекции, поэтому в значительной степени способствуют накоплению вредных примесей в приземном слое атмосферы. В большом городе повторяемость и приземных, и приподнятых инверсий высокая, однако соотношение между видами инверсий противоположно. [2]

1.2 Энергетический источник

В данной курсовой работе Томская ГРЭС-2 является источников загрязнения атмосферного воздуха. На Томской ГРЭС-2 действует десять паровых котлов и пять паровых турбин. Установленная электрическая мощность Томской ГРЭС-2 составляла 331 МВт, тепловая - 815 Гкал/час. Основными видами топлива Томской ГРЭС-2 являются природный газ и каменный уголь. На Томской ГРЭС-2 сжигается уголь Кузнецкого угольного бассейна (например, Моховского и Талдинского угольных разрезов), преимущественно марок «Д», «ДГ». Начиная с 2006 г. на Томской ГРЭС-2 наблюдается устойчивая тенденция увеличения доли используемого угля в общем топливном балансе. В 2012 г. доля угля составляла 48% (341,2 тыс. тонн), против 17,6% (200 тыс. тонн) угля в 2006 г. в корзине энергоресурсов теплоэлектростанции. Средняя зольность угля, сжигаемого на Томской ГРЭС-2, составляет 17,0%, содержание серы - 0,35%, азота - 1,75%, массовая доля минеральных примесей - 2-4%, калорийность угля - 5160 ккал/кг.

Основным видом воздействия Томской ГРЭС-2 на состояние воздушного бассейна города является выброс загрязняющих веществ, образующихся при сжигании в котлах органического топлива. При сжигании твердого топлива с дымовыми газами в атмосферу выделяются: летучая зола углей, диоксид и оксид азота, диоксид серы, бенз(а)пирен, соединения Fe, Mn, Ni, Cr, сероводород, фтористые соединения, углеводороды. Пылегазоулавливающие установки Томской ГРЭС-2 предназначены для улавливания золы и частичного улавливания диоксида серы в мокрых золоуловителях, средняя эффективность очистки составляет 96%. В итоге в атмосферный воздух г. Томска ежегодно выбрасывается до 7,05 тыс. тонн загрязняющих веществ, из которых 32% (2,23 тыс. тонн) - это твердые вещества. Отвод дымовых газов на Томской ГРЭС-2 осуществляется через дымовую трубу высотой 100 м с диаметром устья 6 м. [3]

1.3 Преобразование вредных выбросов ТЭЦ в атмосферном воздухе

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

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

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

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

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

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

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

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

NO₂+УФ=NO+O;

O+O₂=O₃;

NO+O₂=NO₃;

NO₃+O₂=NO₂+O₃.

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

В итоге происходит непрерывное образование озона, который взаимодействуя с оксидом азота, образует снова диоксид азота, т.е.

NO+O₃=NO₂+O₂.

Как показывают исследования, в результате перечисленных реакций происходит постепенное доокисление монооксида NO до диоксида NO₂ по мере удаления дымового факела от дымовой трубы. На выходе из дымовой трубы 85…90% всех оксидов азота представляет NO.

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

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

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

При высоких концентрациях оксидов азота, они частично окисляются под воздействием солнечной радиации до высшего оксида азота N₂O₅, который, взаимодействуя с влагой воздуха, образует азотную кислоту.

Соединения ванадия, аэрозоли бенз(а)пирена, распространяясь в атмосфере вместе с пылью, дождем или снегом, оседают на почву и водоемы.

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

Следует помнить, что почти все выбрасываемые ТЭС вещества не являются инородными для окружающей природной среды и участвуют в круговороте веществ между атмосферой, литосферой и гидросферой.[4]

1.4 Математические модели распространения примеси

1.4.1 Гауссова модель

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

выброс программный моделирование атмосферный

Рис.1 - Идеализированный выброс гауссова типа

Здесь Ф - сечение выброса перпендикулярно направлению движения, h - высота трубы, H - общая высота выброса, которая является суммой высоты трубы и добавочной высоты, на которую поднимается выброс.

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

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

(1)

Где - координаты источника примеси;

Q - мощность источника;

- средние квадратичные отклонения частиц примеси в момент времени t соответственно вдоль координатных осей OX, OY, OZ:

, , (2)

V_x - коэффициент, характеризующий скорость ветра;

h - высота приземного слоя.

Используя принцип суперпозиции, из (1) легко получить формулы для расчета концентрации примеси от точечного источника непрерывного действия, мгновенного линейного источника, линейного источника непрерывного действия, мгновенного площадного источника, площадного источника непрерывного действия, мгновенного объемного источника, объемного источника непрерывного действия, если источник является точечным непрерывного действия, то:

, (3)

где

- координаты источника примеси;

Q - мощность источника примеси;

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

Если источник представляет собой более сложное чем точка тело, то и и - абсциссы начала и конца источника соответственно;

и - ординаты начала и конца источника соответственно;

и - аппликаты начала и конца источника соответственно.

При рассмотрении источников примеси непрерывного действия вводится параметр - время действия источника (предполагается, что источник примеси начал свою работу в момент времени (t = 0).

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

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

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

1.4.2 ANSYS/FLOTRAN

Модуль Flotran конечноэлементного комплекса ANSYS позволяет решать различные гидро- и аэродинамические задачи, включая задачи о распространении примесей, в том числе и задачи об исследовании процессов распространения выбросов в приземном слое атмосферы. Имеется возможность учитывать тепломассоперенос и турбулентность. В этом случае модуль ANSYS/Flotran использует уравнения непрерывности для каждой примеси, уравнение Навье-Стокса и уравнение энергии.

Модуль позволяет рассматривать задачи: с одним источником выброса и одной примесью, с одним источником выброса и двумя примесями, с двумя источниками выброса, решенными при разных граничных условиях и получать графики распределения концентраций примесей в различные моменты времени после выброса. [6]



2. ПРАКТИЧЕСКАЯ ЧАСТЬ

2.1 Постановка физической задачи

Отвод дымовых газов осуществляется через дымовую трубу высотой 100м с диаметром устья 6м. На ГРЭС-2 сжигается уголь Кузнецкого угольного бассейна. Средняя зольность угля, сжигаемого на Томской ГРЭС-2, составляет 17,0%, содержание серы - 0,35%, азота - 1,75%, массовая доля минеральных примесей - 2-4%, калорийность угля - 5160 ккал/кг.

Рассматриваемая область 150x300м. Время расчета 10 мин. Загрязняющее вещество диоксид серы - соединение серы с кислородом состава SO₂. В нормальных условиях представляет собой бесцветный газ с характерным резким запахом. ПДК содержания серы в воздухе рабочей зоны - 10 мг/м³.Необходимо исследовать зависимость изменения концентрации диоксида серы в зависимости от скорости ветра, значения которой задаем в диапазоне от 1 до 5 м/с.

Н - высота трубы; D - диаметр устья; ВГ - верхняя граница; ЛГ - левая граница; НГ - нижняя граница; ПГ - правая граница.

Рисунок 2 - Схема расположения источника выброса

2.2 Постановка математической задачи

В качестве математической модели решения задачи используем дифференциальное уравнение (1):

с+ u + u = (с + () + S (1)

где,

с - плотность загрязняющего вещества;

D_x,y - коэффициент диффузии на оси x, y;

u, - профили скорости ветра на оси х, y;

S - источник загрязняющего вещества;

Начальные условия:

С_t=0 = 0 (С_ф)

Граничные условия:

ЛГ: С_х=0 = 0

ПГ: _х=х = 0

НГ: С_y=0 = 0

ВГ: _y=y = 0

Задаем источник:



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

Наилучшим методом для получения консервативной системы уравнений является метод контрольного объема, смысл которого заключается в следующем: геометрическое и временное пространство разбивается на конечное число объемов, затем для каждого из них записывается баланс субстанции (энергии, импульса, массы и т.д.), т.е. разница величины рассматриваемой вами субстанции поступившей и покинувшей данный объем за определенный промежуток времени (интервал). Что самое важное в расчете, это то, что баланс для каждого контрольного объема в общем случае связан с балансом еще шести смежных с ним объемов. В расчетной области требуется установить граничные условия. Полученная после система алгебраических уравнений решается итерационным способом, далее результаты решения обрабатываются и представляются в удобном для анализа виде. Классический метод конечных (контрольных) объемов (FiniteVolumesMethod (FVM); МКО) основан на интегральных законах сохранения. На первом этапе для любого конечного объема формулируется закон сохранения. Затем расчетная область покрывается сеткой, в узлах которой будут рассчитываться физические характеристики (параметры) моделируемого процесса. Далее выбираются контрольные объемы, чаще всего, с центрами в узлах расчетной сетки и границами, проходящими через центры ребер ячеек сетки. Для каждого полученного контрольного объема записывается дискретный аналог закона сохранения на основе баланса всех потоков через границы рассматриваемого объема. Метод конечных объемов в большинстве случаев позволяет получать консервативные схемы, допускает дискретизацию расчетных областей со сложной геометрией, а также позволяет строить более точные схемы вблизи границ области. Эти достоинства метода обусловлены возможностью использовать нерегулярные сетки, равно как и контрольные объемы произвольной формы.

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

Дифференциальные уравнения, описывающие процессы тепло- и массообмена и гидродинамики подчиняются обобщенному закону сохранения. Если обозначить любую искомую функцию Ф, то обобщенное дифференциальное уравнение принимает в тензорной форме вид (формула 2):

(2)

гдеt, - временная координата, x_i- и пространственная координата, с - плотность, v_i - компоненты вектора скорости, Г_Ц- коэффициент переноса (например, Г_Ц - коэффициент турбулентной вязкости, теплопроводности, диффузии и т.д.), S_Ц - источниковый член. В частности, в S_Цможет входить приток (сток) тепла за счет химических реакций в уравнении энергии или увеличение (уменьшение) концентраций компонент в результате химических реакций в уравнении диффузии. Конкретный вид Г_Ци S_Ц зависит от смысла переменной Ф.

Физически метод контрольного объема можно описать следующим образом: расчетная область разбивается на определенное число непересекающихся контрольных объемов таким образом, что каждая узловая точка содержится в одном контрольном объеме. В двумерном случае мы рассматриваем прямоугольник. Далее, дифференциальные уравнения интегрируются по каждому контрольному объему. Для вычисления интегралов используются кусочные профили, которые описывают изменение функции Ф между узловыми точками. Полученный дискретный аналог выражает закон сохранения для параметра состояния Ф для каждого конечного контрольного объема. Важнейшее свойство метода контрольного объема состоит в том, что при его использовании точно выполняются интегральные законы сохранения таких величин как масса, количество движения, энергия и т.д. в каждом контрольном объеме и для любой группы контрольных объемов и, следовательно, на всей расчетной области. Решение с использованием малого количества контрольных объемов также удовлетворяет точным интегральным балансам во всей области, то есть дискретный аналог (разностная схема) Патанкара-Сполдинга является консервативным. Способ разбиения расчетной области для двухмерного случая на контрольные объемы и типичный контрольный объем представлены на рисунке 3 [7].

Рисунок 3 - Типичный контрольный объем для двухмерного случая

На рисунке изображена область двухмерной сетки, где пунктирной линией обозначен контрольный объем для двухмерной области с центром в точке P. Центры соседних контрольных объемов обозначены буквами N, S, W, E (соответственно «север», «юг», «запад» и «восток»). Соответствующие точки на границах контрольного объема обозначены буквами n, s, w, e. Размеры контрольного объема обозначаются Дx и Дy. Расстояния до центров соседних контрольных объемов обозначены дx_e, дy_n.Таким же образом по аналогии вводятся обозначения трехмерного (пространственного) случая, то есть для контрольного объема, изображенного на рисунке.

Дискретный аналог для двухмерного случая запишется в следующем виде (формула 3):

, (3)

где: _Е = D_eA(|P_e|)+[|-F_e, 0|]; _W= D_wA(|P_w|)+[|-F_w, 0|]; _N= D_nA(|P_n|)+[|-F_n, 0|]; _S= D_sA(|P_s|)+[|-F_s, 0|]; b = SДxДy+; _p = _E+_W+_N+_S+-S_pДxДy; = ДxДy / Дt.

Расходы (F) и проводимости (D) определяются следующим образом (формула 4):

; ; ; ;

; ; ; . (4)

а числа Пекле (формула 5):

; ; ; . (5)



Выражение типа [A,B] означает, что выбирается максимальный элемент из двух. Для функции A(|P|) рекомендуется схема со степенным законом (формула 6).

A(|P|)=[|0,(1-0,1|P|)5|]. (6)

Коэффициенты в соседних точках _Е,_W, _N, _S, учитывают влияние конвекции и диффузии для граней контрольного объема, которые в свою очередь зависят от массового расхода F и проводимости D.Член а_р^оФ_р^о характеризует известную величину Ф для контрольного объема (для времени t), отнесенную к шагу по времени. Оставшиеся члены можно интерпретировать аналогичным образом.

Для решения дискретного аналога был использован алгоритм Томаса или TDMA (Tri-diagonal-Matrix Algorithm - трехдиагональный матрицы алгоритмом). Название TDMA является результатом того, что когда матрица коэффициентов этих уравнений записана, все ненулевые коэффициенты группируются вдоль трех диагоналей матрицы.

Для удобства записи алгоритма введем некоторые обозначения. Присвоим узловым точкам номера 1, 2, 3, ..., n. Номера 1 и n относятся к точкам на границе. Тогда дискретный аналог можно записать в следующем виде (формула 7):

aiTi = biTi+1+ ciTi-1 + di, (7)

где i = l, 2, 3, ..., n. Таким образом, температура Т, связана с соседними значениями T_i-1 и Т_i+1.

Запись уравнений для узловых точек на границе дает (формула 8):

c1 = 0; bn = 0. (8)



Следовательно, температуры Т₀ и Т_n+1 не будут, иметь, смысла (в том случае, когда температуры на границе заданы, уравнения для граничных точек записываются в обычной форме, например, если T₁ задано, имеем а₁b₁ = 0, c₁ = 0 и d₁ равно заданному значению Т₁). [7]

2.4 Анализ полученных результатов

Для достижения поставленной цели были решены следующие задачи:

1. Физическая постановка задачи, т.е. создание физической модели процесса.

2. Математическая постановка задачи: определение математической модели процесса, использование уравнения диффузии с соответствующими граничными и начальными условиями.

3. Выбор численного метода решения.

4. Разработка алгоритма решения задач.

5. Составление программы, которая реализует алгоритм решения задачи на языке MATLAB(приложение А).

6. Анализ полученных результатов: проверка адекватности полученных результатов модели, их физическая интерпретация.

Получены результаты изменения концентрации диоксида серы и динамики выброса в атмосферу с течением времени, которые показывают зависимость данных параметров от скорости ветра (рисунки 4-8).

С распространением загрязняющего облака в заданной области действия источника (с учетом того, что рассматриваемый источник находится на высоте 100 м, в программе источник задается в следующем диапазоне через узловые точки: i = 50:62, j = 200:201) при скоростях ветра от 1 до 5 м/с, мы получаем концентрации диоксида серы в пределах значений от 0,001 до 0,1 кг/м³.



h - высота источника, м; х - длина рассматриваемой области; 1, 2, 3 - концентрации оксида азота 0,1, 0,01, 0,001 соответственно.

Рисунок 4 - Изменение концентрации оксида азота (=1 м/с)



h - высота источника, м; х - длина рассматриваемой области; 1, 2, 3 - концентрации оксида азота 0,1, 0,01, 0,001 соответственно.

Рисунок 5 - Изменение концентрации оксида азота (=2 м/с)



h - высота источника, м; х - длина рассматриваемой области; 1, 2, 3 - концентрации оксида азота 0,1, 0,01, 0,001 соответственно.

Рисунок 6 - Изменение концентрации оксида азота (=3 м/с)



h - высота источника, м; х - длина рассматриваемой области; 1, 2, 3 - концентрации оксида азота 0,1, 0,01, 0,001 соответственно.

Рисунок 7 - Изменение концентрации оксида азота (=4 м/с)



h - высота источника, м; х - длина рассматриваемой области; 1, 2, 3 - концентрации оксида азота 0,1, 0,01, 0,001 соответственно.

Рисунок 8 - Изменение концентрации оксида азота (=5 м/с)

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



ЗАКЛЮЧЕНИЕ

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

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

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

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



СПИСОК ИСПОЛЬЗОВАННЫХ ИСТОЧНИКОВ

1. Влияние ТЭС на окружающую среду // saveplanet.su URL: http://www.saveplanet.su/articles_114.html (дата обращения: 16.04.17).

2. Влияние инверсии температуры на содержание примесей в воздухе // http://gossmi.ru URL: http://gossmi.ru/page/gos1_55.htm (дата обращения: 16.04.17).

3. Эколого-геохимическая обстановка в районах расположения объектов теплоэнергетики по данным изучения нерастворимой и растворимой фаз снега (на примере Томской области) // http://earchive.tpu.ru URL: http://earchive.tpu.ru/bitstream/11683/30541/1/dis00025.pdf (дата обращения: 16.04.17).

4. Преобразование вредных выбросов ТЭС в атмосферном воздухе // refdb.ru URL: https://refdb.ru/look/1448551-p3.html (дата обращения: 16.04.17).

5. Гауссова модель распространения примеси для различных источников загрязнения // ecoloclub.ru URL: http://www.ecoloclub.ru/ecenops-455-1.html (дата обращения: 16.04.17).

6. Особенности моделирования процессов распространения выбросов в приземном слое атмосферы с использованием ansys/flotran // school2003.icmm.ru URL: http://school2003.icmm.ru/book/274.htm (дата обращения: 16.04.17).

7. Патанкар С.В. Численные методы решения задач теплообмена и динамики жидкости. М.:Энергоатомиздат, 1984. - 152 с.



ПРИЛОЖЕНИЕ А

Листинг программы

function y=Maxim(n,m)

%исходные данные, граничные условия, нач. условия

P=1.;

Q=1.85;

%поменять исх. данные

xd=300;

zd=150;

tk=600;

M=500;

V=45000;

tau=0;

dt=0.1;

ue=1;%горизонт профиль скорости для каждой скорости ветра (5 штук)

h=10.;

% нулевые контр. Объемы, вычисляем р-ры КО

dx(1)=0;

dx(n)=0;

for i=2:n-1

dx(i)=xd/(n-2);

end

%расстояние между узловыми точками

for i=1:n-1

hx(i)=(dx(i)+dx(i+1))*0.5;

end

dz(1)=0;

dz(m)=0;

for j=2:m-1

dz(j)=zd/(m-2);

end

for j=1:m-1

hz(j)=(dz(j)+dz(j+1))*0.5;

end

%Задаем координаты для построения графика

x(1)=0;

for i=2:n

x(i)=x(i-1)+hx(i-1);

end

z(1)=0;

for j=2:m

z(j)=z(j-1)+hz(j-1);

end

for j=1:m

for i=1:n

%задаем исходные данные

gam(i,j)=0.4*z(j)+0.00001;

%профиля скоростей ветра

u(i,j)= ue*(z(j)/h);

v(i,j)=0;

r(i,j)=1.2;

%источник разбивается на 2 части (линейную(пост) и остальное (некий коэф.))

sc(i,j)=0;

sp(i,j)=0;

end

end

%основная часть для расчета, рассчитываем источник

for i=50:62

for j=200:201

sc(i,j)=M/V*tk;

end

end

% задаемнач. условия

for i=1:n

for j=1:m

t0(i,j)=0;

t(i,j)=t0(i,j);

tt(i,j)=t(i,j);

end

end

for i=1:n

t(i,1)=0.;

end

% Задаем цикл с условием для времени

while tau<tk

for j=1:m

b(n,j)=0;

c(1,j)=0;

for i=1:n-1

fu=r(i,j)*u(i,j);

du=gam(i,j);

pe=fu/du;

%Полиномиальная функция

pol=(1-0.1*pe)^5;

ap=max(0,pol);

b(i,j)=(du*ap/hx(i)+max(-fu,0))*dz(j);

c(i+1,j)=(du*ap/hx(i)+max(fu,0))*dz(j);

end

end

for i=1:n

e(i,1)=0;

d(i,m)=0;

for j=1:m-1

fv=r(i,j)*v(i,j);

dv=gam(i,j);

pe=fv/dv;

pol=(1-0.1*pe)^5;

ap=max(0,pol);

d(i,j)=(dv*ap/hz(j)+max(-fv,0))*dx(i);

e(i,j+1)=(dv*ap/hz(j)+max(fv,0))*dx(i);

end

end

for i=1:n

for j=1:m

a0=r(i,j)*dx(i)*dz(j)/dt;

a(i,j)=b(i,j)+c(i,j)+d(i,j)+e(i,j)+a0-sp(i,j)*dx(i)*dz(j);

h(i,j)=a0*t0(i,j)+sc(i,j)*dx(i)*dz(j);

end

end

for j=1:m

h(2,j)=h(2,j)+t(1,j)*c(2,j);

c(2,j)=0.;

h(n-1,j)=h(n-1,j)+t(n,j)*b(n-1,j);

b(n-1,j)=0.;

end

for i=1:n

h(i,2)=h(i,2)+t(i,1)*e(i,2);

e(i,2)=0.;

h(i,m-1)=h(i,m-1)+t(i,m)*d(i,m-1);

d(i,m-1)=0.;

end

for i=2:n-1

for j=1:m

a1(j)=a(i,j)*P-b(i,j)*(Q-1.);

b1(j)=d(i,j);

c1(j)=e(i,j);

h11=h(i,j)+b(i,j)*(t(i+1,j)-(Q-1)*t(i,j));

h1(j)=h11+c(i,j)*t(i-1,j)+a(i,j)*t(i,j)*(P-1.);

end

tt=tdma(a1,b1,c1,h1,m-1,2);

for j=2:m-1

t(i,j)=tt(j);

end

end

for j=2:m-1

for i=1:n

a1(i)=a(i,j)*P-d(i,j)*(Q-1);

b1(i)=b(i,j);

c1(i)=c(i,j);

h11=h(i,j)+d(i,j)*(t(i,j+1)-(Q-1)*t(i,j));

h1(i)=h11+e(i,j)*t(i,j-1)+a(i,j)*t(i,j)*(P-1);

end

tt=tdma(a1,b1,c1,h1,n-1,2);

for i=2:n-1

t(i,j)=tt(i);

end

end

tau=tau+dt;

end

lev(1)=0.1;

lev(2)=0.01;

lev(3)=0.001;

tp=t';

contour(x,z,tp,lev,'r');

hold on

up=u';

vp=v';

quiver(x,z,up,vp);

end

%Подпрограмма для нахождения прогоночных коэф. p и q

function [ tf ] = tdma(af,bf,cf,df,nf,n0)

p(n0)=bf(n0)/af(n0);

q(n0)=df(n0)/af(n0);

for i=n0+1:nf

p(i)=bf(i)/(af(i)-cf(i)*p(i-1));

q(i)=(df(i)+cf(i)*q(i-1))/(af(i)-cf(i)*p(i-1));

end

tf(nf)=q(nf);

j=nf-1;

%задаем условия для решения обратного цикла, находим искомую функцию

while j>=1

tf(j)=p(j)*tf(j+1)+q(j);

j=j-1;

end

end

Размещено на Allbest.ru

...