О математическом моделировании экологических процессов в воздушной среде

УДК 519.717:681.326

О МАТЕМАТИЧЕСКОМ МОДЕЛИРОВАНИИ ЭКОЛОГИЧЕСКИХ ПРОЦЕССОВ  В ВОЗДУШНОЙ СРЕДЕ

В.В. ПЕКУНОВ, Ф.Н. ЯСИНСКИЙ

(Ивановский  государственный энергетический  университет)

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

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

Пусть расчетная область содержит источники тепла и загрязняющих веществ, учтем также наличие постоянных воздушных потоков. Введем в области прямоугольные координаты (x₁, x₂, x₃) таким образом, чтобы ось Ox₃ была вертикальной.

Запишем уравнения Навье — Стокса для трех компонент вектора скорости U с использованием эффективной вязкости n_эфф = n_мол + n_турб, где n_мол — молекулярная вязкость, а n_турб — турбулентная вязкость:

;  j = 1, 2, 3, (1)

где F₁ = 0, F₂ = 0, F₃ = bg T; r — плотность воздуха, b — термический коэффициент расширения воздуха, T — “избыточная” температура, g = 9,81 м/c².

Присоединим уравнения  для давления P и температуры T:

, (2)

, (3)

где D_P и D_T — коэффициенты диффузии давления и температуры; c² = a²r, где a — скорость распространения малых возмущений; a_P и a_T — вспомогательные коэффициенты.

Применим модель турбулентности Абрамовича — Секундова, учитывающую такие важные факторы как предыстория потока, конвективный и диффузионный перенос турбулентных пульсаций:

, (4)

, , ,

где k = 2,0, g = 50,0, b = 0,06, L_min — кратчайшее расстояние до твердой стенки.

Запишем уравнения диффузии для N веществ:

, (5)

где  = U₁,  = U₂,  = U₃ + W_j; W_j  — скорость витания j - вещества; — коэффициент диффузии j - вещества; — вспомогательный коэффициент.

Добавим кинетические уравнения:

, (6)

где q — число реакций, R_k - множество номеров веществ, входящих в правую часть k - реакции, L_k — множество номеров веществ, входящих в левую часть k - реакции, A_k = A_k(T) — константа скорости k - реакции, вычисляемая с помощью уравнения Аррениуса [2].

К уравнениям (1)-(6) присоединяются граничные условия 1 и 2 рода. При  необходимости используются также мягкие и циклические граничные условия.

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

, (7)

применим метод расщепления  по физическим параметрам [1], в данном случае заключающийся в последовательном интегрировании трех уравнений вида:

, (8)

где x — одна из осей x₁, x₂, x₃; U_x — проекция вектора скорости на ось x.

К кинетическим уравнениям (6) применим жестко устойчивый метод  Гира [2], что обусловлено жесткостью системы кинетических уравнений. Алгоритм интегрирования на каждой итерации будет выглядеть следующим образом:

1) вычисление коэффициентов  K_H для уравнений (1)-(5);

2) интегрирование уравнений (1)-(5), считая, что ;

3) интегрирование уравнений (6), в результате чего происходит коррекция значений концентраций веществ в соответствии с происходящими химическими реакциями.

Применим для интегрирования уравнений вида (8) метод скалярной  прогонки. Пусть t — шаг интегрирования по времени, t^s+1 = t^s + t. Введем в расчетной области неравномерную сетку узлов, что позволит детально исследовать отдельные участки области. Пусть h_i — размер ячейки сетки по оси Ox между узлами (i, j, k) и (i + 1, j, k). Для повышения вычислительной устойчивости используем при аппроксимации конвективных членов противоточные производные.

Опустив подробности  вывода, запишем прогоночные формулы:

,

, ,

, , ;

,

, .

Для аппроксимации первых производных при вычислении коэффициентов K_H будем использовать следующую схему:

.

Перейдем к вопросу  о распараллеливании вычислений. Заметим, что для задач небольшой размерности (до 10000 узлов расчетной сетки) распараллеливание вычислений не приведет к заметному росту производительности, так как в этом случае основная часть времени будет уходить не на вычисления, а на пересылки данных. В задачах большой размерности распараллеливание оправданно, о чем свидетельствуют результаты экспериментов, приведенные в работе [3].

Будем использовать термин “процессор”, под которым может подразумеваться либо компьютер в сети, либо отдельный процессор многопроцессорной системы. Объединим n процессоров в топологию “труба”. Последовательно пронумеруем процессоры в “трубе” от 0 до n - 1 так, чтобы первый процессор “трубы” имел номер 0, а последний — номер n - 1. Разделим расчетную область по оси Ox₃ на n равных подобластей, последовательно пронумеруем их от 0 до n - 1 так, чтобы нижняя подобласть имела номер 0, а верхняя — номер n - 1. Пусть каждый i - процессор обрабатывает i - подобласть. Пусть подобласти перекрываются на два слоя узлов сетки. Это позволит каждому процессору при прогонке по оси Ox₃ обрабатывать “свою” подобласть независимо от других процессоров, в противном случае каждому i - процессору пришлось бы ждать результатов работы процессоров с номерами от 0 до i - 1, что свело бы к нулю все усилия по увеличению быстродействия. При наличии же перекрытия, граничные ряды узлов подобласти какого - либо процессора либо являются внутренними рядами для соседних подобластей (обрабатываемых процессорами с номерами i - 1 и i + 1), либо являются граничными рядами для всей расчетной области (для процессоров с номерами 0 и n - 1). Поэтому можно поступить следующим образом:

1.  Перед выполнением прогонки по оси Ox₃ процессоры рассчитают значения интегрируемой функции H^s+1 в граничных рядах подобласти по схеме Головичева. При этом им понадобится информация от “соседних” процессоров (из соседних подобластей), то есть необходима организация обмена данными.
2.  Осуществить прогонку по внутренним рядам узлов подобласти, считая, что на верхней (для процессоров с номерами от 0 до n - 2) и нижней (для процессоров с номерами от 2 до n - 1) границах подобласти действуют граничные условия 1 рода.
3.  В качестве окончательных результатов на данном шаге взять значения, вычисленные во внутренних рядах подобластей, так как значения, вычисленные по схеме Головичева, имеют более высокую погрешность.

Приведем формулу схемы  Головичева:

.

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

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

 

ВЫВОДЫ

1.  Сформулирована математическая модель экологических процессов в воздушной среде и предложена методика численного интегрирования.
2.  Предложен алгоритм распараллеливания вычислений для повышения эффективности моделирования.

 

ЛИТЕРАТУРА

1.  Годунов С.К., Рябенький В.С. Разностные схемы.- М.: “Наука”, 1973.- 400 с.
2.  Полак Л.С., Гольденберг М.Я., Левицкий А.А. Вычислительные методы в химической кинетике. - М.: “Наука”, 1984.- 280 с.
3. Ясинский Ф.Н., Чернышева Л.П., Пекунов В.В. Математическое моделирование с помощью компьютерных сетей: Учебное пособие.- Иваново: Изд-во ИГЭУ, 2000.- 201 с.