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

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

2011

Математика и механика

№ 3(15)

УДК 519.6; 001.891.573

А.А. Барт, Д.А. Беликов, А.В. Старченко

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ДЛЯ ПРОГНОЗА КАЧЕСТВА ВОЗДУХА В ГОРОДЕ С ИСПОЛЬЗОВАНИЕМ СУПЕРКОМПЬЮТЕРОВ

Описывается подход, позволяющий при помощи одномерной метеорологической модели и модели переноса примеси прогнозировать качество атмосферного воздуха над урбанизированной территорией. Особенностью метода является использование данных среднесрочного метеорологического прогноза по глобальной модели Гидрометцентра РФ. Модель переноса примеси реализована на суперЭВМ, что позволяет получать прогноз в короткие сроки.

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

Атмосфера Земли - сложная система, в которой ряд физических и химических процессов протекает одновременно. Измерения состояния окружающей среды -это только срез атмосферных условий в конкретное время и для конкретной местности. Часто такие измерения очень сложно интерпретировать без четкой концептуальной модели атмосферных процессов. Одни лишь измерения не могут быть использованы для установления эффективной административной политики решения проблемы загрязнения воздуха, так как информация об отдельных, протекающих в атмосфере процессах (химических, транспортных и др.) не даёт полного представления о системе в целом. Математические модели, в свою очередь, предоставляют необходимый «каркас» для объединения представлений об отдельных атмосферных процессах и их взаимодействия между собой [1].

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

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

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

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

—- +------------------------ +-- +--

дґ дх ду дг

дСг дПЄ1 д¥Сг ЗШСг

(1)

Здесь Ci, с, - осредненная и пульсационная составляющие концентрации i-й компоненты примеси; ось Ох направлена на восток, Оу - на север, t - время, г - вертикальная координата; Ж, м- осредненная и пульсационная составляющие вертикальной компоненты скорости примеси; и, V, и, V - осредненные и пульсацион-ные составляющие вектора горизонтальной скорости; () - осреднение по Рейнольдсу, Si - источниковый член, представляющий выбросы компонентов примеси в атмосферу; Я, описывает образование и трансформацию вещества за счет химических и фотохимических реакций с участием компонентов примеси; о, - скорость влажного осаждения примеси за счет осадков; п - количество химических компонентов примеси, концентрации которых необходимо определить, 0 < t < Т, -Ьх/2 < х < Ьх /2, -Ь/2 <у < Ьу/2, 0 < х < И.

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

Здесь т - временной масштаб турбулентных пульсаций, g - ускорение свободного падения, С10, С2е, С30, Б1С - эмпирические константы, Е - функция, определяющая влияние поверхности на турбулентную структуру течения, 0,0 - осреднённая и

нии, Т - абсолютная температура, Я - газовая постоянная.

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

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

Я

пульсационная составляющие потенциальной

Р0 = 1,013 *105 Н/м2, ср - удельная теплоемкость воздуха при постоянном давле-

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

Для реализации модели переноса примеси необходимо определить для каждого момента времени профили компонент вектора скорости, температуры и турбулентных характеристик. Такого рода данные можно получить усвоением результатов 48-часового метеорологического прогноза, выполняемого с использованием глобальной модели ПЛАВ. ПЛАВ - это полулагранжева модель атмосферы с высоким разрешением (горизонтальная сетка 0,72° - по широте, 0,9° - по долготе, 28 вертикальных слоев), принятая Г идрометцентром в качестве оперативной модели для среднесрочного прогноза с выдачей результатов через 6 часов. В модель ПЛАВ включены параметризации процессов подсеточного масштаба (коротко- и длинноволновая радиация, глубокая и мелкая конвекция, планетарный пограничный слой, торможение гравитационных волн, тепло- и влагообмен с подстилающей поверхностью с учетом растительности), разработанные Метео-Франс и европейским консорциумом LACE (www.rclace.eu) для оперативной региональной модели ALADIN. Прогноз полей метеорологических элементов осуществляется с помощью численного решения уравнений гидротермодинамики в вертикальной сигма-системе координат на сфере. В блоке решения уравнений динамики атмосферы, как и в подавляющем большинстве глобальных моделей численного прогноза погоды, применяется полулагранжево представление адвекции и двухслойная полунеявная схема интегрирования по времени. Оригинальными особенностями блока решения уравнений динамики атмосферы модели ПЛАВ являются применение конечных разностей четвертого порядка на несмещенной сетке для аппроксимации неадвективных слагаемых уравнений и использование вертикальной компоненты абсолютного вихря и дивергенции в качестве прогностических переменных конечно-разностной модели.

В качестве начальных данных модель ПЛАВ-2008 использует поля оперативного объективного анализа на стандартных изобарических поверхностях с горизонтальным разрешением 1,25 градуса по долготе и широте, а также объективные анализы температуры и относительной влажности на уровне 2 м, температуры и влагосодержания поверхностного и глубинного слоев почвы на модельной сетке, полученные на основе усвоения данных только на наземных станциях наблюдений.

Программный комплекс модели распараллелен на основе сочетания технологий MPI и OpenMP, достигнута масштабируемость кода на 512 процессорах [4].

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

(2)

(3)

дQ д . .

— = — (qw), дt дz '

(5)

- компоненты скорости геострофического ветра,

где и, V - компоненты горизонтальной скорости ветра в атмосферном пограничном слое (Ж = 0), Q - влажность, (им), (ум), (0м), (дм) - осредненные по Рейнольдсу корреляции турбулентных пульсаций вертикальной скорости и компонент горизонтальной скорости, температуры и влажности соответственно.

др др' ду ’ дх/

/ = 2О8шу - параметр Кориолиса, у - географическая широта рассматриваемой точки, О - угловая скорость вращения Земли, р - давление.

Для замыкания системы уравнений (2) - (5) применяется трехпараметрическая модель турбулентности, включающая уравнения переноса для энергии, масштаба турбулентных пульсаций и уравнение для дисперсии турбулентных пульсаций потенциальной температуры (02) [3]:

(• Vg )=Pf І"

дк , . ди , .1V g .

— = -< uw)--------(vw)-----1— (wt

дt ' 1z ' 'да ©X

д

~z

а

,-Jkl

дк

дz

Cuk2

cL_

дt

:C"І-{uw) itw IV+t<,,e>)

-+—І а^'/kl д_-1 + CL

k 1z У e 1z ) L

1 -

д(92) дt

= +-

1z

2\Л

C9Jkl (w)2 19-)

1z

- 2(wt)^ - 2^-.

1z Тл

(б)

(7)

(В)

Здесь к = 0,5((ы2)+(у2)+(м’2)) - кинетическая энергия турбулентности; I - интегральный масштаб турбулентности, се = 0,54; Сц = -0,12; С12 = 0,2; Са = 0,19; Се = 3,0 - числовые коэффициенты, те - временной масштаб температурных турбулентных пульсаций.

На верхней границе для определяемых из уравнений (2) - (8) зависимых переменных использовались простые градиентные условия (равенство нулю производной по г, для всех переменных, для потенциальной температуры - равенство нулю второй производной), на нижней границе применялись соотношения теории подобия Монина - Обухова и эмпирические функции параметров подобия [5]. Начальные условия для уравнений (2) - (8) получались в результате обработки данных прогноза, получаемых по модели ПЛАВ.

З

Расчет химических реакций

Моделирование химических и фотохимических реакций в (1) проводилось на основе полуэмпирической кинетической схемы Харли Generic Reaction Set (GRS) [6]. Фотохимическая модель GRS содержит ns = 8 химических компонент (реагирующая часть смога Rsmog, образующиеся в атмосфере аэрозольные частицы APM, органические радикалы RP, перекись водорода H2O2, оксид азота NO, диоксид азота NO2, озон O3, диоксид серы SO2), которые участвуют в 8 химических реакциях. Апробация этой кинетической схемы проводилась для условий образования вторичных загрязнителей в городе Мельбурне и области Пилбара, расположенной на севере Австралии, а также для условий города Томска [2]. Полученные резуль-

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

Члены, входящие в уравнения (1) и отвечающие за химические и фотохимические реакции для схемы GRS, имеют вид [6]

rr т = 0;

L smog J

rapm ] = FCH2 nk1

CRsmog + FHNO3 k7CRPCNO2 + FH2SO4 k8CRPCSO2;

RRP] = k1 CRsmog k2CRPCNO k5CRPCRP k6CRPCNO2 k7CRPCNO2 k8CRPCS‘

RH2O2] =ak5CRPCRP ;

RNO+NO2] = -k6CRPCNO2 - k7CRPCNO2 ;

R NO2] = k2CRPCNO + k4CNOCO3 - k3CNO2 - k6CRPCNO2 - k7CRPCNO2 ;

R O +NO ] k2CRPCNO k6CRPCNO k7 CRPCN

{03+Ш2] Л'2^ЯР^N0 п-6^ЯР^-Ш2 Л7'-'ДР'-'К02 ’

К[302] = -к8СКРСд02 ;

■Рсн = 0,57, = 4,0, -РнЖз = 2,6 - приближенные коэффициенты, оп-

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

' ( [Кт, ] "

золи, п= 0,І, а = max

0,03, exp

-0,026і-

- численные коэффициен-

[N0+N02]

ты [2, 6]. Для расчета компонента примеси С0, который принимался в данной работе химически инертным (К[С0] = 0;), также используется уравнение переноса (1).

Численный метод и его параллельная реализация

Уравнение (1) с учетом представленной выше схемы замыкания можно записать в следующем виде (Ф = С, Б = Б, К=К, с = с,, , - фиксированное значение от 1 до п)

дФ дПФ д¥Ф дЖФ 5 ( дФ ^ дФ ^ дФ^|

---+------+------+------= —I Г„ — + Гху-------+ Гхг---1 +

д( дх ду дг дх \ дх ду дг)

д ( дФ дФ дФ^ д ( дФ дФ дФ^

+—I Гух-------+ гуу-+ Гуг-----|+ —I Ггх---+ Г-----+ Ггг---1-стФ + Б + К. (9)

ду ^ у дх уу дУ дг ^ дг ^ дх ду дг)

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

Для получения конечно-разностного аналога уравнения (9) используется метод конечных объемов. Схема расчетной ячейки представлена на рис. 1. Для аппроксимации диффузионных потоков используется центрально-разностный оператор. Для аппроксимации адвективных членов применяется направленная противопо-токовая схема МЬи [7].

Рис. 1. Вычислительная ячейка с обозначением центров используемых по шаблону ячеек и граней выбранного конечного объема

Дискретный аналог уравнения переноса (9) с учетом принятых аппроксимаций можно записать в виде

Ф р+1 = Б, ФТ+1 + А. Ф к+1 '

'Ъ^Б

где

[аР + аАхДуАгР — (БР + КР Д АхАуА- р + Б, + БЪ

+ \ аЕ Ф Е + аЖ ФЖ + Ф N + аБ Ф Б + аТ Ф Т + аБ Ф Б — аР Ф Р ~\ + Ъ

аЕ = тах (—Ее; 0Д + Ае; а№ = тах ; 0Д + А,; aN = тах(—^;0Д + Ап; а5 = тах(^;0Д + А,;

ат — тах ( ; 0 Д; ав — ^тах (^Ъ; 0 Д ; ар — аЕ + аж + aN + а б + ат + а б ар ;

Ъ = (Рк + Кр Д ДАуАгр + Вшр ; аР = Ах у р

(10)

А/

к+1

Здесь Фкр 1 « Ф(/к + А/, х1, у ^, гР), к - номер слоя по времени, Дх, Ду, Д^ - размеры

конечного объема с центром в точке Р; центры соседних по горизонтали конечных объемов обозначаются по направлениям сторон света - N, Б, Е, Ж, а по вертикали - Т и Б для нижнего и верхнего соответственно; грани конечного объема - п, ,, е, ,, /, Ъ; ¥, Б представляют адвективные и диффузионные потоки через соответствующие грани конечного объема; при аппроксимации источниковых членов

использовалась

тук+1 т)0к . тЛк +1

Лр = Лр + Лр Фр

«линеаризованная» форма

1к о1к

записи:

8І+1 = ^0к + 8р Ф Р+1

Лр

8р < 0 ; Втр содержит смешанные производные на

временном слое к и наклоны в МЬИ-представлении потоков.

Используемая явно-неявная разностная схема (10) позволяет неявным образом провести расчет диффузионных турбулентных потоков вдоль оси Ог и тем самым избежать существенных ограничений на шаг по времени вследствие малых вертикальных размеров вычислительных ячеек, высота которых из-за сгущения вычислительной сетки к поверхности в первом узле составляет всего 2 м. В то же время горизонтальные размеры вычислительной сетки (500 м) позволяют эффективно использовать явную вычислительную схему. Такой подход обеспечивает высокую скорость вычислений вследствие применения экономичного метода прогонки для решения на каждом шаге по времени (Д/ = 6 с) систем линейных уравнений с трехдиагональной матрицей (9).

Расчетная область - равномерная в горизонтальных направлениях (Ох и Оу) и неравномерная, сгущающаяся к поверхности конечно-разностная сетка размером Nx xNy =100x100x100 узлов.

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

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

• р т

!§£.. ..

Ь.І ..

— : ■■■■■

— —

X

Рис. 2. Схема двумерной декомпозиции параллельной реализации явно-неявного метода на примере четырех процессоров (р=4)

На рисунке присутствуют светлые блоки по краям расчетной области, являющиеся фиктивными границами, темные блоки, расположенные со стороны разделения области, являются блоками обмена данными между процессорами. При такой декомпозиции можно использовать большее количество процессоров суперкомпьютера, а доля временных затрат, связанных с необходимостью на каждом шаге по времени передавать данные между процессорными элементами, меньше, чем при одномерной декомпозиции. Расчеты, выполненные на кластере ТГУ СКИФ СуЪепа[8], показали, что при использовании 64 процессоров ускорение вычислений составляет 54, а эффективность соответственно 0,84. При этом один час моделирования на 64 процессорах выполняется в среднем за 20 секунд, в то время как время расчета этой же задачи на одном процессоре - 18 минут.

Сравнение результатов моделирования с данными измерений

Проверка полученных результатов моделирования осуществлялась путем сравнения расчетов с данными измерений, проведенными на ТОР-станции Института оптики атмосферы СО РАН в Томске. Расчетная область представляет собой часть атмосферы размером 50x50x2 км, выбранная таким образом, что ее квадратное основание содержит участок подстилающей поверхности с крупным населенным

о. ^

а &

г \

и м

К ч ' х те

и “ :

и к а

Я “

8 Декабря 2009 г. I ■ • * , 16 Апреля 2010 г.

1 1 1 1 1 1 1 1 1 1 1 1 0 4 8 12 16 20 24 : т т Т т Т : 111111 0 4 8 1 111111 2 16 20 24 ТМ п 1 1.

ПТ I Ш ТВ | х * 1 1 |

0 4 8 12 16 20 24 0 4 8 1 2 16 Я 1 * 0 24

100 315 80 Я 3 60 Ег2 40 §1 20 о 0~* ’ *

0

100 п Ц 80: 60-■§ 3 40-=Г К &

8 12 Время, ч

8 12 Время, ч

0

0

4

8

12

16

20

24 0

4

8

12

16

20

24

4

8

12

16

20

4

16

20

Рис. 3. Сравнение результатов моделирования с данными ТОР-станции 8 декабря 2009 г. и 16 апреля 2010 г.

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

На рис. 3 сплошной линией представлены результаты численного моделирования, символами отмечены результаты измерений, проводимых на ТОР-станции ИОА, для каждого измерения вдоль оси значений приведена погрешность. На графиках дается сравнение измерений и расчетов силы и направления ветра, а также концентрации таких загрязнителей, как монооксид углерода (CO), диоксид азота (NO2) и озон (O3). Из графиков сравнения видно, что математическая модель хорошо предсказывает изменение ветра и значений рассматриваемых концентраций в течение суток.

Заключение

Представленная в работе модель имеет ряд преимуществ, позволяющих использовать ее в реальных условиях в оперативном режиме. Такая возможность появилась благодаря усвоению данных глобального прогноза метеорологических характеристик, а также за счёт параллельной реализации модели переноса примеси при помощи двумерной декомпозиции расчетной области. В настоящее время модель реализована и является ядром программного комплекса, расположенного на выделенном сервере. Ежесуточно программными компонентами комплекса в автоматическом режиме производится получение начальных данных, их обработка, запуск процесса моделирования на кластере ТГУ СКИФ Cyberia [8] и обработка полученных результатов с последующим их представлением в сети Internet в виде полей прогноза распространения в течение суток примеси, наложенных на карту местности[9].

ЛИТЕРАТУРА

1. Seinfeld J.H., Pandis S.N. Atmospheric chemistry and physics. From air pollution to climate change. N.Y.: John Wiley, Sons, Inc., 1998. 1326 p.

2. Беликов Д.А., Старченко А.В. Исследование образования вторичных загрязнителей (озона) в атмосфере г. Томска // Оптика атмосферы и океана. 2005. Т. 18. № 05-06. С. 435-443.

3. Беликов Д.А., Старченко А.В. Численная модель турбулентного переноса примеси в пограничном слое атмосферы // Оптика атмосферы и океана. 2007. Т. 20. № 8, С. 667-673.

4. Толстых М.А. Полулагранжева модель атмосферы с высоким разрешением для численного прогноза погоды // Метеорология и гидрология. 2001. № 4. С. 5-16.

5. Атмосферная турбулентность и моделирование распространения примесей / под ред. Ф.Т.М. Ньистадта и X. Ван Допа. Л.: Гидрометеоиздат, 1985. 350 с.

6. Hurley P.J. The Air Pollution Model (TAPM) Version 2 / CSIRO Atmospheric Research Technical Paper № 55. 2002.

7. Van Leer B. Towards the ultimate conservative difference scheme. Part IV: A new approach to numerical convection // J. Comput. Phys. 1977. V. 23. P. 276-299.

8. Межрегиональный супервычислительный центр ТГУ. URL: http://skif.tsu.ru

9. Барт А.А., Старченко А.В., Фазлиев А.З. Информационно-вычислительная система краткосрочного прогноза качества воздуха над урбанизированной территорией // Материалы Всероссийской молодежной научной конференции «Современные проблемы математики и механики». Томск: Изд-во Том. ун-та, 2010. С. 21-24.

Статья поступила 18.08.2011 г.

Bart A.A., Belikov D.A., Starchenko A.V. SUPERCOMPUTER-BASED MATHEMATICAL MODEL FOR AIR QUALITY PREDICTION IN THE URBAN AREA. A new approach permitting one to predict atmospheric air quality over urban areas using a one-dimensional meteorological model and a 3D pollutant transport model is described. The feature of the method is assimilation of the data of a medium-term meteorological forecast according to the global model of the Hydrometeorological Center of Russia. The pollutant transport model is implemented on a supercomputer, which permits one to obtain the forecast in a short time.

Keywords: transport and formation of secondary air pollutants, mathematical modeling, parallel computing

BARTAndrey Andreevich (Tomsk State University)

E-Mail: bart@math.tsu.ru

BELIKOV DmitryAnatol’evich (Tomsk State University)

E-Mail: dmitry.belikov@nies.go.jp

STARCHENKO Alexandr Vasil’evich (Tomsk State University)

E-Mail: starch@math.tsu.ru