Численный код для расчета аспирационных течений в промышленных цехах

DOI:

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

Имеющиеся в настоящее время инструменты основаны на простых инженерных расчетах, не способных адекватно описывать сложные динамические процессы и реальную геометрию помещений [1; 5]. Более точное решение описанной выше задачи представляется возможным при использовании численного газодинамического моделирования, которые точнее типовых инженерных расчетов и дешевле натурных экспериментов.

1.    Постановка задачи

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

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

р(у) = р о exp

/ ^-M 0^|gL j

I RT ^уу

р (у) = Р о exp

J - Md I s I 1

I KT ^уУ

где р ₀ и р ₀ — значения плотности и давления, соответственно, на уровне у = 0; М _о — молярная масса воздуха; g — ускорение свободного падения, вектор которого направлен противоположно оси у; R — универсальная газовая постоянная; T 0 — равновесная температура.

Динамику среды можно описать с помощью уравнений

U =

^ρ р и р ж ^е

/   ри   \ ри2 + р риж

\ (е + р ^)и /

/    рж   \ риж рж2 + р

\ (е + р)ж /

S =

0 р д \ ржд /

Здесь используются стандартные обозначения: р — плотность; V = (и, ж) — вектор скорости; р — давление; е — объемная плотность энергии; д — величина ускорения свободного падения (полагаем его отрицательным); t — время. Система уравнений (3) замыкается калорическим уравнением состояния

= р ₊ 2 V !

Y - 1 ⁺ 2

где γ — показатель адиабаты.

2.    Численная схема

Одним из наиболее распространенных методов решения системы (3) является применение идеологии MUSCL [4].

В расчетной области {жmin < ж < хтаж, ymin < у < утах} введем неподвижную эйлерову сетку жг  ^min +

(i - 2)^,

i — 1,2,... ,L, Аж =

^{( ж} та:с

-

^ж min ^{) L} ,

^у ,   ^y min +

(j - 2)Ау,

^j — 1,², . . . , ^M , ^{A y} — ^(у тах    ^y min ⁾ / M.

Предполагаем кусочно-линейное распределение переменных внутри ячеек этой сетки. Для интегрирования дискретного аналога системы (3) по времени будем использовать метод Хэнкока [2]. Это двухэтапный метод типа «предиктор — корректор». Предиктор состоит из двух шагов. Первый шаг: предполагая кусочно-линейное распределение параметров внутри расчетных ячеек, вычисляем предварительные значения на новом временном слое:

U*. = un ■ - — If fun • + PA - F f Un ■ - PA - i,J      i'J Аж [ \ i'J + 2Аж/      \ i'J 2АжЛ

• -    ^A [ o(u", + Q^i ) - G^U n - Q^i Yl + At S ( U n,\ (7) А у[_ V ^J   2А у/     V ^J    2А у)\            ’ ^J

На втором шаге предиктора производим усреднение:

u n+ ^1/2 — 2 ( u -y + U"j Y .                         (8)

На этапе корректора вычисляются окончательные значения газодинамических параметров на новом временном слое:

^u "j^{+ 1 —} U ", - аж ( F ⁱ⁺¹A J - F ^-1 / ² , ) -

• -    Ау ( ^G i’J+1/2 - ^G i,J-1/2 ^ + ^A^ ^S ( U ⁿ + ^1/2 ). ⁽⁹⁾

Потоки через границы ячеек можно вычислить различными способами, с использованием точного или приближенного решения задачи Римана:

F i+1/2,J — F ( u f+1/2,J , U” 1/2,J ).                               (10)

Одним из способов нахождения потоков через границы вычислительных ячеек является метод HLLC [6]. Во избежание загромождения формулами запишем его только для потоков в направлении оси ж (индекс j опускаем):

	' Fl ,
	F L + s L ( U * L - U L ),
F i+1/2 — *	F R + s R ( U * R - U R ),
	F r ,

0 <  s ^L , s ^L <  0 <  s, s <  0 <  s ^R , s ^R <  0;

F ^B = F ( U ^B ), В = {L,R};

1                       \

U * ^B = К

I

e^B , 2 +1/2 , о B ^{+ P} i+1/2

S v B Xi+1/2

( ^S - ^{u B} +1/2 ⁾ ^^{S +}

P i+1/2       \

P ^{L ( s B} - M^B 1/2 )

К =

B BB

B ^S ++1/2 ^p i+1/2    _s b - _S ;

^{s L} = min(u f₊1/2 - c f₊1/2 , U^ 1/2 - c f+1/2 ),                       ⁽¹⁵⁾

^S ^ = ^min(U i+1/2 + ^c i+1/2 , ^M R-1/2 ^{+ C} +-1/2 ),                          ⁽¹⁶⁾

P f+1/2 — Pt1/2 + P ^L +1/2 ^{b L} +1/2 ^{( s L} — ^{B L} +1/2 ) — p f+1/2 ^B f+1/2 ( ^{s R} — “ ^B . 2 ¹

^{P l} +1 / 2 ^{( sL} — ^M i+1/2 ⁾ — ^P+- 1 / 2 ^{( s} ^ — ^M +-1/2 )

eL  , -

^C 1+1/2 =

P i+1/2

A ^Y Ll — , \ ^P 1+1/2

r ^R  , -

^C 1+1/2 =

^R , R i+1/2 \ ^Y R^r   ^;

\   Pi+1/2

TT ^L     — тт^п+1/2 I ^Ax p

U i+1/2 = U i ⁺ 2 P^{- 1} ,

tt R _ y-r^n+1/2

U i+1/2 = ^U i+1

-

^{AX P 1 +1} ,

P 1 = r(

TT n

^U l+1

-

Ax

u n u n ^,

-

TT ⁿ U i - 1

Ax

) •

Для ограничения величины наклонов линейного распределения газодинамических переменных внутри вычислительных ячеек был использован простой и эффективный ограничитель ван Лира [3]:

2ab

C = a + b,

0,

ab >  0, ab <  0.

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

Постановка же граничных условий в области фонаря не столь однозначна. Перечислим возможные варианты:

1) условия свободного протекания:

• a)    обычные условия свободного протекания;

• b)    свободное протекание только из расчетной области (в нее ничего втекать не может);

2) гидростатическое равновесие:

• a)    свободное протекание через границу;

• b)    свободное протекание только из расчетной области (в нее ничего втекать не может).

Запишем уравнение гидростатического равновесия в виде

dp

= p g. dy

Дополнительным к условию (22) выберем постоянство температуры вдоль нормали к границе y = у _таж в области фонаря, что соответствует нашей изотермической модели (нижний индекс, соответствующий ж-координате, для простоты опустим):

Т м = Т м +1 = const.

Отсюда следует, что для идеального газа, которым мы моделируем среду, справедливо

соотношение

Р м    Р м +1       ,

— =----= const.

Р м    Р м+1

Запишем уравнение (22) в дискретном виде, заменив производную конечной разно-

стью:

Р м +1 — Р м    р м + р м +1 =             g.                            (25) Ay          2

Тогда из уравнений (24), (25) для значений плотности и давления в фиктивных ячейках нетрудно получить, что

gAy

^{Р м} + — ^{Р м}        Р м +1

^р м+1 = —---- , Р м +1 = ---- Р м

Р м   g^Ay            Р м

-- — --- ^р м ²

Для возможности газообмена с окружающей атмосферой можно задавать значения компонент скорости в фиктивных ячейках следующим образом:

^М м+1 = ^М м , У м +1 = М м .

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

и м +1 = и м ,

(им , и м >  о, [ -и м , ^и м <  0.

4.    Тестирование различных типов граничных условий

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

В прямоугольной расчетной области с размерами 60 х 30 м зададим равномерную сетку из 120 х 60 ячеек. Левую, правую и нижнюю границы области считаем твердыми стенками и задаем на них условия отражения потока. Верхняя граница позволяет газообмен с внешней средой. Для того чтобы результаты зависели только от типа условий на верхней границе, считаем печь выключенной (таким образом, условия одинаковы вдоль всей нижней границы).

Среду моделируем идеальным газом с показателем адиабаты у = 1,4, что соответствует двухатомному газу. Температура у нижней границы Т 0 = 20 ° C, давление — р ₀ = 101325 Па. Молярная масса газа М 0 = 28, 98 х 10 ^{- 3} кг/моль. Эти параметры соответствуют нормальным условиям в атмосфере у поверхности Земли.

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

Результаты тестирования различных типов граничных условий

Тип граничных условий	Полный импульс, кг ^ м/с	Полная кинетическая энергия, Дж
Свободное протекание	- 3 , 008 893 х 10 3	6 , 404 394 х 10 2
Свободное протекание с «диодными» условиями	- 3 , 011761 х 10 - 1	4 , 605 242 х 10 - 4
Изотермический гидростатический баланс, свободное протекание	- 1 , 287 595 х 10 - 1	2 , 302 684 х 10 - 4
Изотермический гидростатический баланс, свободное протекание с «диодными» условиями	- 1 , 287 593 х 10 - 1	2 , 302 681 х 10 - 4

Из полученных данных видно, что наименьшее отклонение состояния газа от равновесного обеспечили граничные условия изотермического гидростатического баланса на верхней границе. При этом разница между двумя вариантами — с «диодными» ограничениями и без них — минимальна.

Заключение

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

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

Авторы благодарны доктору физико-математических наук, профессору А.В. Хопер-скову за ценные замечания и полезные обсуждения.

ПРИМЕЧАНИЕ

• ¹    Работа выполнена при поддержке гранта РФФИ № 14-08-97044 р_поволжье_а.