Параллельное решение задачи численного моделирования распространения загрязнений в воздушном бассейне большого города в окрестности предприятия

( The Parallel Solution for the Task of Numerical Modeling of Pollution Expansion in the Air Space of Cities and Big Plants

Preprint, Inst. Appl. Math., the Russian Academy of Science)

ВВЕДЕНИЕ

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

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

Большинство известных нам отечественных и зарубежных программных продуктов (например, GAS DYNAMICS TOOL, ADAM, CAL3QHC, ISC‑3, PANACHE, REMSAD, UAM‑IV [7, 11]), позволяющих рассчитывать распространение загрязнений, ориентировано на однопроцессорные системы. Аналогичные программные разработки (например, ECOSIM [8] и MAQSIP [4]), предназначенные для многопроцессорных систем, немногочисленны, причем распараллелены лишь отдельные модули. Для обеспечения наивысшей эффективности вычислений необходима более высокая степень распараллеливания, что возможно лишь при организации системы как единого целого. Повышение эффективности вычислений позволит использовать более сложные и точные математические модели и методы решения.

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

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

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

3.  Провести численные эксперименты (в параллельном режиме) по моделированию распространения загрязнений и сравнить полученные результаты с реальными данными.

1. Математическая модель

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

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

                 ,       (1)

где H — одна из основных переменных модели (в их число входят: U_i — проекция вектора скорости воздуха на ось x_i; P — давление; T — температура; n_турб — турбулентная вязкость; Сj — концентрация j-го вещества; U_Pi — проекция вектора скорости пыли на ось x_i; r_p — плотность пыли); x₁, x₂, x₃, t — переменные Эйлера; R_H и a_H — коэффициенты; K_H — функция, вид которой зависит от уравнения:

;  i = 1, 2, 3,

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

; c² = a²r,

где a — скорость распространения малых возмущений;

,

,

, , ,

где n_мол — молекулярная вязкость; g = 50,0; b = 0,06; L_min — кратчайшее расстояние до стенки;

,

где N — число веществ; член  выражает изменение концентрации j‑го вещества в результате химических реакций;

;

;  i=1,2,3;

g₁ = 0, g₂ = 0, g₃ = g;

,

где C_D — коэффициент сопротивления при заданных скоростях потока и частицы; r_pвещ — плотность вещества частицы; a_p — радиус частицы. Коэффициент сопротивления определяется по следующей формуле [5]:

где Re_p — локальное число Рейнольдса для сферической частицы пыли.

Выбор модели турбулентности осуществлен по результатам численных экспериментов с различными моделями (K‑E, Прандтля, Кармана, Риварда, Абрамовича-Секундова) для задачи о турбулентном смешении затопленной плоской незакрученной струи. Результаты экспериментов сравнивались с аналитическим решением Толмиена (в начальном участке) и эмпирическими данными Г.Н. Абрамовича (в переходной зоне и основном участке) [3]. Была выбрана модель Абрамовича-Секундова, показавшая хорошую точность, учитывающая перенос турбулентных пульсаций и предысторию потока, корректно работающая вблизи стенок (без введения дополнительных соотношений).

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

Для уравнений химической кинетики применим: а) метод Гира для жесткой подсистемы кинетических уравнений; б) явно-неявную схему Рожкова, обладающую достаточными точностью и устойчивостью при меньшей вычислительной трудоемкости, в остальных случаях.

2. Распараллеливание вычислений

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

а) если «вычислительная жесткость» [, где Q_j — трудоемкость вычисления K_H для j‑го уравнения вида (1), измеряемая в количестве арифметических операций] высока, то оценивается отношение времен выполнения итерации для методов РПФ и РПП:

                                      ,                            (2)

где N_ур — число уравнений вида (1) [при расчетах по (2) рассматривается N_ур процессоров]. При G₁ < 1 выбирается РПФ, иначе — РПП.

б) иначе анализируется структура орграфа взаимозависимостей переменных подсистемы динамических уравнений вида (1). Если есть подграф, являющийся пустым, или содержащий незначительное число дуг и не связанный с остальной частью графа, то для соответствующих переменных используется РПФ, для остальных — РПП. Если такого подграфа нет, то оценивается отношение максимальных объемов пересылок для РПФ и РПП, приходящихся на долю одного процессора, в системах с индивидуальными каналами связи:

                                        ,                            (3)

где  и — количества обменов данными на итерации для РПП и РПФ соответственно; P_i — степень i-й вершины орграфа взаимозависимостей; N_слоев — число слоев узлов, передаваемых в ходе одной пересылки (РПП); N_z — число слоев узлов по оси, по которой разбивается область (РПП). В системах с разделяемой коммуникационной средой оценивается отношение общих объемов пересылок для РПФ и РПП:

                                  ,                    (4)

где N_П — число процессоров. При G₂ < 1 выбирается РПФ, иначе — РПП.

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

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

Интегрирование уравнений химической кинетики при наличии в расчетной области значительных тепловых неоднородностей требует динамической балансировки загрузки процессоров. Алгоритм централизованный балансировки малоэффективен в связи с высокой загрузкой управляющего процессора. Известные алгоритмы распределенной балансировки [1] обладают невысокой коммуникационной эффективностью в связи с наличием множества «коротких» передач данных. Предлагается новый алгоритм распределенной балансировки, отличающийся минимальным числом пересылок. Алгоритм применим при медленном дрейфе «горячих» участков, когда , где h_min — минимальное расстояние между узлами расчетной сетки; t — шаг интегрирования по времени; V_дрейф — максимальная скорость дрейфа «горячих» участков; k — коэффициент, k > 1. При этом для распределения нагрузки в начале итерации s + 1 с высокой степенью достоверности можно использовать информацию о времени обработки узлов на предыдущей итерации s. Считаем, что время прямо пропорционально числу итераций метода Гира. Алгоритм включает три этапа:

1. Каждый i‑й процессор подсчитывает предполагаемое суммарное время L_i на обработку узлов, выделенных ему при блочном разбиении области, и передает значение L_i остальным процессорам. Вычисляется предполагаемое среднее время выполнения итерации H.

2. Каждый i‑й процессор определяет свой статус для данной итерации. Если  больше допустимого дисбаланса, то i‑й процессор будет участвовать в обменах данных. Если L_i > H, то он должен передать часть своих узлов другим процессорам. В противном случае процессор должен будет принять часть узлов от других процессоров. Число узлов и номера процессоров — партнеров по обмену рассчитываются по единой схеме всеми процессорами, в первую очередь распределяются наиболее длительно обрабатываемые узлы.

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

В следующей таблице приведены сведения об эффективности распараллеливания Q(N_П) и ускорении S(N_П) для кинетического блока на МВС‑1000/16. Очевидно, что при использовании предложенного алгоритма балансировки эффективность повышается на 8% для двух процессоров, на 10% для четырех, на 4% для восьми. При этом даже на восьми процессорах данный алгоритм показал лучшие результаты (на 12%) по сравнению с другим алгоритмом распределенной балансировки [1], также реализованным в экспериментальных целях.

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

На базе предложенных модели и алгоритмов разработаны однопроцессорная и многопроцессорная версии программы моделирования загрязнений, предназначенные для исполнения в операционных системах Windows, Linux и Parix. Многопроцессорный вариант оптимизирован для работы с многопроцессорными системами МВС-1000, Power XPlorer, а также с компьютерными сетями (Windows- и Linux-кластерами). Используются следующие интерфейсы распараллеливания: MPI, TCP Router, Embedded Parix. Показаны хорошая общая эффективность распараллеливания и высокая степень ускорения на МВС‑1000/16 (см. рис. 1, 2). При использовании TCP Router данные показатели имеют несколько более высокие значения, чем при использовании MPI. Это объясняется большей эффективностью средств передачи данных, ориентированных на конкретный класс систем и не расходующих часть ресурсов на поддержание гибкости и универсальности, характерных для MPI.

Рис. 1. Эффективность распараллеливания при различных N_П

Рис. 2. Ускорение вычислений при различных N_П

3. Результаты экспериментов

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

В первом эксперименте объектом моделирования является перекресток двух улиц города Лондона. На нижней границе заданы распределения температуры и концентраций (пыли и первичных загрязнителей): наивысшие значения в районе перекрестка, включая стоп-линии, меньшие — на дорогах вне перекрестка. Распределение концентраций на нижней границе вычислено с помощью программы CAL3QHC [11], для которой средний объем движения рассчитывался по данным отчета об интенсивности движения в Великобритании [10], а интенсивность выбросов в режиме холостого хода определялась с помощью модели MOBILE5 [6]. Использовался набор из 21 реакции, характерный для фотохимического смога: базовые реакции образования озона и пероксиацетилнитрата; реакции с участием оксидов азота, формальдегида и углеводородов C_nH_2n+2, замыкающие цепи циклических механизмов накопления озона; другие реакции, в том числе с образованием паров азотной кислоты. Пыль — частицы песка радиусом до 0,5 мм.

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

Рис. 3. Распределение скоростей во фронтальном сечении вдоль одной из улиц

Рис. 4. Изолинии избыточной температуры во фронтальном сечении вдоль одной из улиц

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

Наиболее загрязненным участком оказался перекресток (наивысшие концентрации первичных загрязнителей, см. рис. 5) и прилегающая к нему зона (умеренные концентрации первичных и наивысшие концентрации вторичных загрязнителей, см. рис. 6).

Рис. 5. Изолинии концентраций угарного газа в горизонтальном сечении на уровне станции наземных измерений

Рис. 6. Изолинии концентраций озона в горизонтальном сечении на уровне станции наземных измерений

На рис. 5 и 6 крестом отмечено местоположение автоматической станции наземных измерений, концентрации здесь и далее указываются в млн^‑1. Произведено сравнение полученных результатов с данными станции наземных измерений. Показано хорошее совпадение по угарному газу CO (рассчитано 0,8¸1,4 млн^‑1, реально 1,1¸1,2 млн^‑1), особенно если сделать поправку на неточность исходных данных. Результаты по CO более точны по сравнению с результатами, выданными известной программой CAL3QHC (2,9 млн^‑1), использующей простую гауссовскую модель. Хорошим является совпадение и по диоксиду азота NO₂ (рассчитано 0,045¸0,08 млн^‑1, реально 0,05¸0,054 млн^‑1). Менее точны результаты по оксиду азота NO, что может быть связано с недостаточно точными исходными данными и неучетом каких-либо реакций с участием NO. Однако по сумме оксидов азота NO_x совпадение также является хорошим. Получены распределения концентраций пероксиацетилнитрата и паров азотной кислоты.

Во втором эксперименте объектом моделирования является окрестность предприятия British Steel complex вблизи города Scunthorpe. Моделируется распространение диоксида серы SO₂ и мелкой пыли, выбрасываемых из трубы предприятия. На верхней границе задан горизонтальный воздушный поток, направленный в сторону города, его скорость выбрана с расчетом, чтобы на высоте выпускного отверстия трубы скорость ветра составляла 2¸4 м/с. Область является полностью открытой. Фактором химической кинетики пренебрегаем, так как влажность воздуха считается низкой, следовательно, развитие влажного смога невозможно. Значения концентрации диоксида серы и плотности пыли на выходе трубы рассчитаны приближенно, на основании годовых выбросов. Точные данные о параметрах частиц пыли также отсутствуют.

В верхней части области образовался горизонтальный поток, уносящий выбросы на значительное расстояние (см. рис. 7). При этом турбулентная вязкость в верхней части имеет значительную величину, а частицы пыли взяты малыми (радиусом около 10^‑4 м), поэтому сила сопротивления частиц потоку велика и оседания пылевого потока на поверхность фактически не произошло. Иная ситуация с диоксидом серы, где за счет диффузии и постоянной скорости оседания загрязненность атмосферы у поверхности более высока (см. рис. 8).

Полученные результаты сравнивались с данными станции наземных измерений. С учетом отсутствия точных исходных данных о величине суточных выбросов, а также искажений под влиянием граничных условий (см. рис. 7, 8), следует признать приемлемым совпадение по SO₂ (рассчитано 3 млн^‑1, реально не выше 0,12 млн^‑1). Значительно менее реальны результаты по пыли (рассчитано 10^‑12 кг/м³, реально 7×10^‑7¸2×10^‑8 кг/м³), что объясняется отсутствием точных исходных данных о составе пыли, а также большим расстоянием от предприятия до станции наземных измерений (около 500 метров), что делает возможным влияние неучтенных источников пыли, например, автомагистралей.

Рис. 7. Распределение скоростей во фронтальном сечении, проходящем через трубу предприятия

Рис. 8. Изолинии концентраций SO₂ во фронтальном сечении, проходящем через трубу предприятия

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

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

2. Сформулирована четкая методика выбора метода распараллеливания при интегрировании уравнений вида (1) на основе оценки «вычислительной жесткости» и количества пересылок (в зависимости от структуры взаимосвязей между переменными).

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

4. Проведено численное моделирование распространения загрязнений в воздушном бассейне большого города (перекресток улиц Лондона) и в окрестности предприятия British Steel complex. Показано хорошее соответствие полученных результатов реальным данным.

Библиографический список

1.  Корнилина М.А., Якобовский М.В. Динамическая балансировка загрузки процессоров при моделировании задач горения // Высокопроизводительные вычисления и их приложения: Тр. Всеросс. науч. конф. — М.: Изд-во МГУ, 2000. — С.34-39.

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

3.  Турбулентное смешение газовых струй/ Абрамович Г.Н., Крашенинников С.Ю., Секундов А.Н., Смирнова И.П. — М.: Наука, 1974.— 272 с.

4.  Brief Description of MAQSIP Science. Draft Version 2.1. — MCNC North Carolina Supercomputing Center. — http://www.emc.mcnc.org/products/maqsip/users/science.html.

5.  Bruse, M. (2002), ENVI-met implementation of the gas/ particle dispersion and deposition model PDDM. — http://www.geographie.ruhr-uni-bochum.de/agklima/envimet/documents/sources.PDF.

6.  Evaluation of MOBILE vehicle emission model. — FHWA-PD-94-38, Federal highway administration environmental analysis division/office of environment and planning, December, 1994. — 204 p.

7.  Guideline on Air Quality Models (Revised) and Supplements. — EPA-450/2-78-027R et seq., published as Appendix W to 40 CFR Part 51 (7-1-99 Edition). U. S. Environmental Protection Agency, Research Triangle Park, North Carolina, 1999.

8.  Mieth, P., Unger, S., Moussiopoulos, N. et al. ECOSIM Telematics Applications Project: ECOSIM Deliverable D05.02. ECOSIM User Documentation. — http://www.ess.co.at/ECOSIM/Deliverables/ D0502.html.

9.  Piva, R., Orlandi, P. Numerical Solutions for Atmospheric Boundary Layer Flows over Street Canyons. Proceedings of the Fourth Int. Conference «Numerical Methods in Fluids Dynamics», ed. R. D. Richtmeyer, Springer-Verlag, New York, pp. 319-325 (1975).

10.  Transport statistics bulletin: Road traffic statistics: 2000 - Statistics report SB (01) 19. — UK Department for transport, local government and the regions, 2001. — 41 p.

11.  User's Guide for CAL3QHC Version 2: A Modeling Methodology for Predicting Pollutant Concentrations Near Roadway Intersections. — EPA-454/R-92-006, U. S. Environmental Protection Agency, Research Triangle Park, North Caroline, 1992.