Численное решение некоторых обратных задач с различными типами источников атмосферного загрязнения Текст научной статьи по специальности «Математика»

ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА

2008 Математика и механика № 2(3)

УДК 519.6:551.551.6

Е.А. Панасенко, А.В. Старченко ЧИСЛЕННОЕ РЕШЕНИЕ НЕКОТОРЫХ ОБРАТНЫХ ЗАДАЧ С РАЗЛИЧНЫМИ ТИПАМИ ИСТОЧНИКОВ АТМОСФЕРНОГО ЗАГРЯЗНЕНИЯ1

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

Ключевые слова: математическое моделирование, обратные задачи, параллельная реализация, численные методы.

Физическая постановка обратных задач

Проблема охраны окружающей среды в настоящее время становится одной из важных задач науки. Это связано с интенсивным развитием промышленности, которое приводит к увеличению доли промышленных выбросов, загрязняющих окружающую среду. Одним из способов регулирования уровня загрязнения атмосферного воздуха является контроль интенсивности выброса вредных веществ. Однако даже разветвленная сеть наземных постов наблюдений не всегда в полной мере способна предоставить необходимую для природоохранных служб информацию. Большую помощь здесь может оказать применение методов математического моделирования [1, 2].

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

Теория обратных задач представляет собой активно развивающееся направление современной математики. Несмотря на то, что широкое исследование обратных задач началось сравнительно недавно, в этой области получено большое число существенных результатов. В теории обратных задач сформировался ряд направлений, обусловленных как различными сферами ее приложений, так и типами математических постановок обратных задач [4]. В настоящее время разработано большое количество постановок и методов решения обратных задач, в число ко-

1 Работа выполнена при финансовой поддержке РФФИ, грант № 07-05-01126 и Программы Союзного государства СКИФ-ГРИД.

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

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

Рис. 1. Размещение постов наблюдения (•) в 30-километровой зоне Сибирского химического комбината [3]

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

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

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

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

Пусть моделирование процессов переноса и диффузии примеси выполняется с использованием трехмерного нестационарного уравнения «адвекции - диффузии» [1, 2]:

дС дис дУС дШС д (гдСЛ д ( гдСЛ д( дСЛ _ п ...

— +-------+-----+------= —I Г— 1 + — I Г— 1 +—I К — 1 + Бг, I > 0, (1)

д( дх ду дх дх V. дх ) ду V ду) дх V 2 дх )

где и, Ж, V - компоненты вектора скорости; С - концентрация примеси; Г, К2 - коэффициенты турбулентной диффузии; Бс - источниковый член,

^ = Оо/(г)5(х - х0 )5(у - Уо )5(г - г0), 5() - дельта-функция,

[ 1, для источника с постоянной мощностью,

[ 8^ - г0), для мгновенного источника,

20 - количество выброшенного вещества (мощность источника), г0, х0, у0, г0 -время (только для мгновенного источника) и координаты выброса соответственно.

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

Пусть начальные и граничные условия для уравнения (1) имеют следующий вид:

г = 0: С (г, х, у, х) = С0 (х, у, г);

0 дС 0 т дС 0

х = 0: — = 0; х = Т : — = 0;

дх х дх

0 дС 0 т дС 0 (2)

у =0:“0 у=ту: “0;

х = 0:дС = 0; X = Т : — = 0.

дг х дг

Математическая простановка сопряженной с (1) - (2) задачи получается следующим образом [1]. Уравнение (1) умножается на некоторую функцию *

С (г, х, у, г) и интегрируется по времени (0 < г < Т ) и пространству:

Т Ьу дС пх Ьу дС

п и С* — сЪсУсЫск + п II С* и—сЪсУсЫск +

0000 д 0000 дх

Т Ьу Lz дС ТЬХ Ьу Lz дС

+п II С* V сЪсУсХск + II | | С *Ш---------dzdydxdt =

0000 ду 0000 д

ТЬх Ьу Lz д дС ТЬх Ьу Lz д дС

= I I || С — (Г—)dzdydxdt + | | | | С — (Г— )dzdydxdt +

0 0 0 0 дх дх 0 0 0 0 ду ду

тьх Ьу Lz * д дС тьх Ьу Lz *

+ || | 1С — (К2-----)dzdydxdt + | | | ^С 5Сdzdydxdt.

0000 д дz 0000

Используя формулы интегрирования по частям, приходим к уравнению

ТЬ 1у Ь

П I I С( -

дС д(иС ) д(УС ) д(ЖС )

дї

дх

ду

дг

д (ГдС ) д (ГдС ) д (К дС ))Ж Ж

- — (Г ^) - Т" (Г ^) - Т" (Кг-Т-^Лу^Ж _ дх дх ду ду дг дг

1Х 1У ^ ,|Т ТЬУ Ь; ТЬХ I; ^У

- І І І С■ С | йгйуйх - | | | и-С ■ С| х ЖгЖуЖї - | | | V -С ■ С| ЖгЖхЖї -000 0 000 0 000 о

ТЬх Ьу

Т ЬУ Ьг

Г'^Ґ'1 Ь% Т ЬУ Ь; ^Ґ'1

-|| | Ж ■С* ■ с[2 ЖуЖхЖї + | | | ГС* — йгйуйі - | } | ГС дС

г г г 0 г г г дх 0 г г г

дх

(3)

сЬЖуск -

Т Ьх ЬГ дС

1-ГГ Г ГС* —

0 0 0 ду

ТЬх Ь; дС*

ЖгЖхЛ - Г [ Г ГС-----------

0 0 0 ду

Т Ьх Ьу

Ьх “У . дС

сЬсЬсА + / І | К2С —

0 0 0 дг

сіусіхсії -

ТЬх ЬУ

- П I К 2 С

дС*

0 0 0

дг

тЬ Ьу Ь *

ЖуЖхЖї + ШІ С Бс ЖгЖуЖхЖї.

0 0 0 0

Пусть С* удовлетворяет уравнениям

дс* дис* дvc* джс*-дГдС*-дГдС*-дК дС* _р

дї

дх

ду

(4)

(5)

дг дх дх ду ду дг 2 дг

і _ 1,...,N, ї < Т, с начальным и граничными условиями:

С *(Т, х, у, г) = 0;

х = 0 :иС* + Г— = 0; х = Ь :иС* + Г — = 0;

дх х дх

- -* дС * дС

у = 0:УС +Г— = 0; у = Ь :УС +Г— = 0;

ду У ду

- -0 дС 0 Ь дС 0

г = 0:-----= 0; г = Ь :------= 0,

дг г дг

где Рі =8 (і - ^) 8 (х - хі) 8 (у - уі) 8 (г - гі), і = 1,..., Ы, N - количество точек наблюдения с координатами (хг-, уг-, гг), а їг- - это момент времени, когда проводи-

лись измерения концентрации Сі = С(гг-, хі, у, гі). Предполагаем также, что выполняется условие С* (0, х, у, 2) = 0 .

С учетом (1), (2) и (4), (5) из (3) можно выписать следующие тождества [1]:

т т

Iйг |рСйО = |йг ІЄІБсёО, і = 1,..,N; О = [0,£х]х[0,1у ]х[0,4],

0 е 0 в

которые с учетом вида функции р можно представить как

0 0 0 0

0

т

С = а {С* (г, Х0, у0, г0)/(г)Ж . (6)

0

Здесь С* (/, х, у, 2) - решение сопряженной задачи (4), (5) при

р = 8(*- *г- Жх - X ЖУ - У Ж2 - ^).

Для случая мгновенного точечного источника (6) принимает вид

с = QoCi (0 , Х0 , у0 , 20 ) ,

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

т

С; = бо {С*({, х0, у0, г0 )Ж . о

Если бы нам были известны аналитические решения задач (4), (5), то при N > 3 мы имели достаточно уравнений для определения неизвестных параметров источников выбросов. Однако решения задач (4), (5), как правило, в важных практических случаях удается получить лишь приближенно с использованием компьютерной техники. Кроме того, в предоставленных значениях С может содержаться погрешность измерений. Поэтому алгоритм решения обратных задач для определения параметров источника по известным данным наблюдений и параметрам приземного слоя воздуха сформулируем следующим образом:

1. Численно решаются N сопряженных задач (4), (5) с различными источниковыми членами р , г = 1,...,N .

2. В процессе решения сопряженных задач при Т > г > 0 ищется минимум следующего функционала:

ф(О, г, х, у, г) = Сг - О |С* (г, х, у, г)/(г

г=1 V 0

при ограничениях 2 > 0, 0 < х < Ьх, 0 < у < Ьу , 0 < г < Ь2 .

Найденный минимум функционала даст приближенное решение обратной задачи, т.е.

Ф(бо> *о > хо ’ Уо > 2о) = шт Ф(б> I, X У,2).

0< х< Ьу

о< у<Ьу

0< г < Ь

о < <т

а >о

Численное решение сопряженной задачи

Для численной реализации задачи (4), (5) использовались метод конечного объёма и явные разностные схемы: противопотоковая схема [6], схема МЬи Ван Лира [7] и схема Ботта [8] для аппроксимации адвективных членов. Все рассмотренные разностные схемы тестировались на двумерном нестационарном уравнении «адвекции - диффузии», описывающем перенос примеси в заданном потоке.

При этом было установлено [9], что наиболее точные результаты дают схема МЬи и схема Ботта. Результаты, полученные на основе этих схем, практически совпадают с точным решением, найденным из (1), (2) в предположении, что и, V, Г - константы, Ж = 0, К2 = 0 (рис. 2).

£

Рис. 2. Сравнение точного решения с приближенным, полученным численно с использованием различных разностных схем (постоянный источник, распределение концентрации в сечении поперек основному потоку, g = л]х2 + у2 )

Распараллеливание метода решения обратной задачи

Следует обратить внимание на то, что при решении обратной задачи на каждом шаге по времени решается не одна, а N независимых сопряженных задач (4), (5), следовательно, они могут решаться одновременно, но на различных вычислительных узлах многопроцессорной вычислительной системы. Однако для вычисления значения функционала в процессе поиска его минимума потребуется собирать рассчитанные значения функции C* с i-го узла.

Такие условия проведения численного моделирования заставляют привлекать высокопроизводительную вычислительную технику, в частности вычислительный кластер ТГУ СКИФ Cyberia, на котором установлена библиотека передачи сообщений MPI (Message Passing Interface). Эта библиотека функций дает возможность осуществлять обмен информацией между процессами параллельной программы, запущенными на различных узлах многопроцессорной вычислительной системы.

Распараллеливание метода численного решения обратной задачи производилось с использованием принципа «master-slave» (рис. 3). При такой реализации управляющий процесс (master) передает каждому slave-процессу значения (U,V,W, Г, Kz, x, y, zt, tt), необходимые для расчетов, а подчиненные процессы (slave), получив исходные данные, в свою очередь, ведут расчеты независимо друг от друга и найденные на каждом шаге по времени решения своей сопряженной задачи возвращают управляющему процессу, который проводит вычисления значений ф(Q, t, x, y, z) и ищет минимум этого функционала в (0,Т) х [0, Ьх] х [0, Ьу ] х [0, Lz ].

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

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

В [10] был рассмотрен способ организации параллельных вычислений, при котором управляющий процесс дополнительно к своей вычислительной работе решает одновременно с другими узлами выделенную ему сопряженную задачу. На рис. 4 представлены характеристики ускорения рассмотренных выше параллельных алгоритмов численного решения обратной задачи, полученного при запуске параллельных программ на вычислительном кластере ТГУ СКИФ СуЪепа. Из рисунка видно, что при выполнении параллельной программы на пяти вычислительных узлах при указанных выше параметрах и при условии, что мастер-процесс только подсчитывает суммы для ф(2, г, х, у, z) и ищет его минимум, удается получить ускорение вычислений почти в 4 раза (на рис. 4 - это кривая с закрытыми значками). В случае же, когда шаБ1ег-процесс не только рассчитывает значения функционала ф(2,г, х, у, z) с целью поиска его минимума, но и решает выделенную сопряженную задачу одновременно с другими процессами, получено ускорение в 3,5 раза (кривая с открытыми значками). Разница в ускорении связана с тем, что в последнем случае неравномерность загрузки процессоров более значительна, поскольку при увеличении числа используемых процессоров возрастает количество простаивающих Б1ауе-процессоров.

Число процессов

Рис. 4. Ускорение для двух способов параллельной реализации и отношение времен работы параллельных программ

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

Заключение

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

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

ЛИТЕРАТУРА

1. Берлянд М.Я. Прогноз и регулирование загрязнения атмосферы. Л.: Гидрометеоиздат, 1985. 270 с.

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

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

4. Денисов А.М. Введение в теорию обратных задач. М.: Изд-во МГУ, 1994. 206 с.

5. http://gis.green.tsu.ru/Website/askro/Yiewer.htm

6. Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. М.: Наука, 1984. 152 с.

7. Bott A. A Positive definite advection scheme obtained by Nonlinear Renormalization of the Advective Fluxes // Monthly Weather Review. 1988. V. 117. P. 1006 - 1115.

8. Noll B. Evaluation of a bounded high-resolution scheme for combustor flow computations // AIAA Journal. 1992. V. 30. № 1. P. 64 - 68.

9. Панасенко Е.А. Численное исследование переноса примеси в атмосфере // Третья Всероссийская конференция молодых ученых «Фундаментальные проблемы новых технологий в 3-м тысячелетии». Томск: Изд-во ИОА СО РАН, 2006. С. 582 - 586.

10. Панасенко Е.А., Старченко А.В. Численное решение некоторых обратных задач переноса примеси на многопроцессорных вычислительных системах // Четвертая Сибирская школа-семинар по параллельным и высокопроизводительным вычислениям. - Томск: Дельтаплан, 2008. С. 139 - 148.

Статья принята в печать 18.06. 2008 г.