CC BY
30
Шафран Ю.В.1, Бутенко М.А.2, Кузьмин Н.М.,3 Хоперсков А.В.4
1 Волгоградский государственный университет, г. Волгоград, аспирант кафедры информационных систем и компьютерного моделирования, yury.shafran@yandex.ru
2 Волгоградский государственный университет, г. Волгоград, старший преподаватель
кафедры информационных систем и компьютерного моделирования,
butenkoma@gmail.com
3 Волгоградский государственный университет, г. Волгоград, доцент, к.ф.-м.н., доцент
кафедры информационных систем и компьютерного моделирования,
nmkuzmin@gmail.com 4 Волгоградский государственный университет, г. Волгоград, профессор, д.ф.-м.н., заведующий кафедрой информационных систем и компьютерного моделирования,
ahoperskov@gmail.com
Программное обеспечение для оптимизации системы вентиляции крупных промышленных цехов
КЛЮЧЕВЫЕ СЛОВА
Информационная модель, системы вентиляции, аспирационные течения, газодинамика, управление и оптимизация.
АННОТАЦИЯ
В работе обсуждаются возможности авторского научного программного обеспечения для решения задач оптимизации и управления системами вентиляции в крупных промышленных цехах на основе прямого газодинамического моделирования с применением параллельных технологий. Компьютерная модель учитывает различные режимы работы разнообразных вентиляционных и аспирационных устройств, металлургических печей и других источников тепла и воздуха. Предусмотрена возможность воздухообмена между помещением цеха и окружающей атмосферой. Программный комплекс предназначен для решения задачи оптимального размещения технических устройств, образующих вентиляционную систему, и управление ею. Оптимальное размещение вентиляционных устройств и эффективный режим их работы могут обеспечить наилучшие санитарно-гигиенические условия в рабочей зоне.
Введение
При построении системы вентиляции необходимо обеспечить соответствие определенным санитарным нормам большого набора параметров (поля температур и скоростей воздуха, его химический состав, фракционный состав аэрозолей, влажность). В реальных производственных помещениях имеется большое число разнообразных источников тепла, загрязняющих примесей, устройств воздухообмена. Возникает задача
проектирования и оптимизации системы вентиляции промышленных помещений с целью улучшения условий в рабочей зоне по температуре, подвижности воздуха, его химическому составу. Особую актуальность данная задача приобретает в случае крупных цехов, содержащих металлургические печи, прокатные станы, химическое производство. Решение этой задачи важно для химической, деревообрабатывающей, металлургической, цементной и многих других отраслей промышленности. Имеющиеся в настоящее время инструменты оптимизации основаны на простых инженерных расчётах [1], не способных адекватно описывать сложные динамические процессы и реальную геометрию помещений [2, 3, 4].
Решение описанной выше задачи представляется возможным при использовании численного газодинамического моделирования. Методы численного моделирования позволяют получать результаты точнее по сравнению с типовыми инженерными расчётами и дешевле по сравнению с натурными экспериментами. Математические модели и их компьютерная реализация для изучения течений воздуха с учётом переноса тепла и примесей внутри вентилируемых помещений на основе прямого газодинамического расчёта являются эффективным методом оптимизации и управления системами вентиляции [5, 6, 7, 8]. Предлагаемое программное обеспечение может использоваться для широкого круга задач управления микроклиматом в офисных и жилых зданиях, медицинских учреждениях с высокими требованиями к качеству воздуха и температурному режиму. Результаты вычислений позволяют получить достаточно полную и наглядную картину распределения параметров газа, в том числе, и загрязняющих примесей. Выполненные с помощью данного программного комплекса расчёты, в отличие от традиционно применяемых инженерных методик, основанных на интегральных оценках и полуэмпирических формулах, позволяют проанализировать работу сложных систем вентиляции и кондиционирования в трёхмерной постановке. Информационная модель
При проектировании программного комплекса использовались модульный подход и объектно-ориентированное программирование. Программный комплекс состоит из девяти основных модулей и блоков, подробное описание которых приведено в работе [9]. Здесь опишем только логику работы расчётного модуля нашего программного комплекса (рис. 1).
Расчётный модуль состоит из двух блоков: блока управления и вычислительного блока. Блок управления ответственен за задание начальных условий, расположение источников тепла, загрязнений и воздуха. С одной стороны, в нём учитываются возможность тепло- и массообмена с окружающей атмосферой и характеристики вентиляционных и аспирационных устройств. С другой стороны, в этом же блоке учитываются факторы, связанные с геометрией задачи: граничные условия, геометрическая форма производственного помещения,
расположение в нем различных технологических устройств и их возможные перемещения. Здесь же реализовано построение расчётных сеток.
Блок управления - main
Загрузка входных данных
" Граничные условия
Построение произювд ПОМСи геометрии тненного цения
¥ч#Т р;ьс';пгу:ж|!г|ин ТАШопо™ч«>ма устройств Учет рек***» райоты вцмиицрт yinpflAcT* пеиьчиенир;
Построение расчетных се™
Начальные условия
Источники теппэ
Источники загрязнений
Источнивм геппймаеенсбмена
тмпо-
млсал&лн ( J-UOCtWPÍU
ЦенГЧПЯииСин&О и аЯ1и|Мцнвнчыв
Выходные данные
Поле температур
Турбулизация воздуха
Поле скоростей
Динамика примесей
Химические примеси
Аэрозоли и пыль
Рис.1. Диаграмма, описывающая принцип работы расчётного модуля Вычислительный блок обеспечивает прямое численное интегрирование уравнений нестационарной газодинамики с учётом источников тепла, загрязнений и воздуха. Реализованная в нём численная схема описана ниже. Она достаточно эффективна для данного класса задач и хорошо поддаётся распараллеливанию. На выходе расчётный блок выдаёт результаты в виде пространственного распределения газодинамических параметров, таких как температура, плотность воздуха, поле векторов скоростей, распределение химических примесей и т.п.
Отметим ключевые особенности разработанного программного комплекса:
1. Расчётный модуль можно запускать отдельно на любом локальном компьютере при наличии непосредственного физического или удалённого доступа;
2. Количество расчётов для отдельно взятой машины ограничено только системными ресурсами;
3. Возможность работать с изменяющимися во времени трёхмерными данными аналогично 4D-подходу, рассмотренному в работе [10].
Математическая модель
В основе математической модели лежит использование полной системы уравнений динамики сжимаемого вязкого газа с учётом силы тяжести. Используемые уравнения могут быть представлены в виде:
ш+ат(и) + тс(и) + аи(и) = ^ док (и)+^ )
^ дх Ту & ^ дк '
где х,у и г - оси координат в декартовой системе, к е (х, у г) и
( Р Л
и =
Р»х Р^у
Е
(2)
у
где р - плотность, их, иу и иг - компоненты скорости воздушного потока и Е -полная энергия единицы объёма. Операторы F(U), G(U), H(U), D(U) и S(U) в уравнении (1) задаются векторами
( по Л ( РП Л ( по Л
Р(и) =
Р»х
Р»х2 + Р
Р»у»х Р»Пх К»х (Е + Р)
, С (и) =
7
Р»у
Р»х»у
2 , Р»у + Р
Р»Пу
°у (Е + Р)
, Н(и) =
7
Р»х»г
РПу»2
2 , Р°, + Р
Ко, (Е + РI
(3)
О к (и) =
( 0 ^ ( о 1
(кх 0
(ку , 8(Ц) = 0
Рg
°х(х +оу(ку +о,(к, у
(4)
где р - давление, д - модуль ускорения свободного падения и Ту - элементы тензора вязких напряжений:
2 (дох ТПу до, Л
дх ду д,
[
(
дох до у Л х + у
ду дх
у
( дПу дох Л —- +—х
дх ду
(до* дох [ —- + —х
К дх д,
2 (дПу дох до, Л
ду дх д,
3
(о , до, , д, ( доу
( до, доу Л
(
2
дх
до,
д, ду
до, дох
доу Л
д, дх ду
(5)
ду д,
П - коэффициент сдвиговой вязкости.
Для замыкания системы уравнений используется уравнение состояния идеального газа в виде
Р
£ = ■
(6)
р{У- 1)'
где % - внутренняя энергия и у - показатель адиабаты. Используя уравнение (6), полную энергию Е в уравнениях (2) и (3) можно записать следующим образом:
т
(Г~ 1) 2 1
Численная схема
Для упрощения расчёта система уравнений (1) с использованием техники расщепления по физическим процессам может быть разделена на две задачи [11]:
Ти + 5Р-и) +ТС-и) +ТН-и) = 0 (8)
д1 дх ду & ( )
и
дИ ^дБ к -И) _/ттч
Покроем расчетную область регулярной сеткой, состоящей из прямоугольных ячеек. Тогда для решения системы уравнений (8) можно использовать численную схему М^^, имеющую второй порядок точности [12]. Для трёхмерного случая новые значения состояния могут быть вычислены следующим образом [11]
Ип+2 = р —Н—Н д.с— (Ип), (10)
где № - состояние на шаге п и ЛЬ - шаг по времени, определяемый из стандартного для газодинамики условия Куранта-Фридрихса-Леви. Указанный в уравнении (10) порядок применения операторов по осям позволяет сохранить второй порядок точности при обобщении исходной одномерной численной схемы на многомерный случай.
В качестве примера промежуточного шага схемы М^^ по одной из
осей приведём формулы вычисления нового состояния по оси у:
( \
- (11)
И= Ик + — 1 1 —у
С ! - О !
и 1—, к г,. н—, к V 2 ' 2 У
где Ду - размер ячейки по оси у, индексы i, j и к показывают номер ячейки в трёхмерной сетке и Gi>¡+%>k - потоки вещества, импульса и энергии между ячейками с индексами (Цк) и (1]'+1,к), вычисленные посредством решения задачи Римана методом HLLC
С, 1+1 к = Снььс -И,'-м;Ии1+1к), (12)
2'
по приближённым значениям на границе двух ячеек. Эти значения вычисляются следующим образом:
1 И -Ду .)-С. (13)
1 И 2ДУ О-И .)-С-И1)], (14)
2
И1= Ик+~., (15)
и1 = и__
и г,],к и г,/,к -
где - приближённая производная значений состояния для ячейки с
индексами (Цк), скорректированная ограничителем MINMOD
( Л
А . }.к = MINMOD
..! к;А
г,/—, к г,/ н—, к
V 2 ' 2 У
А. .. = <
г,} ,к
тах
min
А 1, = и - и I. ,к
( Л
0;min А 1 ; А 1
г, /—,к г, /н— ,к
V 2 2 У
с
V
0;тах А 1 ; А 1
г,/—,к г, /н—,к
. V 2 ' 2 У
А
г, /н—,к 2
(17)
(18)
(19)
< 0
г, /н—,к 2
Для решения уравнения (9) используется схема Рунге-Кутты второго порядка [11]. При последовательном вычислении производных скоростей и элементов тензора вязких напряжений используются центральные разностные аппроксимации производных, например:
(дту, Л (ху \.к )
*У 'г, /--,к
ду ( ди* Л
Уг, / ,к
ду
"У
У г, /н—,к 2
"У
(20)
(21)
Для обеспечения второго порядка точности по времени используется следующая последовательность при вычислении нового временного слоя:
& Дt & Дt
и п+2 = R 2 F " G" Н " R 2 R 2 н" G" F " R 2 (ип)' где R - решение уравнения (9) методом Рунге-Кутты. Пример моделирования
При построении модели использовались параметры электросталеплавильного цеха №2 ВМЗ «Красный Октябрь» [13, 14]. В качестве примера рассмотрим конвективный столб, формирующийся от печи ДСП-200. Для очистки воздуха и улучшения температурного режима в рабочей зоне рассмотрим различные режимы работы воздушного зонта (рис. 2), который вытягивает воздух. Для отсечения горячего конвективного воздушного столба рассмотрим аэратор, создающий горизонтальные воздушные потоки, обеспечивающие более эффективную работу зонта.
Сталеплавильная печь моделировалась цилиндром радиуса 5 м и высоты 4 м. Активная зона радиусом 1,1 м расположена на верхней поверхности. Температура активной зоны 1250°С, скорость восходящего воздушного потока 1 м/с. Граничные условия пассивной зоны соответствуют модели без скольжения: все компоненты скорости
(22)
0
воздушного потока равны нулю.
Активная зона аэратора состоит из трёх компонент: к центральной части высотой 0,5 м и шириной 3 м с обеих сторон примыкают боковые, высота и ширина которых равны 1,5 м.
Активная зона вытяжного зонта квадратной формы со стороной 4 м втягивает воздух со скоростью 2,5 м/с. Пассивные части аэратора и вытяжного зонта моделируются аналогично пассивной зоне печи.
Вытяжной шит
Аэрнгор
и
о
1,5
9 к'ч ь
10
Л1ратор } 1,5
ЙЫ1ЯЖНОИ
МЖТ
(а)
(6)
Рис. 2. Расположение объектов в расчётной области: (а) - в плоскости YZ; (б) - в плоскости XY Числами обозначены соответствующие размеры в метрах Эффективность работы воздушного зонта будет характеризоваться интегральной величиной внутренней энергии воздуха, отводимой вытяжным зонтом:
* = *, (23)
где ип - скорость потока, нормальная к активной поверхности вытяжного зонта и 5п - площадь сечения ячейки плоскостью, параллельной поверхности вытяжного зонта. Суммирование производится по расчётным ячейкам, граничащим с активной зоной вытяжного зонта. Энергия Е сильно зависит от параметров остальных устройств.
В качестве примера, рассмотри результаты расчётов для различных скоростей потока и„, создаваемого аэратором: 2,5 м/с, 5 м/с и 10 м/с (рис. 3).
При малых скоростях влияние аэратора незначительно и мы имеем небольшое отклонение горячего столба воздуха в сторону зонта (рис. 4). С увеличением скорости иа эффективность работы зонта растет, поскольку увеличивается величина Е вместе с перехватываемым горячим воздухом. Однако, с ростом иа увеличиваются расходы на вентиляцию, и растёт подвижность воздуха в рабочей зоне, приводя к нарушениям санитарных норм. Эти требования (пунктирная линия на рис. 3) приводят к ограничениям на параметры вентиляционных устройств. Таким образом,
использование данного программного обеспечения помогает находить наиболее оптимальные режимы работы вентиляционных систем.
5.25 ^Е, МДж/с
4.75--
4.5-
5-
4.2$
V, м/с
2.5
5
7.5
10
Рис. 3. Зависимость внутренней энергии, поглощаемой вытяжным зонтом в единицу времени от скорости потока, создаваемого аэратором
Работа выполнена при поддержке Российского фонда фундаментальных исследований (грант 14-08-97044 р_поволжье_а).
1. Хрусталев Б.М. Теплоснабжение и вентиляция. - М.: АСВ, 2008. - 784с.
2. Хоперсков А.В., Азаров В.Н., Хоперсков С.А., Коротков Е.А., Жумалиев А.Г. Формирование нестационарных режимов при моделировании аспирационных течений: неустойчивость Кельвина-Гельмгольца // Вестник ВолГУ, Сер. 1: Математика. Физика. - 2011. - Т14, №1. -С.151-155.
3. Логачев И.Н., Логачев К.И. Аэродинамические основы аспирации. - Санкт-Петербург: Химиздат, 2005. - 659с.
4. Аверкова О.А., Логачев И.Н., Логачев К.И. Аэродинамика противопылевой вентиляции. -2012. - 370с.
5. Raji A., Hasnaoui M., Bahlaoui A. Numerical study of natural convection dominated heat transfer in a ventilated cavity: Case of forced flow playing simultaneous assisting and opposing roles // International Journal of Heat and Fluid Flow. - 2008. - Vol.29. - Pp.1174-1181.
6. Zhao Y., Liu D., Tang G.-F. Multiple steady fluid flows in a slot-ventilated enclosure // International Journal of Heat and Fluid Flow. - 2008. - Vol.29. - Pp.1295-1308.
7. Conceifao E.Z.E., Gomes J.M.M., Lucio M.J.R. Numerical simulation in a complex topology building using HVAC systems equipped with PMV control // International Journal of Advanced Engineering Applications. - 2014. - Vol.7, Iss.4. - Pp.53-63.
8. Farea T., Ossen D., Alkaff S., Ghadikolaei F., Dodo Y. Numerical Simulation of Convective Natural Ventilation in a Light Well Connected to a Horizontal Void: Comparison of Various RANS-Based Turbulence Models // Journal of Basic and Applied Scientific Research. - 2014. - Vol. 4, Iss. 7. -Pp.218-230.
9. Бутенко М.А., Бурнос Д.В., Хоперсков С.А., Холодков В.С., Морозов А.Г. Информационная модель программного комплекса для оптимизации и управления системами вентиляции на основе прямого газодинамического моделирования // Вестник ВолГУ. Сер. 10: Инновационные технологии. - 2012. - №6. - C.31-37.
10. Храпов С.С., Кобелев И.А., Писарев А.В., Хоперсков А.В. 4D-модели в задачах экологического моделирования: проектирование информационной системы // Вестник ВолГУ Сер. 10: Инновационные технологии. - 2011. - №5. - C.119-124.
11. Toro E.F. Riemann Solvers and Numerical Methods for Fluid Dynamics. A Practical Introduction. -
Литература
2006.
1999. - 624p.
12. Berton C. Why the MUSCL-Hancock scheme is Ll-stable // Numerische Mathematik. Vol.104. - Pp.27-46.
13. Хоперсков А.В., Шафран Ю.В., Бутенко М.А. Численное моделирование вентиляционных течений в промышленных помещениях // Южно-Сибирский научный вестник. - 2014. -№2(6). - С.98-102.
14. Butenko M., Shafran Y., Khoperskov S., Kholodkov V., Khoperskov A. The optimization problem of the ventilation system for metallurgical plant // Applied Mechanics and Materials. - 2013. -Vol.379. - Pp.167-172.
4J-4Л-
35-Я»-
2S-»
№
5»
id 25
Ш
JO-
Ji 30-ii' 5№ I S-Kb
s-
О
J(t-35' 3U-25' 2U-15' 10-i' 0-
лк-
>S' .uv
25 20-li Lfr
15 M (Г)
15 JO ic)
Рис. 4. Распределение температуры при разных скоростях потока аэратора: (а) и (б) при 2,5 м/с, (в) и (г) при 5 м/с, (д) и (е) при 10 м/с; (а), (в) и (д) - в плоскости YZ, проходящей через центр печи, (б), (г) и (е) - в плоскости ХУ, проходящей через центр
аэратора
