Перенос и рассеяние воздушным потоком тяжелого газа, эмитированного точечным источником

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

В настоящей работе рассматривается задача о движении автомобиля, моделируемого как точечный источник, эмитирующий поток газа с интенсивностью m и мощностью £ , содержащий загрязняющую воздух пассивную газовую примесь ( q = 0,1 m , то есть 10 % эмитируемой массы). Решение такой задачи обусловлено необходимостью исследования переноса и рассеяния отработанных газов от потока транспортных средств, моделируемого как множество отдельных подвижных точечных источников загрязняющей примеси различной интенсивности в произвольной пространственной области со сложной геометрией.

Использование численных методов в настоящее время является основным способом построения решений краевых задач газовой динамики. При исследовании подобных задач большее внимание уделяется конечно-разностным методам [1-4], требующим относительно малых вычислительных ресурсов. Значительный интерес вызывают схемы факторизации эволюционных дифференциальных уравнений [5, 6], по- зволяющие «расщеплять» многомерные задачи на последовательности одномерных задач, что приводит к существенному повышению эффективности вычислительных алгоритмов. В ряде исследований [7-13] для решения задач движения сплошной среды используются методы конечных и граничных элементов с различными видами аппроксимации полей скорости, перемещения, плотности, давления, температуры и других характеристик.

Для решения прикладных задач механики сплошных сред интенсивно развивается и успешно применяется метод крупных частиц (метод Давыдова), позволяющий выполнять расчеты вихревых структур с учетом отрывных явлений, исследовать фильтрационные и струйные [14] течения, газо- и гидродинамические потоки с большими переме щениями и соударяющимися поверхностями раздела, движение многокомпонентных [15, 16], сыпучих и деформируемых сред [17], течения сквозь проницаемые объекты [18] и многие другие процессы. С использованием системы уравнений Эйлера построены трехмерные вычислительные модели взаимодействия струй с поперечными потоками [19, 20], истечения газа из отверстий и каналов [21, 22], воздействия ударных волн на препятствия и границы раздела сред [23, 24, 25], влияния локализованных источников энергии и массы на газодинамические характеристики газовых потоков [26, 27].

Для моделирования движения газовоздушной смеси используется система дифференциальных уравнений Эйлера (в дивергентной форме [17, 18, 29]), описывающая движение сплошной среды в рассматриваемой области, включающая эволюционные уравнения неразрывности

°P" + V(pV ) = q 5( rs),(1)

О t движения

°^ + V (pVV ) + V P = pg + q V q 5 (rs),(2)

О t полной удельной энергии

°(рУ)+ V(p UV ) + V( P V ) = 0,(3)

состояния (адиабатический процесс)

p( k -1)( U - V • V/ 2)-P = 0(4)

и концентрации (в рамках дифференциального приближения [28])

— + V ( C V ) = V ( XV C ) ₊ q CO 2 s ( r S ) .

Здесь р - плотность газовоздушной смеси; V - вектор скорости потока с компонентами V x , V y , V z ; U - полная удельная энергия; P - давление; С - концентрация загрязняющей примеси; к - показатель адиабаты; 5 () - дельта-функция Дирака; r _s - радиус-вектор положения точечного источника; X - коэффициент диффузии; q - интенсивность выброса смеси отработанных газов; q _CO 2 - интенсивность выброса диоксида углерода; V_q - скорость выброса смеси отработанных газов из выхлопной трубы.

Точечный источник с координатами x_S , y_S , z_S моделирует положение отверстия выхлопной трубы автомобиля в расчетной области G (рис. 1).

Рис. 1. Схема расчетной области и обозначения ее границ

В начальный момент времени в области G известны распределения компонент вектора скорости, плотности, энергии, концентрации и давления:

V x ( °, x , У , z ) = V о , V y ( °, x , У , z ) = ° , V z ( 0, x , y , z ) = 0 , C ( 0, x , y , z ) = 0 , p ( 0, x , y , z ) = р 0 , U ( 0, x , y , z ) = U 0 , P ( 0, x , y , z ) = P 0 x , y , z e G .

Граничные условия имеют следующий вид:

v x = о, v y = о, v z = o, ap/a z = о, a и/ a z = о, a с / a z = o, x , y , z ea G 1 ;

a vx/az = о, a vy/az = о, a vz azz = o, ap/az = о, а и/ az = о, a Cl az = о , x, y, z e aG 2;                      (6)

v x = v» a v y /a x = о, a v z la x = о, p = p о , и = и о , с = о, x , y , z ea g 3 ;

a vx axx = о, a vy /ax = о, a vz axx = о, ap/ax = о, a и axx = о, a с/ az = о, x, y, z e a g 4;

a v x a y = о, v y = о, a v z a y = о, ap/a y = о, a и a yy = о, a с] a z = о, x , y , z e a g ₅;

a v x a yy = о, a v y /a y = о, a v z /a y = о, ap/a y = о, a и a yy = о, a cl a z = о, x , y , z e a g ₆.

Для учета плавучести загрязняющей газовой примеси используется приближение Буссинеска [о], согласно которому в предположении о малости концентрации C используется представление плотности газовоздушной смеси в виде

p ( с ) = p ( с о ) + p' c ( с о )( C - с о ) = P o + a C ,            (7)

где p ₀ =p ( C о ) ; a = p' _C ( C _о ) - концентрационный коэффициент плотности; начальная концентрация загрязняющей примеси отсутствует, то есть C о = о.

Пусть в некотором объеме Q содержатся воздух, имеющий плотность p о, и пассивная газовая примесь с плотностью pп, занимающая объем Оп. Массы воздуха и примеси соответственно равны mо =p0 (Q-Qп ) и mп =pпQп. Плотность газовоздушной смеси в рассматриваемом объеме Q определяется выражением m о + m п _p о (^-^ п l.p^ п_п   p о p п^ п p ==1= p

Q       Q Q ^о p _п Q

= P o

- ^p ° C + C = P o +f 1 -^ ^p п               V ^p п>

C ,

где С - концентрация примеси в объеме Q. Сравнение полученного выражения и соотношения (7) позволяет определить величину концентрационного коэффициента плотности:

a = 1 --^P 0 .

^p n

Для корректного учета скоростных граничных условий целесообразно полное давление газо-воздушной смеси Р в уравнениях (2)-(4) представить в виде суммы гидростатического давления Р _г и отклонения p давления от гидростатического,

P = Рг + Р , причем Рг = Р0 - рgz, Р0 - давление воздуха на уровне поверхности земли; g - ускорение свободного падения.

С учетом приближения выражений (7) и (8) уравнения движения (2), полной удельной энергии (3) и состояния (4) принимают вид

^^PV + V(pVV ) + Vp = a Cg + q V q 5(rs ),

^^P^^ + V(p UY ) + V( p V ) = -V( Рг V),(10)

p(k-1)(U- V • V/2)-p = Рг.(11)

В соответствии с идеей метода крупных частиц система уравнений (1), (5), (9)-(11) расщепляется по физическим процессам. Расчет каждого временного слоя разбивается на три этапа. Для реализации первого (эйлерова) этапа считается, что изменяются лишь величины, относящиеся к ячейке в целом, а сплошная среда предполагается заторможенной. Система уравнений (1), (3), (5) и (11) представляется в виде

'p = const,

C = const,

а , .

—(p V ) ^+V p = ^a Cg , а t

- ( p U ) + V ( p V ) = -V ( Р г V ).

Уравнения этой системы записываются в виде явных конечноразностных схем:

V xijk ^- V xijk   ( pi + 1/2 jk    p i - 1/2 jk ) ^A t /p ijk h x

V yijk ^- V yijk  ( P ij + 1/2 k   p j - 1/2 k ) ^A t /^pijk h y ,

V zijk ^- V zijk   ( P ijk. + 1/2    p ijk - 1/2 ) ^A t /^pijk h z   ^U C ijk g ^A t /^p ijk ,

U i^jk U ijk ^A t | [ ( ^pV x ⁾ i + 1/2 jk ^{( P}V ^{x )} i - 1/2 jk J / h ^{x + ( P}V y ⁾ ij _{+ 1/2 k} ^{( P}V y ⁾ ij _{- 1/2 k} / h y ⁺

⁺ ( ^pVz ⁾ ijk + 1/2 ^{- ( P} V ^z j - 1/2 / h ^{z +} _⁽ P T V ^{х )} i + 1/2 jk ^{- (} P T V ^{х )} i - 1/2 jk ^hx ⁺

+ ( P_rv J -( P_rv J /h _v + Г ( P_rV _z)..,„-( P_rV_z L _1ZJ/ h_z Pa_k .

_V ^{1 y} ‘ ij +1/2 k V ¹ y ‘ ij -1/2 k J/ ^y L^{V 1} zj ^+1/2 ' ¹ z ^^lj^{k -1/2} J/ z^^ijk

В приведенных соотношениях величины с дробными индексами определяются как

^Л i ± 1/2 jk

^- ( ^Л ijk ^{+ Л} i ± 1 jk ) /² , ^Л у ± 1/2 k ^- ( ^Л ijk ^{+ Л} д ± 1 k ) /2 ,

^Л ijk ± 1/2

^{- (} ^ ijk ⁺ ^ ijk ± 1 ⁾ /² ,

где Л принимает значения p , pV_x , PV y , pV_zP_t V_x , P_TV_y , P_TV_z .

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

^w                        ^^

(Л V )

V ^x)i + 1/2 jk

^Л ijk V xi + 1/2 jk ,    V xi + 1/2 jk - 0,

^w                     ^w

^Л i + 1 jk ^Vxi + 1/2 jk ,    ^Vxi + 1/2 jk ^< °:

(л V J    =)

V     ^x)i - 1/2 jk

^w

^Л i - 1 jk V xi - 1/2 jk ,

.w

V xi - 1/2 jk ^- °,

(Л V y)

\    y )ij + 1/2 k

(Л V y)

\     y/ij - 1/2 k

^w ^Л ijk^Vxi - 1/2 jk ,

^w

^Л Ц1Уyij + 1/2 k ,

^w

^Л ij + 1 k V yij + 1/2 k ,

^w

^Л у - 1 k V yij - 1/2 k ,

.w

^Л ijk^Vyij - 1/2 k ,

^w

V xi - 1/2 jk <  °;

^w

V yij + 1/2 k ^- °,

^w

V yij + 1/2 k <  °;

^w

V yij - 1/2 k ^- °,

^w

V yij - 1/2 k <  °;

( Л V -	-. Л jk + 1/2	7 Л ijkVzijk + 1/2 ,	^w V zijk + 1/2 - 0,
		Л ijk + 1 V zijk + 1/2 ,	^w V zijk + 1/2 <  0;
( Л V -	=< jjk - 1/2	17 Л Ц - 1 k V zijk - 1/2 ,	7 V zijk - 1/2 - 0,
		^ Л ijkVzijk -1/2 ,	7 V zijk - 1/2 <  0 •

В этих выражениях Л принимает значения р , р V x , р V y , р V z , р U , *^           *^            *^

C , CV x , CV y , CV - .

На третьем, заключительном, этапе окончательные значения массы, импульса и энергии в момент времени t = t + A t определяются за-

конами сохранения массы, импульса и энергии, записанными с учетом промежуточных значений параметров потока:

¥^{+v (p V )-} о, а ( ^{p V} H^VV ) = 0 - £ ( р U И( р U V ) = ° CC + V ( C V ) = V(W C ) .

Разностная аппроксимация уравнений приводит к системе разрешающих соотношений

^р j ^- р j ^-A t <

(рVJ .  -(рVJ  , h + x/i+12jk   '    x/i-12jkJ/ x

⁺ [ ^{( р} V У ⁾ . . k -^{( р} V y ^:     . k J /h ■ '■’ V zz ⁾ jk . . ^-(р V z ⁾ jk - 1/2 J Ih-

V ... =v -.p-. o..,-M\ xjk     xijk ^р ijk| ^р ijk

С р VxVx) ,   -(р VxVx) ,    I Pithx + x x/i+12 jk        x x/i-12 jk J] v4jk x

^{+ (} P V x V y ⁾ j + 12 k ( p V x V y ⁾ у _ ij2 k ^ ^u^jkh/ ■ I p V x V z j + 12 ( p V x V z j - 12 J?ijk h z

A ~

V yijk = V yijk ^p ijk I ^p ik ^{- A} t *

^{( p} V y V x ⁾ i ₊ 12 jk ^{( p} V y V x ⁾ i - 12 jk _ /^p ijk h x ⁺

⁺ I P ^{V y V y )} _y + 12 k ^{(p V y V y )} „ - 12 k ^ \^{p i}j^kh^ ⁺ | p V y V z L + 12 ^(p V y V z j - 12    Zkjjh.

V zjjk = V zjjk p yk /p ijk ^-A t *

^(p V z V x ⁾ i + 12 jk ^(p V z V x ⁾ i - 12 jk J^PHk^ h x ⁺

^ ^                ^ ^                            ^ ^               ^ ^

⁺ ^ ( ^p V z V y ⁾ j + 12 k Ip V z V y ⁾ j - 12 k \l^pij^{k h} y ' J^p V z V z j + 12 ^-(p V z V z ⁾ jk - 12 У ^p j^k h ^z

A ~           / _

^U ijk = ^U ijk ^p У^к j ^p У^{к - A} t *

^{( p UV}x ⁾ _{z 4}12 jk  ^{( p UV}x ⁾ i - 12 jk \/ p ,kh x ⁺

■ p UU V ■ ) j + 12 k "I UU V ⁷. y ⁾ j - 12 k \/^p ik h y ■ ' p UU ^V'- ) jk ■       p UU ^V'z ⁾ jk - 12    ^A jk h z • .

C ijk = C ijk - A t | _ ( CV x ⁾ i + 11 jk ^-( CV x ⁾ i - 12jk \ Jk ⁺ ( CV y ⁾ i + 12 k ^- ( CV y ⁾ j - 12 k ^h ^{+ +} _ ^{( C} V ^{z )} jk + 12 ^{- ( C} V ^{z )} jk - 12 _ / h ^- ' ^{+ A} t * [⁽^^{V C )} 1 + 12^jk ' ^AV c ⁾ ^ - 1/2 kk - \/ h x ⁺

+ [ (№ C ) i + 12 jk (^{XV C} ) , - 12. Ik ]/ h y       ^{V C} ) i + 12 jk ( X ^{V C )} i - 12 jk \/ h z ^ .

В последнем выражении использованы обозначения (^ ^C )_M2k = ( ^Х i + 1 jk ^+Х jk )( C i + 1kk - C ijk )/ ² h x . ^(XV C ⁾ i - 12 jk ^{= ( X} jk ^+X i - 1 jk ⁾⁽ C ijk ^- C i - 1 jk ⁾ /2 ^hx ,

(^ C ) j . 1/2 k = (^X i . 1 jk ⁺^ k )( C i.4k ^- C ijk )/² h y • (^ C ) ₈ - V2k = (^X ijk ^{+ X} i-4k )( C ijk - C i-4k ) /² h y . ( ^XV C j + 1/2 = (^X i + 1 jk *4 ) ( C i . 1 jk - C ijk )/² h z . ( ^{XV C} ) jk - 1/2 = (^X ,k + ^X . - )( C j, - C , - 1,k )/² h z .

Поле давления вычисляется из уравнения состояния (11)

р.            /и _о 5(72 +jz2 +к2 + pijk (       )р ijk \ ijk ’ ) xijk yijk ' zijk

^{+ P} ijk

(у...

xijk

J2+(n xijk          yijk

J2+(^

yijk           zijk

^+ xzjk

^^^^H

Р г .

По окончании третьего этапа в ячейке, где расположен точечный источник, в соответствии с уравнениями (1), (5) и (10) пересчитывают ся значения плотности, энергии и концентрации с учетом заданных величин q, qCO2 и Vq:

V xijk ( ' ) = ^Vxijk + qV qx ^A t /p ijk h xhy h z , P yk ( t ) = P ijk + q ^A tlh x h y h z ,

C ijk ⁽ t ^{) =} C ijk ⁺ q CO₂ ^A t Ihx h y h z •

С использованием описанного алгоритма разработан комплекс программ для численного исследования пространственного переноса и рассеяния примесей, выбрасываемых точечным источником. Тестирование программного комплекса выполнено с использованием точного решения частной задачи об эмиссии сжимаемого газа точечным источником [31].

Для выполнения вычислительных экспериментов принято: интенсивность источника отработанных газов q = 1,2 г/с; интенсивность эмиссии диоксида углерода q _CO2 = 0,12 г/с; скорость V _q = 10 м/с; показатель адиабаты k = 1,67; плотность воздуха р ₀ = 1,29 кг/м ³ , плотность диоксида углерода (содержание в отработанных газах автотранспорта достигает 10 %, [0]) р _п = 1,98 кг/м ³ . Вычислительный эксперимент выполнялся с учетом ветра, изменяющегося по высоте рассматриваемой области линейно от 0 до 1 м/с (см. рис. 1).

^19 62^=                  15.79^ 375
а
^19 75^= _       =
"                           4.03
б
^19.85^^ ___ :    '---
3 76

в

^20.21 =-- _

4.40

г

Рис. 2. Распределение концентрации диоксида углерода (г/м³) в плоскости у = 0,5 м после 1,0 с ( а ), 1,25 с ( б ), 1,5 с ( в ) и 1,75 с ( г )

На рис. 2 показаны распределения концентрации загрязняющего вещества в плоскости у = 0,5 м после 1, 1,25, 1,5 и 1,75 секунд после начала рассматриваемого процесса. Приведенные на рисунке распределения концентраций показывают, что за счет силы тяжести относительно тяжелая примесь диоксида углерода опускается вниз, нарушая симметричное первоначальное распределение загрязняющего вещества в составе отработанных газов автомобильного транспорта.

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