CC BY
84
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Говорухин В.H., Моргулис А.Б.,Тютюнов Ю.В. Медленный таксис в модели хищник-жертва // Докл. РАН. - 2000. - С. 372, N 6. - С.730 - 732.
2. Латун B.C. Устойчивость системы фитопланктон-зоопланктон-рыба // Экологи-
ческая безопасность прибрежной и шельфовой зон и комплексное использование ресурсов шельфа. - Севастополь, 2004. - Вып. 10. - С. 211 - 218.
3. . .
фитопланктон-зоопланктон-рыба // Морской экологический журнал. 2005. N 4. T. IV. 2005. С.49-60.
Никитина Алла Валерьевна
Технологический институт федерального государственного образовательного учреждения высшего профессионального образования «Южный федеральный университет» в г. Таганроге;
E-mail: [email protected]
347928, Россия, Таганрог, ГСП 17А, пер. Некрасовский, 44 Тел.: 8(8634) 37-16-06
Nikitina Alla Valerievna
Taganrog Institute of Technology - Federal State-Owned Educational Establishment of Higher Vocational Education “Southern Federal University”
E-mail: [email protected]
44, Nekrasovsky, Taganrog, 347928, Russia, Ph.: +7(8634)37-16-06
519.86
А.И. Сухинов, B.K. Гадельшин, Д.С. Любомищенко
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ РАСПРОСТРАНЕНИЯ ВРЕДНЫХ ПРИМЕСЕЙ В АТМОСФЕРЕ ГОРОДА НА ОСНОВЕ
IOCV МЕТОДА2
В данной работе рассмотрен подход моделирования загрязнения окружающей среды от автотранспортных средств на основе уравнений, записанных в универсальной форме, позволяющий существенно упростить процесс компью-.
метода контрольного объема. Вычислительный эксперимент выполняется на кластере распределенных вычислений.
SIMPLE - метод; моделирование распространения загрязнения; параллель.
A.I. Sukhinov, V.K. Gadelshin, D.S. Lyubomishchenko
MATHEMATICALAL MODELING OF DISTRIBUTION OF POLLUTION IN CITY ATMOSPHERE BASED ON IOCV-METHOD
2 Работа выполнена по гранту РНПВШ 2.1.1/6584
In the paper an approach for pollution modeling from motor transport is developed. It based on equations written at universal form what allows to simplify the modeling process. The numerical realization is based on control volume method and performed on distributed memory cluster.
SIMPLE-method; pollution modeling; parallel programming.
Задача оценки загрязнения воздушной среды города становится все более актуальной из-за постоянного роста количества загрязнителей. Автомобильный транспорт и промышленные предприятия вносят наиболее существенный вклад в неблагоприятную картину загрязнения современного города с развитой промышленной и автодорожной инфраструктурой. В связи с этим возникает необходимость оперативной оценки загрязнения атмосферы города для контроля и выработки рекомендаций по нормализации уровня выбросов. Математическое моделирование является одним из наиболее приемлемых способов решения поставленных
.
В данной работе предложена математическая модель расчета уровня загряз. - -нии конвекции-диффузии, записанных в форме уравнений переноса универсально. -пользования универсальных численных алгоритмов для всех уравнений системы.
Запишем базовую трехмерную систему уравнений в форме универсального уравнения переноса для переменной Ф [1, 2, 6]:
д(р„,Ф)
dt
ф)
дх
_д_
дх
eff
дФ
дх
(1)
+
Значения коэффициентов для разных уравнений приведены в табл. 1:
Таблица 1
Значения коэффициентов для различных типов уравнений
Ф г eff № формулы
1 0 0 (2)
и 1 Veff -дР/ дхі + Su (3)
T keff/ / Pm Ql cpm (4)
В формулах (1) - (4) р - плотность среды, - эффективная вязкость, р
- давление, - компоненты вектора скорости в направлении трех координатных
осей, и^ - компоненты вектора скорости в направлении трех координатных осей,
удовлетворяющие дискретному аналогу уравнения неразрывности (2).
Далее для построения численного [6] метода будем рассматривать универсальное уравнение (1). Данное уравнение удобно представить в виде (5): д(рФ) дд дд2 дд
9t
- + —і + —^ + дх1 дх2 дх
3 = S,
3
где qj совместный конвективно-диффузионный поток в направлении j :
дФ ■
qj =pu^-reff’j =1,3.
(5)
(5) , -
вании по контрольному объему (ЮСУ).
(рФ) -(рф) AxAyAz+Г() -(J+
J ) -(J
y/n \ у
А+[(Jz), -(Jz )b] А‘-
(7)
=SuAtAxAyAz+SpOPAtAxAyAz
В левой части уравнения (7) J представляют собой потоки через соответствующие грани контрольного объема. Для простоты опустим индексы координатных направлений у J. Определим J в точках w , e, s, n, b и t:
Jw = Dw [[w - АФР ], Je = De [Op - AeO, ], (8)
J = Ds [фS - AOp ] Jn = Dn [Op - AnO„ ]
Jb = Db [Фв - AbOp ], Jt = Dt [Фp - AtOr ].
На рис. 1 представлено расположение узлов и граней соседних контрольных .
Коэффициенты A и B в (8) определяют специфику потока через грани контрольного объема и определяются следующим образом:
A (р)=A (И) + max (- P,0) (9)
B (P ) = A (|P|) + max (P,0)
В формуле (9) P - число Пекле, A (PI) - , -
вий движения среды в рассматриваемой задаче. Число Пекле определяется по
:
P = F/D,
F - , D -
.
Значения F и D для точек w , e , s, n, b и t представлены формулами:
= [Р11),,АУ^' Fe = (Ри)е A-vAz'
= AVAZ' Я = T/sxe 4vAz •
Fs = (p)s AxAz, Fn = (pu)n AxAz , Ds = %s AxAZ , Dn = ^ AxAZ , (10)
Fb = (Pu )b AxAy,
Ft = (pu )t AxAy,
Db =Yb/sxu AxAy,
D, =
AxAy.
Рис. 1. Схема расположения узлов и граней соседних контрольных объемов
Функция А (|Р|) вводится в
схему для корректной аппроксимации конвективных членов исходного уравнения. В табл. 2 представле-
ны возможные варианты для А (Р|).
2
Способы задания А (р|)
Схема аппроксимации Значения А (р|)
Центрально-р^ностная 1 - 0.5 р|
Противопоточная 1
Степенная тах (0, (1 - 0.5 |р| )5)
Ф£
После подстановки в (9) значений (10) получим выражение вида ррДхАуДг_ + аЕ + ^ + а^ + ав + БрДхДуДг + Ре - Р„ + Рп - ^ + Р,- Ръ
Аґ
— Фшаш + ФЕаЕ + Фхах + Фмам + Фвав + Фтат + SuЛАyЛ + Ф^
Рр АхАуА Аґ
■ (11)
(2)
центром в точке Р:
Рр -Рр Аґ
(12)
Теперь вычтем из (11) (12), умноженное на Фр. После приведения подобных и введения новых обозначений получим:
арФр = ажФж + аЕФЕ + а5Ф5 + анФн + авФв + атФг + Ъ , (13)
а^ А
+ тах ,0), аЕ = А (|р |) + тах (-Рв ,0):
+ тах (,0), аы = А (|рп|) + тах (-Гп,0),
ав = А (|рь |) + тах (ь,0) ат = А (|р |) + тах (-р ,0)
а5 — А
а
о
, Рр АхАуА.
Аґ
а
аш + аЕ + а5 + а^ + ав + ат + а0Р - ЯрАхАуАг .
р ' ^Е ' ^ ^ ' “В ' ^Т ' ^Р
Ь — БиАхАуА + а0р ФР.
Или в другой форме:
[ар + 5р] Фр = ашФш + аЕФв + аБФх + аыФ^ + авФв + атФт + Su ,
где
ар — аф + аЕ + а8 + а^ + ав + ат
Яр — Рр АхАуАг + БрАхАуАг, Аґ
Яи — ЯиАхАуАі + Ф0
рр АхАуЛг Аґ
(14)
(15)
(16)
(17)
(18)
(19)
(20) (21)
В выражениях (20) и (21) члены Su и Sp, стоящие в правой части, сформированы в ходе линеаризации источникового члена исходного уравнения, a Su и
Sp . -
ют большую роль при постановке граничных условий, а также в задании итераци-
(20).
На рис. 2 представлена модель области проведения численного эксперимента. На нем обозначены основные грани, на которых требуется задание граничных условий. Грани front и back представляют соответственно области втекания и истечения атмосферного воздуха, а остальные внешние грани считаются непроницаемыми. Грани поверхностей зданий (l, r, f, bc, b, t) являются также непроницаемы.
Опре делим математические формулировки граничных условий [2].
На грани front: u = uin, v = 0, w = 0 ;
, , d 2u
на грани back: ---------= 0, v ■
dx2
0,
w = 0:
du
на гранях left, right, top: — = 0 , dn
^=0, dn
dw
dn
= 0, где n - вектор внеш-
Puc. 2. Схема области моделирования
ней нормали к грани;
на грани bottom: — = -uu , dn ^
-dv = -juv, — = 0, где и - коэффици-dn дп
ент трения о подстилающую поверх;
Ф = u, v, w, где u , v и w - компоненты вектора скорости в направле-
ниях Ox, Oy и Oz соответственно.
, ,
высокой вязкостью. Счет в них осуществляется так же, как в областях с нормальной .
Учет граничных условий в численной модели происходит с помощью коэффициентов Sp и Su. Например, в уравнении
(18) при Ф = u полу чаем:
на грани front: Su = Su + aWun,
Sp = Sp + aW , aW = 0, на грани back: aE = 0.
Задание граничных условий на гранях left, right, top, bottom с учетом трения рассмотрим более подробно. На рис. 3 пред-
Рис. 3. Профиль вектора скорости вблизи стены
ставлена схема поведения компоненты скорости и вблизи стены
Предполагается, что движение среды вблизи стены ламинарное. Тогда коэффициент напряжения имеет вид:
Тогда сила трения выражается в виде:
ау.
р. = —тА = -ц—^л , * " АУр
(19)
(20)
где Ц - коэффициент трения о стену, ир - значение скорости в пристеночном .
(19) (20)
скорости от дистанции до стены.
Тогда в разностной схеме учет граничного условия осуществляется с помощью введения фиктивного источника:
Бр = -
Ц
АУр
л .
Согласно 81МРЬБ-методу, в задаче определения поля скоростей требуется модифицировать универсальную схему.
На рис. 4 приведен пример контрольного объема для компоненты скорости и . Запишем универсальное уравнение (13) для Ф = и в точке w:
= X ал +ъ+((г- Рр )А, (21)
где апЬ и 11пЬ - значения коэффициентов в соседних узлах, Ъ -правая часть соответствующего уравнения для Ф = и, рш и рр -
значение давления в точках Р и Ам - площадь поверхности, на которую действует перепад давления (Ам = ЛуЛ).
Запишем (21) для приближен*
ного поля давления р и компо-
Рис. 4. Пример контрольного объема для компоненты скорости и
ненты скорости ич
(22)
аши1 = X апъи1ъ + Ъ +( Р*Ш - РР ) А .
Предположим, что истинное поле давления связано с приближенным соотношением р = р* + р'. Необходимо выяснить, как будут изменяться соответствующие компоненты скорости:
и = и* + и , V = V* + V , м = мм + М . (23)
Вычтем из (24) (25), тогда при и получим соотношение:
аши'м = X апЪи'пЪ + Ъ + {р'ш - РР )А . (24)
= <
Член ^ апЬыпЬ согласно методу SIMPLE можно отбросить. В итоге получим выражение для поправки скорости u :
u'w = dw (p'w - p'p ) (25)
где d = Aw/ .
w / aw
Тогда поправочная формула для скорости может быть переписана в виде:
uw = u*w + dw p'w - p'p). (26)
Подставляя соответствующие компоненты скорости в (12), можно получить
уравнение для отыскания поправки к давлению:
aPpP = awp'w + aEp'E + aSpS + aNp'N + aBpB + aTpT + Ь, (27)
где
aw = Pwdw^y^z, aE = pedeAyAz , aS = psdsAxAz , aN = pndnAxAz, (28)
aB = pbdb AxAy, aT = ptdtAxAy , aP = aw + aE + aS + aN + aB + aT,
Ь =(Pp PpAtAxAyA + ^pu*)w )J AyA + [pv)s _(pv*)n]AzAx + . (29) +\[pw* )Ь ~{pwt )J AxAy
P,
контрольных объемов необходимо прибегнуть к соответствующей интерполяции.
Г раничные условия для уравнения (27) формулируются в виде p' = 0.
Задача распространения вредной примеси формулируется на базе универсальной формы записи уравнений переноса при Ф = (р.
[aP + Sp]pP = awpw + aEpE + aSpS + aNpN + aBpB + aTpT + Su . (30)
p образует контрольный объем с центром в точке P.
К уравнению присоединяется начальное условие р = р0 (фоновое загрязнение перед началом моделирования).
:
др = 0 - условие беспрепятственного ухода примеси за пределы области мо-дп
делирования, др = ар — условие поглощения подстилающей поверхностью при-dz
.
Sp Su -
личия источников загрязнения в исследуемой области [4, 5].
[3] -
.
кластере распределенных вычислений. На рис. 5 представлена схема деления области на подобласти.
В каждой подобласти соответствующим узлом кластера реализуется расчет необходимого дискретного уравнения. Между границами подобластей производятся обмены приграничными элементами из-за специфики шаблона разностной схемы. Функции, реализующие обмен данными между различными компьютер, MPI.
Вычислительный эксперимент был проведен для разного размера сеток и количества задействованных процессоров. В табл. 3 приведены результаты эффективности алгоритма.
Численный эксперимент проводился для участка уличнодорожной сети города Таганрога с захватом близлежащих кварталов. На рис. 6 представлены результаты моделирования при следующих условиях: интенсивность подвижных источников (автотранспортного потока) 1200 единиц/час, загрязняющая примесь - однокомпонентный газ СО, скорость движения 30 км/ч, мощность выбросов 5,5
/( / ), 5 / ,
температура 150 С, давление 750 мм. р. ст. Номерам изолиний соответствуют следующие кон-1 .:
1 - 6,2 мг/м3, 2 - 5,1 мг/м3,
Рис. 5. Схема деления области на подобласти
3 - 2,1 мг/м3, 4 - 0,7 мг/м3.
Рис. 6. Результаты моделирования
В работе предложен алгоритм решения задачи распространения загрязняющих веществ от автотранспортных средств в условиях реальной городской за.
Таблица 3
Эффективность алгоритма_________________________
p/N 151x646x20 800x600x20
2 0,72 0,851
4 0,73 0,823
6 0,712 0,8
8 0,67 0,74
Численная реализация основана на использовании универсальной дискретной схемы, которая позволяет в процессе вычисления подбирать оптимальные аппроксимации конвективным членам для усиления устойчивости алгоритма. Данный подход позволяет упростить процесс моделирования и увеличить точность за счет использования информации о дополнительных факторах, поведение которых описывается с помощью универсального уравнения переноса. Построен параллельный алгоритм и представлены результаты моделирования для одного из перекрестков улично-дорожной сети города Таганрога.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Берлянд М. Е. Прогноз и регулирование загрязнений атмосферы. - Л.: Гидроме-тиоиздат, 1985. - 271 с.
2. Марчук Г. И. Математическое моделирование в проблеме окружающей среды.
- М.: Наука, 1982. - 319 с.
3. . . . - .: , 1999.
4. Методика расчетов выбросов в атмосферу загрязняющих веществ автотранспортом на городских магистралях. - М.: НИИАТ, 1997. - 54 с.
5. МатвеевJJ.T. Физика атмосферы. - СПб.: Гидрометеоиздат, 2000. - 779 с.
6. Патанкар С.В. Численные методы решения задач обмена и динамики жидкости. -М.: Энергоатомиздат, 1984. - 152 с.
Сухинов Александр Иванович
Технологический институт федерального государственного образовательного учреждения высшего профессионального образования «Южный федеральный университет» в г. Таганроге;
E-mail: [email protected]
347928, Россия, Таганрог, ГСП 17А, пер. Некрасовский, 44 Тел.: 8(8634)371606
Гадельшин Валерий Камельянович E-mail: [email protected]
Любомищенко Денис Сергеевич E-mail: [email protected]
Suchinov Alexsander Ivanovich
Taganrog Institute of Technology - Federal State-Owned Educational Establishment of Higher Vocational Education “Southern Federal University”
E-mail: [email protected]
44, Nekrasovsky, Taganrog, 347928, Russia, Ph.: +7(8634)37-16-06
Gadelshin Valeriy Kameljanovich E-mail: [email protected]
Lubomishenko Denis Sergeevich E-mail: [email protected]
519.8:533
Б. В. Сидоренко
MRT LATTICE BOLTZMANN МЕТОД В МОДЕЛИРОВАНИИ ГИДРОДИНАМИКИ МЕЛКОВОДНЫХ ВОДОЕМОВ
В данной работе рассматриваются MRT Lattice Boltzmann модели в вычис-.
для реальных водоемов. При численном моделировании была задействована D3Q19 модель, которая была модифицирована для некубических сеток, то есть с преобладанием какого-то шага по пространств}’. MRT модели показали высокую скорость вычислений и хорошую устойчивость к входным параметрам.
Решеточное уравнение Больцмана; Решеточный метод Больцмана; одиночная релаксационная модель; многовременная релаксационная модель; Граничные .
B.V. Sidorenko
MRT LATTICE BOLTZMANN METHOD IN FLUID DYNAMICS FOR SHALLOW WATER BASINS
In the work MRT Lattice Boltzmann model in CFD are considered. With their help numerical experiments for real basins are made. At numerical modeling used D3Q19 model which has been modified for not cubic grids, that is with prevalence of one of step in space. MRT models have shown: high speed of calculations and stability against entrance parameters.
LBE; LBM; LB-BGK; Lattice Boltzmann; MRT; Bounce Back BC.
Введение
Lattice Boltzmann метод (LBM), использующий минимальные дискретные кинетические модели для решения задач в механике жидкости и других областях физики, привлекает большое внимание в последние годы [1-4]. Вместо прямого
- LBM
уравнения Больцмана (LBE), которое описывает развитие распределения ансамбля частиц на решетке, коллективное поведение которых асимптотически представля-
. LBE,
получаются течения жидкости, представляемые слабо сжимаемыми уравнениями -.
LBE -
бой релаксационные модели. Одна из наиболее общих - это одновременная релак-186
