Численный код для расчета аспирационных течений в промышленных цехах Текст научной статьи по специальности «Физика»

www.volsu.ru

001: http://dx.doi.Org/10.15688/jvolsu1.2015.5.4

УДК 004.94, 519.63 ББК 22.193

ЧИСЛЕННЫЙ КОД ДЛЯ РАСЧЕТА АСПИРАЦИОННЫХ ТЕЧЕНИЙ В ПРОМЫШЛЕННЫХ ЦЕХАХ1

Николай Михайлович Кузьмин

Кандидат физико-математических наук, доцент кафедры информационных систем и компьютерного моделирования, Волгоградский государственный университет nmkuzmin@gmail.com, infomod@volsu.ru

просп. Университетский, 100, 400062 г. Волгоград, Российская Федерация

Мария Анатольевна Бутенко

Старший преподаватель кафедры информационных систем и компьютерного моделирования, Волгоградский государственный университет butenkoma@gmail.com, infomod@volsu.ru

просп. Университетский, 100, 400062 г. Волгоград, Российская Федерация

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

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

ю

о

сч

^ Введение

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

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

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

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

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

р(у) = Ро ехр

»4

(1)

р(у)

Г ~Мо\ё\ \

Ро ехр\ ~кг0Т у),

(2)

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

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

ди д¥ д С з

дЬ дх ду

(3)

где

и

( р \ ри ри

V е

¥

( ри \

ри2 + р рии V (е + р)и У

С

/

ри

рии 2

\

ри2 + р

\ (е + р)и У

Б

0 0

рд

\ р^ у

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

Р

У -1

+

рV2

(4)

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

е

2

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

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

В расчетной области {хтЫ < х < хтах, утЫ < у < утах] введем неподвижную эйлерову сетку

** = + (г - 1г = 1 2,■■■,L, АХ = - Хтт)/1, (5)

У] = УтЫ + - ^ Ау, Э = 1, 2,...,М, Ау = (утах - Ушгп)/М. (6)

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

и*. = Цп.- —

^ ^ Аж

Е +

1,3

р

г,3

2Ах) V г'3 2Ах)

\

АЛ

А

ЧЦ-+2А£) - с(ц« - Ш

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

и:

га+1/2

1,3

ЦЦ*. + им.

2 \ ьи ~ ьи I

+ А1Ъ(Ц11). (7)

(8)

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

1 = Ц'Ъ - А (Fi+1/2,j - ^г-1/2,з

)

- А/ {^^+1/2 - ^ЬЗ-1/2

^+АЩЦпг+1/2). (9)

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

Ег

г+1/2,3

Е(Ц

и

К

1+1/2,з, Ц г+1/2,3

).

(10)

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

Е

+1/2

Еь, 0 < вь,

Еь + вь(Ц*ь - Ц), зь < 0 < 5, ЕК + 8К(Ц*К - ЦК), в < 0 < вК,

¥К, вК < 0;

(11)

ь

рБ = ), в = {LR}.

тВ

(12)

U

*В

К

В

еВ / рв \ I Р^ + (s -<1/2)vs+Р -Ъ))

г+1/2

Pi+1/2

^ = РВ

оВ _ 7<В

В й "г+1/2.

-1/2 В

SB - S

S = тт(ит/2 — Cj+1/2 ,«¿+1/2 — Сг+1/2) = min(wf+1/2 + cf+1/2, «¿+1/2 + cf+1/2)

(13)

(14)

(15)

(16)

1

s

Pi+1/2 Pi+1/2 + Pi+1/2Mi+1/2(s Мг+1/2) Pi+1/2Mi+1/2 (s Мг+1/2) Pi+1/2(sb — Mi+1/2) — Pi+1/2(si — Mr-1/2)

(17)

rL

4+1/2

\

^г+1/2 Р +1/2

U i

U +1/2

иГ+'/2 + ^ Рг,

nR

4+1/2

U i

U +1/2

\

^¿+1/2 Р +1/2

Un+1/2 — P U г+1--Г" P

Аж,

+1 ,

(18) (19)

Рг

£('u?+1 - и? и? - U-

Аж

Аж

О-

(20)

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

2аЬ

С = < а"ТЪ, 0, (21)

0, аЬ < 0.

s

3. Особенности постановки граничных условий

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

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

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

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

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

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

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

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

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

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

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

% = П. (22)

dy

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

Тм = Тм+i = const. (23)

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

^ = = const. (24) Рм Рм+1

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

Рм+1 — Рм Р м + Р м+1 /огч

-;-=-Q. (25)

Ay 2 v '

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

, 9АУ„

Рм + тРм Рм+1 (26)

Р м+1 = —-Рм+1 =-Рм. (26)

Рм 9АУ Рм

Рм 2

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

им+1 = им, Ум+1 = им. (27)

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

им+1 = им, Ум+1 = < (28)

- Ум, Ум < 0.

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

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

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

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

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

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

Тип граничных условий Полный импульс, кг-м/с Полная кинетическая энергия, Дж

Свободное протекание -3,008 893 х 103 6,404 394 х 102

Свободное протекание с «диодными» условиями -3,011761 х 10-1 4,605 242 х 10-4

Изотермический гидростатический баланс, свободное протекание -1,287595 х 10-1 2,302 684 х 10-4

Изотермический гидростатический баланс, свободное протекание с «диодными» условиями -1,287593 х 10-1 2,302 681 х 10-4

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

Заключение

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

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

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

ПРИМЕЧАНИЕ

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

СПИСОК ЛИТЕРАТУРЫ

1. Программное обеспечение для оптимизации системы вентиляции крупных промышленных цехов / Ю. В. Шафран, М. А. Бутенко, Н. М. Кузьмин, А. В. Хоперсков // Современные информационные технологии и ИТ-образование. — 2014. — Т. 1, № 1 (9). — C. 509-517.

2. Leer, B. van. On the relation between the upwind-differencing schemes of Godunov, Engquist-Osher and Roe / B. van Leer // SIAM Journal of Scientific Statistic Computing. — 1984. — Vol. 5, № 1. — P. 1-20.

3. Leer, B. van. Towards the ultimate conservative difference scheme II. Monotonicity and conservation combined in a second order scheme / B. van Leer // Journal of Computational Physics. — 1974. — Vol. 14, № 4. — P. 361-370.

4. Leer, B. van. Towards the Ultimate Conservative Difference Scheme V. A Second Order Sequel to Godunov's Method / B. van Leer // Journal of Computational Physics. — 1979. — Vol. 32, № 1. — P. 101-136.

5. The optimization problem of the ventilation system for metallurgical plant / M. A. Butenko, Yu. V. Shafran, S. A. Khoperskov, V. S. Kholodkov, A. V. Khoperskov // Applied Mechanics and Materials. — 2013. — Vol. 379. — P. 167-172.

6. Toro, E. F. Restoration of the contact surface in the HLL-Riemann solver / E. F. Toro, M. Spruce, W. Speares // Shock Waves. — 1994. — Vol. 4, № 1. — P. 25-34.

REFERENCES

1. Shafran Yu.V., Butenko M.A., Kuzmin N.M., Khoperskov A.V. Programmnoe obespechenie dlya optimizatsii sistemy ventilyatsii krupnykh promyshlennykh tsekhov [Software for Optimization of System of Ventilation of Big Industry Workshops]. Sovrеmеnnyе informatsionnyе tеkhnologii i IT-obrazovaniе, 2014, vol. 1, no. 1 (9), pp. 509-517.

2. Leer B. van. On the Relation Between the Upwind-Differencing Schemes of Godunov, Engquist-Osher and Roe. SIAM Journal of Scientific Statistic Computing, 1984, vol. 5, no. 1, pp. 1-20.

3. Leer B. van. Towards the Ultimate Conservative Difference Scheme II. Monotonicity and Conservation Combined in a Second Order Scheme. Journal of Computational Physics, 1974, vol. 14, no. 4, pp. 361-370.

4. Leer B. van. Towards the Ultimate Conservative Difference Scheme V. A Second Order Sequel to Godunov's Method. Journal of Computational Physics, 1979, vol. 32, no. 1, pp. 101136.

5. Butenko M.A., Shafran Yu.V., Khoperskov S.A., Kholodkov V.S., Khoperskov A.V. The Optimization Problem of the Ventilation System for Metallurgical Plant. Applied Mechanics and Materials, 2013, vol. 379, pp. 167-172.

6. Toro E.F., Spruce M., Speares W. Restoration of the Contact Surface in the HLL-Riemann Solver. Shock Waves, 1994, vol. 4, no. 1, pp. 25-34.

NUMERICAL CODE FOR COMPUTATION OF ASPIRATION FLOWS IN INDUSTRIAL BUILDINGS

Nikolay Mikhaylovich Kuzmin

Candidate of Physical and Mathematical Sciences, Associate Professor, Department of Information Systems and Computer Simulation, Volgograd State University nmkuzmin@gmail.com, infomod@volsu.ru

Prosp. Universitetsky, 100, 400062 Volgograd, Russian Federation

Mariya Anatolyevna Butenko

Senior Lecturer, Department of Information Systems and Computer Simulation,

Volgograd State University

butenkoma@gmail.com, infomod@volsu.ru

Prosp. Universitetsky 100, 400062 Volgograd, Russian Federation

Abstract. The design of numerical code for gas-dynamics simulation of aspiration flows was presented. The choice of method for integration of gas-dynamics equations was made. Several possibilities for boundary conditions of hydrostatic case was described.

For simulation of aspiration flows two-dimensional approach was used. In this case, computational domain is rectangle, that corresponding to vertical slice of industrial building. Left, right and bottom boundaries are solid borders. For them we used reflected boundary conditions. On the top boundary gas free flow is available.

There are several possibilities of boundary conditions on the top boundary:

1) free flow:

(a) obvious free flow;

(b) free only outflow from computational domain (but not inflow into);

2) hydrostatic equilibrium:

(a) free flow;

(b) free only outflow from computational domain (but not inflow into).

We perform numerical simulations for all listed variants. The results showed that minimal change of initial hydrostatic balance was given by isothermal hydrostatic boundary conditions. And difference between two their variants is very small.

The two-dimensional numerical code on the base of MUSCL approach and isothermal hydrostatic balance satisfying boundary conditions was developed by us. It is possible to support hydrostatic equilibrium in computational domain with very high accuracy. We think, that this code may be used for tasks of projecting and optimisation of ventilation systems of industrial building with success.

Evolution of this code may be in two directions: the full three-dimensional code and isothermal hydrostatic balance boundary conditions of two order accuracy.

Key words: aspiration flows, industrial building, gas-dynamics, numerical simulation, boundary conditions.