Supercomputer-based mathematical model for air quality prediction in the urban area.pdf Атмосфера Земли - сложная система, в которой ряд физических и химическихпроцессов протекает одновременно. Измерения состояния окружающей среды -это только срез атмосферных условий в конкретное время и для конкретной мест-ности. Часто такие измерения очень сложно интерпретировать без четкой концеп-туальной модели атмосферных процессов. Одни лишь измерения не могут бытьиспользованы для установления эффективной административной политики реше-ния проблемы загрязнения воздуха, так как информация об отдельных, проте-кающих в атмосфере процессах (химических, транспортных и др.) не даёт полно-го представления о системе в целом. Математические модели, в свою очередь,предоставляют необходимый «каркас» для объединения представлений об от-дельных атмосферных процессах и их взаимодействия между собой [1].Успешное решение задач прогноза уровня загрязнения воздуха основано наиспользовании математических моделей, учитывающих физические особенностираспространения примесей в атмосфере, связи между концентрациями примесей иметеорологическими параметрами: скоростью и направлением ветра, инверсией(повышение температуры воздуха с высотой), осадками, туманами.Целью работы является разработка и апробация математической модели дляпрогноза качества атмосферного воздуха над городской территорией, способнойоперативно и c достаточной точностью предсказывать перенос и образование вто-ричных загрязнителей.Математическая модель переноса примесиДля расчета концентрации компонентов примеси с учетом химических реак-ций применяется эйлерова модель турбулентной диффузии, включающая транс-портные уравнения с описанием адвекции, турбулентной диффузии и химическихреакций [2]:; 1, 2,...,i i i ii i ii i i i sC UC VC WCt x y zc u c v c wC S R i nx y z   + + + =     = − − − − ƒ + + =  . (1)Здесь Ci, ci - осредненная и пульсационная составляющие концентрации i-й ком-поненты примеси; ось Ox направлена на восток, Oy - на север, t - время, z - вер-тикальная координата; W, w- осредненная и пульсационная составляющие верти-кальной компоненты скорости примеси; U, V, u, v - осредненные и пульсацион-ные составляющие вектора горизонтальной скорости; 〈〉 - осреднение по Рей-нольдсу, Si - источниковый член, представляющий выбросы компонентов приме-си в атмосферу; Ri описывает образование и трансформацию вещества за счет хи-мических и фотохимических реакций с участием компонентов примеси; ƒi - ско-рость влажного осаждения примеси за счет осадков; ns - количество химическихкомпонентов примеси, концентрации которых необходимо определить, 0 ≤ t ≤ T, -Lx /2 ≤ x ≤ Lx /2, -Ly/2 ≤ y ≤ Ly /2, 0 ≤ z ≤ h.Корреляции 〈ciu〉, 〈civ〉, 〈ciw〉, представляющие турбулентную диффузию, также как концентрации Ci, являются неизвестными, что делает уравнения (1) не-замкнутыми. Для определения неизвестных корреляций необходимо привлекатьдополнительные замыкающие соотношения. В данной модели переноса примесиони выводятся в предположении равновесного приближения для дифференциаль-ных уравнений переноса турбулентных потоков массы в условиях локальной од-нородности атмосферного пограничного слоя и имеют вид [3]( ) 2211 i i i;i ic u C c w U u C vuC wuCC z x y z ƒƒ− = ƒ⎛⎜⎝−  +  +  + ⎞⎟⎠( ) 2211 i i i;i ic v C c w V uvC v C wvCC z x y z ƒƒ− = ƒ⎛⎜⎝−  +  +  + ⎞⎟⎠( ) 231 11 i i i.i iCc w C g c uwC vwC w CC D F x y z ƒƒ− = +ƒ ⎛⎜⎝− − ƒ ƒ +  +  + ⎞⎟⎠Здесь ƒ - временной масштаб турбулентных пульсаций, g - ускорение свободногопадения, C1ƒ, C2ƒ, C3ƒ, D1C - эмпирические константы, F - функция, определяющаявлияние поверхности на турбулентную структуру течения, ƒ,ƒ - осреднённая ипульсационная составляющие потенциальной температуры: 0 pRP cTPƒ = ⎛⎜ ⎞⎟⎝ ⎠,P0 = 1,013*105 Н/м2, cp - удельная теплоемкость воздуха при постоянном давле-нии, T - абсолютная температура, R - газовая постоянная.Для нахождения напряжений Рейнольдса 〈uiuj〉, турбулентных тепловых пото-ков 〈uiƒ〉 и корреляции 〈ci ƒ 〉 используются алгебраические выражения и диффе-ренциальное уравнение [3].В качестве граничных условий на боковых границах для уравнений (1) исполь-зовались простые градиентные условия для отклонения концентраций от фоновыхзначений (фоновые значения рассчитывались на основе боксовой фотохимиче-ской модели без учета источников примеси). На верхней границе также применя-лись простые градиентные условия, а на нижней учитывалось сухое осаждениекомпонент примеси за счет их взаимодействия с подстилающей поверхностью [2].Начальные условия для уравнений (1) соответствовали фоновым значениям кон-центрации либо использовались результаты прогноза за предшествующий периодвремени.Задание параметров атмосферного пограничного слояДля реализации модели переноса примеси необходимо определить для каждо-го момента времени профили компонент вектора скорости, температуры и турбу-лентных характеристик. Такого рода данные можно получить усвоением резуль-татов 48-часового метеорологического прогноза, выполняемого с использованиемглобальной модели ПЛАВ. ПЛАВ - это полулагранжева модель атмосферы с вы-соким разрешением (горизонтальная сетка 0,72° - по широте, 0,9° - по долготе, 28вертикальных слоев), принятая Гидрометцентром в качестве оперативной моделидля среднесрочного прогноза с выдачей результатов через 6 часов. В модельПЛАВ включены параметризации процессов подсеточного масштаба (коротко- идлинноволновая радиация, глубокая и мелкая конвекция, планетарный погранич-ный слой, торможение гравитационных волн, тепло- и влагообмен с подстилаю-щей поверхностью с учетом растительности), разработанные Метео-Франс и ев-ропейским консорциумом LACE (www.rclace.eu) для оперативной региональноймодели ALADIN. Прогноз полей метеорологических элементов осуществляется спомощью численного решения уравнений = − , (5)где U,V - компоненты горизонтальной скорости ветра в атмосферном погранич-ном слое (W ≅ 0), Q - влажность, 〈uw〉, 〈vw〉, 〈ƒw〉, 〈qw〉 - осредненные по Рей-нольдсу корреляции турбулентных пульсаций вертикальной скорости и компо-нент горизонтальной скорости, температуры и влажности соответственно.( g,g) 1 ,U V p pf y x=ƒ ⎛⎜⎝−  ⎞⎟⎠- компоненты скорости геострофического ветра,f = 2ƒsinƒ - параметр Кориолиса, ƒ - географическая широта рассматриваемойточки, ƒ - угловая скорость вращения Земли, p - давление.Для замыкания системы уравнений (2) - (5) применяется трехпараметрическаямодель турбулентности, включающая уравнения переноса для энергии, масштабатурбулентных пульсаций и уравнение для дисперсии турбулентных пульсацийпотенциальной температуры 〈ƒ2〉 [3]:32D ;ek uw U vw V g w kl k C kt z z z z l = −  −  + ƒ ƒ +  ⎛⎜⎝ƒ  ⎞⎟⎠−(6)2l1 e l2 1 ;l C uw U vw V gw l kl l C k lt z z k z z z = ⎜⎝⎛−  −  +ƒ ƒ⎟⎠⎞ + ⎜⎝⎛ƒ  ⎟⎠⎞+ ⎡⎢⎣ −⎜⎝⎛ ƒ ⎟⎠⎞ ⎤⎥⎦(7)2 2 2C kl w2 2w 2 .t z z z ƒƒ〈ƒ 〉 ⎛ 〈ƒ 〉⎞ ƒ ƒ = + ⎜⎝ 〈 〉  ⎟⎠− 〈 ƒ〉 − ƒ(8)Здесь k = 0,5(〈u2〉+〈v2〉+〈w2〉) - кинетическая энергия турбулентности; l - инте-гральный масштаб турбулентности, ƒe = 0,54; Cl1 = −0,12; Cl2 = 0,2; CD = 0,19;Cƒ = 3,0 - числовые коэффициенты, τθ - временной масштаб температурных тур-булентных пульсаций.На верхней границе для определяемых из уравнений (2) - (8) зависимых пере-менных использовались простые градиентные условия (равенство нулю произ-водной по z, для всех переменных, для потенциальной температуры - равенствонулю второй производной), на нижней границе применялись соотношения теорииподобия Монина - Обухова и эмпирические функции параметров подобия [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]smog0; R⎡⎣R ⎤⎦=R[APM]=FCH2ƒk1CRsmog +FHNO3k7CRPCNO2 +FH2SO4k8CRPCSO2 ;R[RP]=k1CRsmog−k2CRPCNO−k5CRPCRP−k6CRPCNO2−k7CRPCNO2−k8CRPCSO2 ;R[H2O2] = ƒk5CRPCRP ;R[NO+NO2] = −k6CRPCNO2− k7CRPCNO2;R[NO2]=k2CRPCNO+k4CNOCO3−k3CNO2−k6CRPCNO2 −k7CRPCNO2 ;R[O3+NO2]=k2CRPCNO−k6CRPCNO2− k7CRPCNO2;R[SO2] = −k8CRPCSO2 ;FCH2 = 0,57 , FH2SO4 = 4,0 , FHNO3 = 2,6 - приближенные коэффициенты, оп-ределяющие переход устойчивых негазообразных устойчивых соединений в аэро-золи, ƒ= 0,1, [ ]smog2max 0,03,exp 0,0261NO+NOƒ = ⎛⎜⎜ ⎛⎜⎜− ⎡⎣R ⎤⎦ ⎞⎟⎟ ⎞⎟⎟⎝ ⎝ ⎠ ⎠- численные коэффициен-ты [2, 6]. Для расчета компонента примеси CO, который принимался в даннойработе химически инертным ( R[CO] = 0; ), также используется уравнение перено-са (1).Численный метод и его параллельная реализацияУравнение (1) с учетом представленной выше схемы замыкания можно запи-сать в следующем виде (ƒ = Ci, S = Si, R=Ri, ƒ = ƒi, i - фиксированное значение от1 до ns):.xx xy xzyx yy yz zx zy zzU V Wt x y z x x y zS Ry x y z z x y zƒ+ƒ+ƒ+ƒ=⎛⎜⎝ƒ ƒ+ƒ ƒ+ƒ ƒ⎞⎟⎠++⎛⎜⎝ƒ ƒ + ƒ ƒ + ƒ ƒ⎞⎟⎠+ ⎛⎜⎝ ƒ ƒ + ƒ ƒ + ƒ ƒ ⎞⎟⎠ − ƒƒ + +(9)При построении разностной схемы для уравнения (9) использовалась сетка спеременным по z шагом, сгущающаяся к поверхности Земли. При замене диффе-ренциального уравнения (9) его разностным аналогом применялись явно-неявныеразностные схемы с первым порядком аппроксимации по времени и вторым покоординатам. Все члены уравнения (9) аппроксимируются на k-м слое по времени,за исключением вертикальной диффузии и членов, отвечающих за влажное осаж-дение примеси, поступление примеси от источников и химические реакции, кото-рые берутся на k+1-м слое по времени.Для получения конечно-разностного аналога уравнения (9) используется методконечных объемов. Схема расчетной ячейки представлена на рис. 1. Для аппрок-симации диффузионных потоков используется центрально-разностный оператор.Для аппроксимации адвективных членов применяется направленная противопо-токовая схема MLU [7].Рис. 1. Вычислительная ячейка с обозначением центровиспользуемых по шаблону ячеек и граней выбранногоконечного объемаДискретный аналог уравнения переноса (9) с учетом принятых аппроксимацийможно записать в виде0 (1 1 ) 1 1 1,k k k k kP P P P P t b P t T b Bk k k k k k kE E W W N N S S T T B B P Pa xyz S R xyz D D D Da a a a a a a b⎡⎣ +ƒƒ ƒ ƒ − + ƒ ƒ ƒ + + ⎤⎦ƒ + = ƒ + + ƒ + ++⎡⎣ƒ + ƒ + ƒ + ƒ + ƒ + ƒ − ƒ⎤⎦+ (10)где aE=max(−Fe;0)+De ; aW=max(Fw;0)+Dw ;aN=max(−Fn;0)+Dn ; aS=max(Fs;0)+Ds ;aT=max(−Ft;0); aB= max(Fb;0); oaP =aE+aW+aN+aS+aT+aB−aP;( 0k 0k) *b= SP +RP ƒxƒyƒzP+DmP ; 0 PPx y zatƒ ƒ ƒ=ƒ.Здесь k 1 ( , , , )ƒP+ ≈ƒtk+ƒtxl yjzP , k - номер слоя по времени, ƒx, ƒy, ƒzp - размерыконечного объема с центром в точке P; центры соседних по горизонтали конеч-ных объемов обозначаются по направлениям сторон света - N, S, E, W, а по верти-кали - T и B для нижнего и верхнего соответственно; грани конечного объема - n,s, e, w, t, b; F, D представляют адвективные и диффузионные потоки через соот-ветствующие грани конечного объема; при аппроксимации источниковых членовиспользовалась «линеаризованная» форма записи: : k1 0k 1k k1SP+ ≅SP +SPƒP+ ,k1 0k 1k k1RP+ ≅RP +RPƒP+; 1k, 1k 0RP SP ≤ ; Dmp* содержит смешанные производные навременном слое k и наклоны в MLU-представлении потоков.Используемая явно-неявная разностная схема (10) позволяет неявным образомпровести расчет диффузионных турбулентных потоков вдоль оси Oz и тем самымизбежать существенных ограничений на шаг по времени вследствие малых верти-кальных размеров вычислительных ячеек, высота которых из-за сгущения вычис-лительной сетки к поверхности в первом узле составляет всего 2 м. В то же времягоризонтальные размеры вычислительной сетки (500 м) позволяют эффективноиспользовать явную вычислительную схему. Такой подход обеспечивает высокуюскорость вычислений вследствие применения экономичного метода прогонки длярешения на каждом шаге по времени (ƒt = 6 с) систем линейных уравнений стрехдиагональной матрицей (9).Расчетная область - равномерная в горизонтальных направлениях (Ox и Oy) инеравномерная, сгущающаяся к поверхности конечно-разностная сетка размеромNx .Ny .Nz =100.100.100 узлов.При моделировании переноса примеси используются среднегодовые данныеантропогенной эмиссии загрязнителей городского воздуха. Источники делятся на3 категории: точечные, линейные (дороги) и площадные (крупные предприятия).Для каждой из ячеек поверхности задается расход выбросов 17-ти различных посоставу химических веществ.Параллельная версия алгоритма была построена с использованием двумернойдекомпозиции сеточной области с перекрытием (рис. 2). Декомпозиция осущест-влялась по направлениям Ох и Оу.xyzРис. 2. Схема двумерной декомпозиции параллельной реализацииявно-неявного метода на примере четырех процессоров (p=4)На рисунке присутствуют светлые блоки по краям расчетной области, являю-щиеся фиктивными границами, темные блоки, расположенные со стороны разде-ления области, являются блоками обмена данными между процессорами. При та-кой декомпозиции можно использовать большее количество процессоров супер-компьютера, а доля временных затрат, связанных с необходимостью на каждомшаге по времени передавать данные между процессорными элементами, меньше,чем при одномерной декомпозиции. Расчеты, выполненные на кластере ТГУСКИФ Cyberia[8], показали, что при использовании 64 процессоров ускорениевычислений составляет 54, а эффективность соответственно 0,84. При этом одинчас моделирования на 64 процессорах выполняется в среднем за 20 секунд, в товремя как время расчета этой же задачи на одном процессоре - 18 минут.Сравнение результатов моделирования с данными измеренийПроверка полученных результатов моделирования осуществлялась путем срав-нения расчетов с данными измерений, проведенными на ТОР-станции Институтаоптики атмосферы СО РАН в Томске. Расчетная область представляет собой частьатмосферы размером 50.50.2 км, выбранная таким образом, что ее квадратное ос-нование содержит участок подстилающей поверхности с крупным населенным16 Апреля 2010 8 Декабря 2009 г. 0 г.Скоростьветра, м/с86420Направлениеветра, град.3602401200КонцентрацияCO, мг/м310.80.60.40.20КонцентрацияNO , мкг/м2 3100806040200КонцентрацияO , мкг/м3 3100806040200Время, ч0 4 8 12 16 20 24Время, ч0 4 8 12 16 20 240 4 8 12 16 20 24 0 4 8 12 16 20 240 4 8 12 16 20 24 0 4 8 12 16 20 240 4 8 12 16 20 24 0 4 8 12 16 20 240 4 8 12 16 20 24 0 4 8 12 16 20 24Рис. 3. Сравнение результатов моделирования с данными TOР-станции8 декабря 2009 г. и 16 апреля 2010 г.пунктом в центре. В расчете учитывается эмиссия примеси от линейных, площад-ных и точечных источников, находящихся на поверхности или приподнятых надземлей. Область покрывается вычислительной сеткой с размерами 100.100.100 уз-лов. Период прогностического моделирования обычно составляет 24 часа.На рис. 3 сплошной линией представлены результаты численного моделирова-ния, символами отмечены результаты измерений, проводимых на ТОР-станцииИОА, для каждого измерения вдоль оси значений приведена погрешность.На графиках дается сравнение измерений и расчетов силы и направления ветра,а также концентрации таких загрязнителей, как монооксид углерода (CO), диок-сид азота (NO2) и озон (O3). Из графиков сравнения видно, что математическаямодель хорошо предсказывает изменение ветра и значений рассматриваемых кон-центраций в течение суток.ЗаключениеПредставленная в работе модель имеет ряд преимуществ, позволяющих ис-пользовать ее в реальных условиях в оперативном режиме. Такая возможностьпоявилась благодаря усвоению данных глобального прогноза метеорологическиххарактеристик, а также за счёт параллельной реализации модели переноса приме-си при помощи двумерной декомпозиции расчетной области. В настоящее времямодель реализована и является ядром программного комплекса, расположенногона выделенном сервере. Ежесуточно программными компонентами комплекса вавтоматическом режиме производится получение начальных данных, их обработ-ка, запуск процесса моделирования на кластере ТГУ СКИФ Cyberia [8] и обработ-ка полученных результатов с последующим их представлением в сети Internet ввиде полей прогноза распространения в течение суток примеси, наложенных накарту местности[9].

Скачать электронную версию публикации