CC BY
162
ISSN 1992-6502 (P ri nt)_
2015. Т. 19, № 1 (67). С. 213-223
Ъыьмт QjrAQnQj
ISSN 2225-2789 (Online) http://journal.ugatu.ac.ru
УДК 519.86
Математическая модель распространения примеси в приземном слое атмосферы и ее программная реализация на многопроцессорной вычислительной системе
а. и. Сухинов , д. с. Хачунц , а. е. Чистяков
1 [email protected], 2 [email protected], 3 [email protected]
1 Таганрогский институт имени А. П. Чехова (филиал) ФГБОУ ВПО «Ростовский государственный экономический университет» (РИНХ) 2,3 ФГАОУ ВО «Южный федеральный университет» (ЮФУ)
Поступила в редакцию 1 февраля 2015 г.
Аннотация. Целью работы является разработка математической модели движения многокомпонентной воздушной среды, учитывающей транспорт загрязняющих веществ (ЗВ) и тепла, фазовые переходы, а также влияние растительного покрова (лесных насаждений) на распространение ЗВ в атмосфере; построение и исследование разностных схем, аппроксимирующих исходную задачу; программная реализация разработанных параллельных алгоритмов и проведение численных экспериментов по моделированию движения многокомпонентной воздушной среды, в том числе применительно к рекреационной среде прибрежных систем. Поскольку задачи прогнозирования распространения загрязняющих веществ должны решаться в реальном или ускоренном масштабах времени, на сетках, включающих 106-109 узлов, необходима разработка параллельных алгоритмов задач гидродинамики на системах с массовым параллелизмом. Предложена новая математическая модель аэродинамических процессов, учитывающая повышенную влажность воздуха, изменчивость атмосферного давления и температуры и др., что особенно свойственно прибрежным районам. Разработаны параллельные алгоритмы исследования этих моделей, реализованные в виде комплекса программ.
Ключевые слова: воздушная среда; двумерная дискретная модель; загрязняющие вещества.
ВВЕДЕНИЕ
Одним из практически важных разделов наук, связанных с изучением математического моделирования, является изучение процесса распространения примесей в приземном слое атмосферы. Математическое моделирование является надежным и эффективным методом анализа и оценки состояния воздушной среды, которое позволяет определить распределение концен-
Статья рекомендована к публикации программным комитетом международной научной конференции «Параллельные вычислительные технологии 2015».
Работа выполнена при частичной поддержке проектов Программы № 43 фундаментальных исследований Президиума РАН по стратегическим направлениям развития науки «Фундаментальные проблемы математического моделирования».
трации примеси в атмосфере, учитывающее многокомпонентный характер среды.
В области математического моделирования движения загрязнений в атмосфере и разработки численных методов для него, в настоящее время, сложилась ситуация, при которой проводимые в мире работы рассматривают отдельные явления и не охватывают их комплекса. Работы, посвященные математическому моделированию приземной аэродинамики, впервые в Советском Союзе появились у академика А. С. Монина [1]. Далее в Ленинграде профессор М. Е. Берлянд перехватил научную эстафету. Более поздние исследования такого типа начались в Новосибирске академиком Г. И. Марчуком [2, 3]. Рассеяние примеси, которое учитывает турбулентность, и к ним примыкают дополнительные условия (распространение тепла, влажности и т. д.), которые потом начинает воздействовать как единая система. Конечно, эта система -сложная, поэтому требует привлечения числен-
ных методов. Для решения проблем, отвечающие поставленной задаче, необходима разработка новых математических моделей, которые базируются на уравнениях газовой динамики и законах сохранения вещества.
Однако сложность и взаимосвязанность процессов распространения и рассеяния примеси, происходящих в турбулентном атмосферном пограничном слое, делают модели прогнозирования качества воздуха громоздкими в математической записи и весьма требовательными к вычислительным ресурсам. Перспективным способом решения этих проблем является применение эффективных численных схем высокого порядка точности и использование компьютеров с параллельной архитектурой при проведении вычислений. Параллельному программированию и методам распараллеливания численных алгоритмов посвящены работы В. В. Воеводина [4, 5], Б. Н. Четверушкина [6], Ю. И. Шокина, В. П. Гергеля [4, 5], В. И. Ма-лышкина, И. Б. Петрова [7], М. В. Якобовского [8], В. А. Вшивкова, D. Dubdub и др.
ОСНОВНЫЕ УРАВНЕНИЯ МОДЕЛИ ПРИЗЕМНОЙ АЭРОДИНАМИКИ
Атмосфера представляет собой сложную динамическую систему, в которой протекают различные динамические и физико-химические процессы. Эти процессы обусловлены как атмосферной циркуляцией, так и трансформацией газовых и аэрозольных примесей.
Движение воздушной среды и распространение в ней загрязняющих веществ проходит в четыре стадии [9]: модель движения многокомпонентной воздушной среды, которая предназначена для поля вектора скорости, учитывающая турбулентный обмен, переменную плотность, зависимость плотности воздушной среды от давления; модель, описывающая процессы переноса примеси, учитывающая переход воды из жидкого в газообразное состояние и наоборот, осаждение вещества, модель притока тепла, которая описывается уравнениями теплопроводности газа и конденсата, учитывающая теплообмен между жидкими и газообразными состояниями и транспорт тепла; модель расчета давления, учитывающая сжимаемость среды, тепловое расширение, источники вещества, связанные с переходом воды из жидкого состояния в газообразное и обратно, а также турбулентное перемешивание многокомпонентной воздушной среды.
Исходными уравнениями модели движения многокомпонентной воздушной среды являются
уравнение движения:
/ Лг = -(1 / р)ф / дх. + div(^гаа(у^))
ёг
» уравнение транспорта вещества [10]: др / Ы + с1гу(ру) = с1гу([1§гас1(р)) + /р;
» уравнение состояния (аналог уравнение Менделеева-Клапейрона):
Р = Х РДТ / М ;
г
уравнение транспорта примеси:
Лщ / Лг = div (^гаа (ф)) -
• уравнение притока тепла:
Лд / Лг = &у (д)) + div (^гаа (Т)) + /е;
• уравнение модели турбулентности:
^ = (С*Л)2 ^>
где ф - объемные доли 1-й фазы (/ = 0 - воздух,
1 - вода в газообразном состоянии, 2 - газ на источнике, 3 - вода в жидком состоянии, 4 -сажа); - проекции компонентов скорости на оси Оxj, j = 1, 2, 3; g - ускорение свободного падения; р - давление; Я - универсальная газовая постоянная; М - молярная масса; I - функция, описывающая распределение и мощность источников примесей; Т - температура газовой фазы; 2 - тепловая энергия, X - коэффициент теплопроводности; ц - коэффициент турбулентного обмена; р - плотность воздушной среды. Здесь и далее символ «^^» означает
полную производную по времени от функции, зависящей в общем случае от времени ^ и трех пространственных координат х, у, г.
Переход от трехмерных моделей к двумерным. С целью получения нетрудоемких для вычисления дискретных моделей осуществляется переход от трехмерных моделей к двумерным. Рассмотрим трехмерное уравнение диффузии-конвекции-реакции:
др д(ри) д^) д(р^)
дг дх ду дг
дх
др
+ -
д
дх) ду 1 дУ
др
Л
д_
дг
др
+/
Будем рассматривать уравнение при следующем граничном условии: рцр^ = -т, где т -параметр, описывающий наличие источника на боковых поверхностях; п - вектор нормали.
В результате преобразований получим:
+ д(ври) | д(ер^_ дГ дрр |
дг дх дг дх I дх)
8 ( 8р ^ т 82 ^ 82) р
+ е/„
где е - параметр, описывающий относительную величину объема, свободного от растений, (е, п ) - коэффициент турбулентного обмена, зависящий от проницаемости и видового состава лесного насаждения, который задается параметром щ I = 1, 2, ..., Ь, Ь - общее количество видов, составляющих данное насаждение.
ДВУМЕРНАЯ ДИСКРЕТНАЯ МОДЕЛЬ
РАСПРОСТРАНЕНИЯ ПРИМЕСИ В ПРИЗЕМНОМ СЛОЕ АТМОСФЕРЫ
Математическая модель аэродинамики приземного слоя атмосферы. Рассмотрим основные уравнения динамики воздушной среды:
• система уравнений Навье-Стокса
еи, + иеи, + уеи , =
: =-1 (еР)',
+
р
еж + иеч + чем.
+(РеиХ )х + (реи2 )г + е/х,
' =-1 (ер у, +
+()х + (ре< )г + е/; О)
• уравнение неразрывности ер;+(ери)' х+(ер^) ' г =
= (ерр :) ' х +(ерр:) ' г + е/р; (2)
уравнение состояния
р = Ур-яг,
; М
(3)
где е - параметр, описывающий относительную величину объема, свободного от растений; р = р (е, п ) - коэффициент турбулентного обмена, зависящий от проницаемости и видового состава лесного насаждения, который задается параметром 1=1,2,.,Ь; Ь - общее количество видов, составляющих данное насаждение.
Предположив, что воздушная среда, изначально, находится в состоянии покоя, начальные условия будут иметь вид: и = 0, ч = 0, Р = р , где V = {и, ы} - значение компонент вектора скорости, р - атмосферное давление.
Система уравнений (1), (2) рассматриваются при следующих граничных условиях: • на непроницаемой границе: рж = Т;, (г) р^П = тг, , ^ = 0, Р = о;
на боковых проницаемых границах:
и: = 0, Ч = о, Р = о;
и = и, ч = Ж, Р = 0
• на источнике: где V = {и, у} - значение компонент вектора
скорости; Р - давление; и, Ж- компоненты вектора скорости на источнике; р - плотность;
(е,: ) - коэффициент турбулентного обмена, зависящий от проницаемости и видового состава лесного насаждения; М - молярная масса; Я - универсальная газовая постоянная, тх, тг - составляющие касательного тангенциального напряжения; Т - температура.
Согласно методу поправки к давлению [1, 12] исходная модель гидродинамики (1), (2) разбивается на три подзадачи. Первая подзадача представлена уравнением диффузии-конвекции-реакции, на основе которого рассчитываются компоненты поля скорости на промежуточном слое по времени й — и ' '
к
■ +
иггТх + м>гТ?2 =(\ю?х)х
\v-\v
Для аппроксимации по временной переменной уравнения диффузии-конвекции-реакции использованы схемы с весами. Здесь й = ай + (1 -о)и, ое[0,1] - вес схемы [13].
Опишем граничные условия системы (4): на непроницаемой границе:
рЖ = т;,ьОХ рЖ = тг,ь(0, на боковых проницаемых границах: и'п = 0, м>'п = 0, на источнике: и = и, ч = Ж, р = 0. Вторая подзадача позволяет рассчитать распределение давлений
о——— = + г^р;: + о——— - (врй)' -
р а » » , т ы \ V ) х
-(вР#)'г+(вцР: )'х+(вцР: )'г+в1Р.
(5)
Третья подзадача позволяет по явным формулам определить распределение скоростей на верхнем временном слое
й — й 1
Ь -ЦгР)' , г^ = -1-(гР)' , (6) к р ; к р где к - шаг по временной координате; и - значение поля скорости на предыдущем слое по времени; й - значение поля скорости на промежуточном слое по времени; и - на текущем слое по времени.
Следует отметить, что при расчете давления учитывается сжимаемость среды, тепловое расширение, источники вещества, связанные с пе-
е
реходом воды из жидкого состояния в газообразное и обратно, а также турбулентное перемешивание многокомпонентной воздушной среды. При расчете полей плотности вещества и температур источники, вызванные расширением (сжимаемостью) среды, присутствует в операторе конвективного переноса.
Аппроксимация задачи аэродинамики приземного слоя атмосферы. Для построения дискретной модели движения воздушной среды, рассмотрим расчетную область, которая вписана в прямоугольник. Покроем область равномерной прямоугольной расчетной сеткой:
ш,
ш.
= (гя = пК,0< п < N -1,Ь = К(N-1)},
= (х, = 1Их, 0 < г < N -1,1х = Их (Nx -1)} ,
шг = (г. = ,0 < ] < N -1, ¡2 = К (N -1)}, где п, г, . - индексы по временной координате и пространственным координатным направлениям Ох , Ог соответственно; К, К, К - шаги по временной координате и пространственным координатным направлениям Ох, Ог соответственно; N, N, N - количество узлов по временной координате и пространственным координатным направлениям Ох, Ог соответственно; I, К, К - длина расчетной области по временной координате и пространственным координатным направлениям Ох, Ог соответственно.
Для аппроксимации задачи движения воздушной среды будем учитывать заполненность ячеек. В связи со сложной геометрией расчетной области учет заполненности ячеек повысит точность решения.
Дискретные аналоги операторов конвективного ие'х и диффузионного (ЦсП)^. переноса
в случае частичной заполненности ячеек с граничными условиями третьего рода с'п (х, г, г) = апс + Ри могут быть записаны в следующем виде [14, 15]:
( 4о), ,ис'х = ( 41)
с , . — с .
1+1, ] г,) и иг+\12,3 +
2И
+ ( 42 ),, }иг -1/2,}
с. . - с. , .
',1 '-1,1
2И
(4о),, (Цсх )Пх = (41, Ц,+1/2,
с. , . - с. .
'+1,1 ', ]
К2
-(42),/ ц
с. . - с. , .
1 ,-1,.
1 ^-1/2,1
где я , я , я - коэффициенты заполненности контрольных областей.
Аппроксимация первого уравнения системы (4) запишется в виде [16, 17]:
( 4о),
А,
(4:)
,, 1 иг+1/2,1
и. , . - и. .
Ж
( 42 )
и. . — и._х . 2 Л,- ' ^ + (4з )
, 1-2,1 Иг
г, +12
и. -и. . 2И
+ ( 44 )
и. . - и.
и.,, . -и.
4 ^^ -1/2 2К
-(42),, Ц--( 44 )1., 1 Ц,
1/2,./
и и-1 =( 41X, 1 Ц,+1/2,1
и • ',1+1/2 К2
и . -и. , .
г,1 г-1,1
+ (4з),, Ц,
и. -и. .
и. . - и. . ,
!-1/2 7.2
((4з)„-(44)„.)(-0„ •
(7)
Параметр щ принимает значение 1 в случае, если узел (', 1) принадлежит множеству граничных узлов, находящихся в придонной области, в противном случае щ = 0. Аналогичным образом получим дискретный аналог для второго уравнения системы.
С учетом несложных выкладок аппроксимация первого уравнения системы (6) примет вид:
( 4о),
ии ~ии
1 К
Р - Р
1,1 2кхр,
Р - Р -(42),, "
1,1 ^р.
(8)
Для решения сеточных уравнений используем адаптивный попеременно-треугольный итерационный метод.
Аппроксимация уравнения (5) примет вид:
р Р - Р
(4о00 =к
и* Р К
Р - Р
(Я),, "
К
Р - Р Р - Р Р - Р
-(42 + (4, -(Я4
Т, - Т, (41>,^(р" ),„r.j-М,! (ри ),-уХ1
г,1 г,1 г,1
т К
г,1 <
К
( 4з ) (р,,; +1/2 -( 44),., (р
)г,}-1/2
, тЛ (42),,1 -(41 ),,1 (44X, 1 -(4з),,1
-(ри ). ,1-^-1 - (рЖ ). 1-^-1 +
(41),,; -(42),,
а хс,, 1 + Р х
К
ш = ш х ш х ш :
г
+
К
г
К
х
+ (41 ),,, Цг'+1/2,
р'+11 -р,,1
р
+1 Г'',.;
-(42 ), Ц,--(44),,Ц
+ (43),,, Цу+1/2" (44),,, Цг', j-'/2
+(4о),,/р, (9)
где {и, Ж} - значение компонентов вектора скорости в граничных узлах расчетной области. В случае непроницаемой границы параметры и и Ж равны нулю. В случае проницаемой границы параметры и, Ж равны и, w соответственно. Также компоненты и, Ж могут быть заданы.
Дискретная модель транспорта ЗВ в воздушной среде. Для описания процессов переноса ЗВ в работе используются уравнения, учитывающие такие факторы, как переход воды из жидкого состояния в газообразное и обратно, а также в процессе транспорта примесей взвешенные частицы осаждаются.
Запишем систему уравнений, описывающую процессы переноса примеси в многокомпонентной воздушной среде, в следующем виде [16]:
д Г ,дФо
е ^ + и е ^ + * е дфо
дг
дх дФо
д Г
+—I Це
дг I дг
дг дх
г = о, 2;
це-
дх
+
(10)
дф дф дф д е—1— + и е—— + * е—— = —
дг дх дг дх
д Г дф ^ ^ +¥Г Це аГ)+7.8;
дф дф дф, е—112 + ие—112 + м>е-
д
дг
+-
дг
Це
дх
дф2 дг
дг дх
Це 1Г ' +
дх
Це дф2, +
(11)
дх
+ /е - е (1 - е) а2 ф2;
(12)
д
дф дф дф д е—3 + ие—3 + (* - ) е—3 = —
дг дх о дг дх
Ц£
дфз дх
+
д Г дф,
+—I це—-
дг 1 дг
|—- е - е (1 - е) а3ф3; р,
дф4 дф4 , ч дф4 д г
е—114 + ие—+ (* - ) е—114 = — дг дх дг дх
(13)
це
дф4
дх
+
+-
&
це■дф± |-е(1 -е)а4ф4,
V
дг
(14)
где V = / (ря - р1) - массовая скорость испарения; е - параметр, описывающий проницаемость; ц = ц (е,. ) - коэффициент турбулентного обмена, зависящий от проницаемости и видо-
вого состава лесного насаждения; рп - плотность насыщенных паров; ф - объемные доли 7-й фазы (7 = 0 - воздух, 1 - вода в газообразном состоянии, 2 - газ на источнике, 3 - вода в жидком состоянии, 4 - сажа); щ - коэффициент
скорости абсорбции.
Соответствующие системе уравнений (10)-(14) граничные условия имеют вид
(ф' )'п (г, х, г)
= 0.
(х, г )еГ
Уравнение для расчета скорости абсорбции имеет вид
да, дг
= в1,, (ао - а,' ) - Р2,г'ф,а,' , ■ = 2,33,44
где Р1 , - коэффициент, описывающий рост скорости абсорбции; Р21 - коэффициент, описывающий уменьшение скорости абсорбции за счет насыщения; а0 - максимальное значение коэффициента скорости абсорбции; щ - коэффициент скорости абсорбции.
Дискретный аналог системы уравнений (10)-(14) транспорта ЗВ на промежуточном временном слое с учетом частичной заполенно-сти ячеек будет иметь вид
«п+1 п п+о я+о
(4 \ ф,+ -фп1 +и ф+Ц -ф,+ (4 \ + (4оЛ, К +и'+1/2,1 2К (41Л, +
„я+о
и+о
+ы.. --(</.-).. +^+1/2Фц+' +
i-'/2,]
2К
2К
Фц ~Фу-1/ ч _ ° ° ---(Я ) -
-1/2
2К
и+о ,лп+о
п+о _ п+о
II ^ ф-7л ^ III фi,]^+1 ^ ^
-Цг-■-~2-(42 А, + Ц+т-~2-(43 -
ф',]+' -', 1
—(42)',, + Цi^ь/+'/2
+о ^п+о ф^ - ф^-1 ^ Ч ^ , ^ Ч
"Ц,,!-1/2--( 44),', 1 + ¿,1 ( 4о),', 1
(15)
где м> = м> для уравнений (10)—(12) и для (13), (14), о е [0,1] - вес схемы.
Исследование консервативности модели движения многокомпонентной воздушной среды и устойчивости модели транспорта примеси. Выполнена проверка устойчивости разностных схем, предназначенных для математического моделирования транспорта ЗВ. Для исследования устойчивости был использован сеточный принцип максимума [18].
Таким образом, получены: • ограничения на шаг по пространству: Их <2т1и(ц/\и\);
д
Рис. 1. Поле распределения температуры
• ограничения на шаг по времени (получено из условия положительности коэффициентов сеточных уравнений согласно условиям применимости принципа максимума):
И, < ((1 - о)(4шт (р) / И] - (и + ^) / 2))-1.
Для функций концентрации имеют место
:
оценки: Цф^Ц <1 \фЧ + к£11 /к\ , г = 1..4 .
Для модели движения многокомпонентной воздушной среды справедливо балансовое соотношение:
Т . Т .
£ (ч0 ), р,= £ (ч0 )<,,+
ге[1,^;-1] т ¡, - ¡е[1,N-1] т ¡,-
-е[1, N-1] -е[1, N-1]
г (ч2) -(ч)
+к £ (40),Л-(ри) -(Ч2л' (41 )г-'
N -1] . /е[1,Nz -1] v
h
-(рж)., - л
2 У
Из полученного выражения следует, что количество вещества £ (ч0) р i - на
¿е[1,^; -1] ,
-е[1,Кг -1]
текущем временном слое зависит от количества вещества на предыдущем временном слое, от интенсивности внутренних и граничных источников, и с точностью до членов порядка О ( И2 + т ) выполняется закон сохранения и превращения вещества.
СОПОСТАВЛЕНИЕ РЕЗУЛЬТАТОВ ЧИСЛЕННЫХ ЭКСПЕРИМЕНТОВ И НАТУРНЫХ ИЗМЕРЕНИЙ ПО МОДЕЛИРОВАНИЮ ДВИЖЕНИЯ АЭРОЗОЛИ
Описание натурного эксперимента. Цель проведения натурного эксперимента заключает-
ся в сравнении исходных данных с результатами работы разработанного программного комплекса «AeroEcology», моделирующего физические процессы, рассматриваемые в работе. При проведении эксперимента измеряется распределение поля температуры нагретой аэрозоли от источника при помощи тепловизора. На рис. 1 приведены две картины распределения температуры, полученные экспериментальным путем.
Для обработки результатов натурных измерений использована гипотеза о том, что поле температуры распределено по нормальному закону. При этом для описания распределения температуры необходимо иметь информацию о следующих параметрах, зависящих от вертикальной координаты: максимальное и минимальное значении температуры, расположение центра трубки тока аэрозоли (описывает математическое ожидание), ширину трубки тока (описывает среднее квадратичное отклонение).
На рис. 2 представлено поле температуры со средними статистическими параметрами среди всех проведенных наблюдений: температуры в центре трубки тока, отклонение центра трубки тока от вертикали, ширины трубки тока. Поле температуры, приведенное на рис. 2, являлось эталоном при верификации численной модели.
Описание логической структуры разработанного программного комплекса «Aero-Ecology» и его параллельная реализация.
На основе построенных алгоритмов был создан комплекс программ «AeroEcology», предназначенный для численного моделирования движения многокомпонентной воздушной среды. Для верификации построенного комплекса программ «AeroEcology» был проведен численный эксперимент и выполнено сопоставление натурных данных с результатами численных экспериментов.
k=0
Рис. 2. Поле температуры со средними статистическими параметрами
Рис. 3. Поле температуры аэрозоли
При решении задачи были использованы расчетные сетки размерами 300 х 200 . Шаги по пространственным переменным равны 1 см., скорость движения воздушной среды на левой границе задавалось равной 1,4 м/с. Для решения задачи применены схемы с весами, при этом вес схемы задавался равным 0,5. Шаг по временной переменной равен 0,001 с, расчетный временной интервал составлял 5 с. На рис. 3 палитрой показана температура аэрозоли.
При моделировании транспорта аэрозоли (результаты численных экспериментов приведены на рис. 3) в граничных условиях использованы параметры, полученные на основе натурных экспериментов (рис. 2).
Для решения задачи транспорта многокомпонентной воздушной среды использован адаптивный модифицированный попеременно-треугольный метод (МПТМ) минимальных поправок. При параллельной реализации использованы методы декомпозиции сеточных областей для вычислительно трудоемких задач диффузии-конвекции, учитывающие архитектуру и параметры многопроцессорной вычислительной системы РОЦ НИТ ЮФУ. Пиковая производительность МВС составляет 18.8 TFlops. В качестве вычислительных узлов используется 128 однотипных 16-ядерных Blade-серверов HP ProLiant BL685c, каждый из которых оснащен
четырьмя 4-ядерными процессорами AMD Opteron 8356 2.3GHz и оперативной памятью в объеме 32ГБ. Схему итерационного двухслойного модифицированного попеременно-треугольного метода запишем в следующем виде [19]:
/+1 = У" - тn+w", (D + mR) D l (D + MR) w" = г", г" = Ay" - f,
где г" - вектор невязки; w" - вектор поправки.
Для параллельной реализации адаптивного МПТМ использованы методы декомпозиции области по одному направлению. Наиболее трудоемким расчетом с точки зрения построения параллельной программной реализации является расчет вектора поправки. Он выполняется в два шага:
1) (D + mR ) y" = г" ,
2) (D + mR ) w" = Dy" .
На первом шаге рассчитываются элементы вспомогательного вектора y" снизу вверх, а затем, зная вектор y , рассчитываются элементы вектора поправки w" сверху вниз.
Идея расчета вектора поправки представлена на рис. 4. Стрелками показано направление счета и передачи. В схеме для расчета вектора ym только первый процессор не требует дополнительной информации и может независимо от других процессоров посчитать свою часть области, остальные процессоры ждет от предыдущего, пока он не передаст элементы, стоящие вначале строки. Передача по одному элементу не оптимальна, т.к. появляются временные затраты, связанные с организацией передач, суммарное время на накладные расходы можно уменьшить путем увеличения объема передач. Данные рассуждения используются для расчета вектора поправки.
Результаты использования многопроцессорных технологий для расчета полей течений приведены в таблице.
Теоретически получены оценки ускорения и эффективности параллельного алгоритма, зависящие от времени выполнения арифметической операции, времени передачи данных и ла-тентности. Практически максимальное ускорение достигалась на 64 узлах и равнялась 39,5.
Таблица 1
Количество ядер Время, с Ускорение Эффективность
1 1447.41 1 1
5
2 734.728 1.97 0.985
4 387.009 3.74 0.935
8 199.643 7.25 0.906
16 109.653 13.2 0.825
32 62.659 23.1 0.722
64 36.643 39.5 0.617
Верификация математической модели приземной аэродинамики. Для верификации разработанной математической модели аэродинамики проведен ряд численных экспериментов и выполнено сопоставление результатов расчетов с расчетами, выполненными в среде
На рис. 5, 6 показана интенсивность движения воздушной среды, реализованных в среде разработанного программного комплекса «Aero-Ecology» и в среде ANSYS. Из увеличенных областей наблюдается качественное совпадение результатов расчета, что свидетельствует о ва-лидности комплекса программ «AeroEcology».
Результаты численного эксперимента по моделированию транспорта аэрозоли, приведенные на рис. 3, на качественном уровне совпадают с результатами натурных наблюдений (рис. 2). Ниже приведены результаты верификации математической модели и количественные оценки.
На рис. 7 показано сопоставление для следующих параметров: а) температуры в центре трубки тока; б) траектории движения аэрозоли; в) ширины трубки тока.
Рис. 4. Схема расчета вектора поправки (слева показана схема расчета вектора уп,
справа - вектора поправки)
Рис. 5. Интенсивность движения воздушной среды (м/с) при горизонтальном сечении области
Рис. 6. Расчет скорости движения воздушной среды в среде ANSYS, горизонтальное сечение
Рис. 7. Сопоставление численных экспериментов (кривая 1) с натурными данными (кривая 2)
Отклонение результата расчета от натурных данных составило: для температуры в центре тока 15%, для траектории движения аэрозоли 5% и для ширины трубки тока 9,7%.
МОДЕЛИРОВАНИЕ РАСПРОСТРАНЕНИЯ ПРИМЕСИ
В ПРИБРЕЖНОЙ ЗОНЕ ВОЗДУШНОЙ СРЕДЫ
Проведен вычислительный эксперимент по моделированию распространения примеси от автотранспорта в прибрежном районе города Таганрог. Исходными данными являются: плотность воздушной среды 1,29 кг/м3; плотность выброса 1,4 кг/м3; температура окружающей среды 20° С; температура выброса 120° С; скорость течения воздушной среды 1 м/с в направлении прибрежной зоны; удельная мощность выброса 5 л/с.
При решении модельной задачи были использованы расчетные сетки размерами 300 х 80. Шаги по пространственным переменным равны 1 м, скорость движения воздушной среды на левой границе задавалось равной 1 м/с. Для решения модельной задачи применены схемы с весами, при этом вес схемы задавался равным 0,5. Шаг по временной переменной равен 0,1 с, расчетный временной интервал составлял 300 с. На рис. 8 приведены результаты численного эксперимента по моделированию движения ЗВ без учета лесных насаждений.
На рис. 9 приведены результаты численного эксперимента по моделированию движения ЗВ при наличии лесных насаждений. Коэффициент проницаемости воздушной средой лесных насаждений составляет 50%. Высота насаждения составляет 30 м, ширина области насаждения - 50 м. При решении задачи использованы расчетные сетки размерами 300 х 80. Расчетный интервал составляет 300 с.
Сопоставив рис. 8, 9, можно сказать, что при наличии лесных насаждений наблюдается
более интенсивное перемешивание примесей в воздушной среде.
Рис. 8. Поле концентраций ЗВ через 300 с после начала выброса
Из рис. 8, 9 также виден подъем разогретых ЗВ, обладающих более низкой плотностью по сравнению с окружающей средой в районе источника. Однако ЗВ осаждаются за счет уменьшения их температуры и увеличения плотности, значительная часть которых попадает в зону отдыха.
ЗАКЛЮЧЕНИЕ
Разработана непрерывная двумерная математическая модель движения многокомпонентной воздушной среды, которая учитывает такие
Рис. 9. Поле концентраций ЗВ через 300 с после начала выброса
факторы, как транспорт ЗВ и тепла; влияние растительного покрова; изменение коэффициента турбулентного обмена; переход воды из жидкого в газообразное состояние; осаждение вещества; изменение температуры за счет конденсации и испарения аэрозоли; турбулентное перемешивание многокомпонентной воздушной среды; теплообмен между жидкими и газообразными состояниями; наличие распределенных источников вещества и температуры; силу Архимеда; тангенциальное напряжение на границах раздела сред; переменную плотность, зависящую от концентрации загрязняющих веществ, температуры и давления; сжимаемость среды за счет: изменения температуры, испарения и конденсации жидкости, изменения давления, наличия источников. Отличительной особенностью разработанной математической модели движения воздушной среды является учет влияния растительного покрова и турбулентного перемешивания в уравнении неразрывности среды.
Предложены консервативные разностные схемы для модели многокомпонентной воздушной среды, учитывающие такие физические процессы, как турбулентное перемешивание многокомпонентной воздушной среды, наличие распределенных источников вещества и температуры, силу Архимеда, тангенциальное напряжение на границах раздела сред, переменную плотность, сжимаемость среды, изменение давления, наличие источников. Выполнено исследование погрешности аппроксимации, устойчивости и консервативности разработанных разностных схем.
Разработан и реализован комплекс программ «AeroEcology», учитывающий процессы тепломассообмена в прибрежной зоне воздушной среды и предназначенный для построения турбулентных потоков поля скорости многокомпонентной воздушной среды на сетках с высокой разрешающей способностью, а также для расчета концентрации ЗВ и транспорта тепла. На основе разработанных моделей с использованием построенного комплекса программ «AeroEcology» выполнены численные эксперименты по прогнозированию распределения ЗВ в многокомпонентной воздушной среде прибрежной зоны г. Таганрога, результаты которых согласуются с натурными данными.
СПИСОК ЛИТЕРАТУРЫ
1. Монин А. С., Яглом А. М. О законах мелкомасштабных турбулентных движений жидкостей и газов // УМН. 1963. Т. 18, № 5 (113). С. 93-114. [[ A. S. Monin and A. M. Yaglom, "On the laws of small-scale turbulent flow of
liquids and gases," (in Russian), Russian Mathematical Surveys, vol. 18, no. 5 (113), pp. 93-114, 1963. ]]
2. Марчук Г. И. Математическое моделирование в проблеме окружающей среды. М.: Наука, 1982. [[ G. I. Marchuk, Mathematical modeling in environmental issues, (in Russian). Moscow: Nauka , 1982. ]]
3. Алоян А. Е. Динамика и кинетика газовых примесей и аэрозолей в атмосфере. Курс лекций. М.: ИВМ РАН, 2002. 201 с. [[ A. E. Aloyan, Dynamics and kinetics of trace gases and aerosols in the atmosphere, Lectures, (in Russian). Moscow: IVM RAN, 2002. ]]
4. Воеводин В. В., Гергель В. П. Суперкомпьютерное образование: третья составляющая суперкомпьютерных технологий // Вычислительные методы и программирование: новые вычислительные технологии. 2010. Т. 11, № 2. С. 117-122. [[ V. V. Voevodin, V. P. Gergel, "Supercomputer Education: The third component of supercomputer technologies," (in Russian), Vychislitel'nye Metody i Programmirovanie, vol. 11, no. 2, pp. 117-122, 2010. ]]
5. Антонов А. С., Артемьева И. Л., Буханов-ский А. В., Воеводин В. В., Гергель В. П., Демкин В. П., Коньков К. А., Крукиер Л. А., Попова Н. Н., Соколин-ский Л. Б., Сухинов А. И. Проект "Суперкомпьютерное образование": 2012 год // Вестник Нижегородского университета им. Н. И. Лобачевского. 2013. № 1-1. С. 12-16. [[ A. S. Antonov, I. L. Artemyeva, A. V. Bukhanovsky, Vl. V. Voevodin, V. P. Gergel, V. P. Demkin, K. A. Konkov, L. A.Krukier, N. N. Popova, L. B. Sokolinsky, A. I. Sukhinov, "Supercomputing education project: year 2012," (in Russian), Vestnik of Nizhny Novgorod University, no 1-1, pp. 12-16, 2013. ]]
6. Chetverushkin B., Gasilov V., Iakobovski M., Polyakov S., Kartasheva E., Boldarev A., Abalakin I., Minkin A. Unstructured mesh processing in parallel CFD project GIMM // Parallel Computational Fluid Dynamics 2005. Elsevier B. V., 2006. С. 501-508. [[ B. Chetverushkin, V. Gasilov, M. Iakobovski, S. Polyakov, E. Kartasheva, A. Boldarev, I. Abalakin, A. Minkin, "Unstructured mesh processing in parallel CFD project GIMM," in Parallel Computational Fluid Dynamics 2005, Elsevier B. V., 2006, pp. 501508. ]]
7. Беклемышева К. А., Петров И. Б., Фаворская А. В.
Численное моделирование процессов в твердых деформируемых средах при наличии динамических контактов с помощью сеточно-характеристического метода // Математическое моделирование. 2013. Т. 25, № 11. С. 3-16. [[ K. A. Beklemysheva, I. B. Petrov, A. V. Favorskaya, "Numerical simulation of processes in solid deformable media in the presence of dynamic contacts using the grid-characteristic method," in Mathematical Models and Computer Simulations, vol. 6, Issue 3, pp. 294-304, May 2014. ]]
8. Якобовский М. В. Инкрементный алгоритм декомпозиции графов // Вестник Нижегородского университета им. Н. И. Лобачевского. Серия: Математическое моделирование и оптимальное управление. 2005. № 1. С. 243. [[ M. V. Iakobovski, "Incremental graph decomposition algorithm," (in Russian), in Vestnik of Nizhny Novgorod University, no. 1. p. 243, 2005. ]]
9. Сухинов А. И., Чистяков А. Е., Хачунц Д. С. Математическое моделирование движения многокомпонентной воздушной среды и транспорта загрязняющих веществ // Известия ЮФУ. Технические науки. 2011. № 8 (121).
C. 73-79. [[ A. I. Sukhinov, A. E. Chistyakov, and
D. S. Khachunts, "Mathematic modelling of multicomponent air motion and pollutants transportation," (in Russian), in Izvestiya SFedU. Engineering Sciences, no 8. pp. 73-79, 2011. ]]
10. Четверушкин Б. Н. Пределы детализации и формулировка моделей уравнений сплошных сред // Матем. моделирование. 2012. Т. 24, № 11. С. 33-52. [[ B. N. Chetverushkin "Resolution limits of continuous media mode and their mathematical formulations," in Mathematical Models
and Computer Simulations, vol. 5, Issue 3, pp. 266-279, May 2013. ]]
11. Белоцерковский О. М., Гущин В. А., Щенни-
ков В. В. Метод расщепления в применении к решению задач динамики вязкой несжимаемой жидкости // Журнал вычислительной математики и математической физики. 1975. Т. 15, № 1. С. 197. [[ O. M. Belotserkovskii, V. A. Gushchin, V. V. Shchennikov, "Use of the splitting method to solve problems of the dynamics of a viscous incompressible fluid," in USSR Computational Mathematics and Mathematical Physics, 15:1, 190-200, 1975. ]]
12. Гущин В. А., Матюшин П. В. Математическое моделирование и визуализация трансформации вихревой структуры течения около сферы при увеличении степени стратификации жидкости // Журнал вычислительной математики и математической физики. 2011. Т. 51. № 2.
C. 268-281. [[ V. A. Gushchin, P. V. Matyushin, "Numerical simulation and visualization of vortical structure transformation in the flow past a sphere at an increasing degree of stratification," in Computational Mathematics and Mathematical Physics, 51:2, 251-263, 2011. ]]
13. Сухинов А. И., Чистяков А. Е., Шишеня А. В. Оценка погрешности решения уравнения диффузии на основе схем с весами // Математическое моделирование. 2013. Т. 25. № 11. С. 53-64. [[ A. I. Sukhinov, A. E. Chistyakov, A. V. Shishenya, "Error estimate for diffusion equations solved by schemes with weights," in Mathematical Models and Computer Simulations, vol. 6, Issue 3, pp. 324-331, May 2014. ]]
14. Сухинов А. И., Чистяков А. Е., Фоменко Н. А. Методика построения разностных схем для задачи диффузии-конвекции-реакции, учитывающих степень заполненности контрольных ячеек // Известия Южного федерального университета. Технические науки. 2013. № 4. С. 87-98. [[ A. I. Sukhinov, A. E. Chistyakov, and N. A. Fomenko, "Method of construction difference scheme for problems of diffu-sion-convectionreaction, takes into the degree filling of the control volume," (in Russian), in Izvestiya SFedU. Engineering Sciences, no. 4, pp. 87-98, 2013. ]]
15. Сухинов А. И., Чистяков А. Е., Тимофеева Е. Ф., Шишеня А. В. Математическая модель расчета прибрежных волновых процессов //Математическое моделирование. 2012. Т. 24, № 8. С. 32-44. [[ A. I. Sukhinov, A. E. Chistyakov, E. F. Timofeeva, A. V. Shishenya, "Mathematical model for calculating coastal wave processes," in Mathematical Models and Computer Simulations, vol. 5, Issue 2, pp. 122-129, April 2013. ]]
16. Сухинов А. И., Хачунц Д. С. Программная реализация двумерной задачи движения воздушной среды // Известия Южного федерального университета. Технические науки. 2013. № 4. С. 81-86. [[ A. I. Sukhinov and
D. S. Khachunts, "Multicomponent air motion task including condensation and evaporation," (in Russian), in Izvestiya SFedU. Engineering Sciences, no. 4, pp. 81-86, 2013. ]]
17. Чистяков А. Е., Хачунц Д. С. Задача движения многокомпонентной воздушной среды с учетом парообразования и конденсации // Известия Южного федерального университета. Технические науки. 2013. № 4. С. 15-21. [[ A. E. Chistyakov and D. S. Khachunts, "Program realization of air motion two-dimensional task," (in Russian), in Izvestiya SFedU. Engineering Sciences, no. 4, pp. 15-21, 2013. ]]
18. Самарский А. А. Теория разностных схем. М. Наука, 1989. [[ A. A. Samarskii, Difference Scheme Theory, (in Russian). Moscow: Nauka, 1989. ]]
19. Самарский А. А. Николаев Е. С. Методы решения сеточных уравнений. М.: Наука, 1978. 588 с. [[ A. A. Samar-skii and E. S. Nikolaev, Numerical Methods for Grid Equations, (in Russian). Moscow: Nauka, 1978; Basel: Bikhauser Verlag, 1989. ]]
ОБ АВТОРАХ
СУХИНОВ Александр Иванович, декан фак-та физики, математики, информатики. Дипл. по специальности «Электронные вычислительные машины» (ТРТИ, 1977). Д-р физ.-мат. наук (ИММ РАН, 1996). Иссл. в обл. построения и исследования математических моделей шельфовых систем и мелководных водоемов.
ХАЧУНЦ Дианна Самвеловна, инж.-иссл. НИИ многопроц. выч. систем им. акад. А. В. Каляева. Дипл. Учитель математики и информатики (ТГПИ, 2009). Канд. физ.-мат. наук (ЮФУ, 2014). Иссл. в обл. вычислительной математики ЧИСТЯКОВ Александр Евгеньевич, доц. каф. интел. и многопроц. выч. систем. Дипл. Математик, системный про-грммист (ТРТУ, 2005). Канд. физ.-мат. наук (ЮФУ, 2010). Иссл. в обл. вычислительной математики
METADATA
Title: Mathematical model of impurities in the atmospheric boundary layer and its program implementation on a multiprocessor computer system. Authors: A. I. Sukhinov1, D. S. Khachunts2, A. E. Chistyakov3 Affiliation:
1 Taganrog Institute of AP Chekhov (branch) VPO "Rostov State University of Economics" (RINH), Russia.
2,3 FGAOU VO "Southern Federal University" (SFU), Russia. Email: 1 [email protected], 2 [email protected],
3 [email protected]. Language: Russian.
Source: Vestnik UGATU (scientific journal of Ufa State Aviation Technical University), vol. 19, no. 1 (67), pp. 213-223, 2015. ISSN 2225-2789 (Online), ISSN 1992-6502 (Print). Abstract: The aim is to develop a mathematical model of the motion of a multicomponent air environment, taking into account the transport of contaminants (pollutants) and heat, phase transitions, as well as the influence of vegetation (forest plantations) on the distribution of pollutants in the atmosphere; construction and study of difference schemes approximating the original problem; software implementation of the developed parallel algorithms and conducting numerical experiments modeling the motion of a multicomponent air, including in relation to recreational environment of coastal systems. Since the problem of forecasting the spread of pollutants should be dealt with in real or accelerated time scales, on grids including 106-109 nodes need to develop parallel algorithms hydrodynamics problems on massively parallel systems. A new mathematical model of aerodynamic processes, taking into account the high humidity, the variability of atmospheric pressure and temperature, etc., Which is especially characteristic of coastal areas. Developed parallel algorithms for the study of these models, implemented as a set of programs. Key words: air; two-dimensional discrete model; pollutants. About authors:
SUKHINOV, Alexander Ivanovich, Dean of the Faculty of Physics, Mathematics, Computer Science. Dipl. in "Electronic Engineering" (TRTI, 1977). Dr. Sci. Sciences (IMM RAS, 1996).
KHACHUNTS, Dianna Samvelovna, Research Engineer Research Institute of multiprocessor computer systems Academician AV Kalyaeva "Southern Federal University". Dipl. Teacher of Mathematics and Computer Science (TSPI, 2009). Cand. Sci. Sciences (SFU, 2014). Inst. in the region. computational mathematics.
CHISTYAKOV, Alexander Evgenevich, Associate Professor, Dept.. "Intelligent and multiprocessor computer systems" (IMS) "Southern Federal University". Dipl. Mathematician, system progrmmist (TSURE, 2005). Cand. Sci. Sciences (SFU, 2010).
