CC BY
5
УДК 533.6.011.7
DOI 10.33764/2618-981X-2022-8-1-155-166
Исследование характеристик поднимающегося в газе нагретого плюма
И. Б. Палымский1,2*, В. И. Палымский3, И. В. Фролов4 1 Сибирский государственный университет телекоммуникаций и информатики, г. Новосибирск, Российская Федерация 2 Сибирский государственный университет геосистем и технологий, г. Новосибирск,
Российская Федерация 3 Институт теплофизики им. С. С. Кутателадзе СО РАН, г. Новосибирск,
Российская Федерация 4 Новосибирский государственный университет, г. Новосибирск, Российская Федерация
* e-mail: palymsky@yandex.ru
Аннотация. Путем численного моделирования исследуются характеристики поднимающегося нагретого плюма в сжимаемом, вязком и теплопроводном газе. Известно, что при докри-тической высоте области реализуется изобарный режим конвекции, а при закритической -адиабатический и сверхадиабатический. Показано, что при закритической высоте области определяющую роль играют адиабатические, а при докритической диффузионные процессы.
Ключевые слова: конвекция Рэлея-Бенара, газ, устойчивая стратификация, неустойчивая стратификация, температурный градиент
Investigation of the characteristics of a heated plume rise in a gas
I. B. Palymskiy1,2*, V. I. Palymskiy3,I. V. Frolov4 1 Siberian State University of Telecommunications and Information Sciences, Novosibirsk,
Russian Federation
2 Siberian State University of Geosystems and Technologies, Novosibirsk, Russian Federation 3 Kutateladze Institute of Thermophysics SB RAS, Novosibirsk, Russian Federation 4 Novosibirsk State University, Novosibirsk, Russian Federation
* e-mail: palymsky@yandex.ru
Abstract. The characteristics of an ascending heated plume in a compressible, viscous and heat-conducting gas are studied by numerical simulation. It is known that at the subcritical height of the region, the isobaric convection regime is realized, and at the supercritical one, adiabatic and superadiabatic. It is shown that adiabatic processes play a decisive role at the supercritical height of the region, and diffusion processes at the subcritical height.
Keywords: Rayleigh-Benard convection, gas, stable stratification, instable stratification, temperature gradient
Введение
Конвективные течения чрезвычайно распространены в технике и природе, поэтому конвекция и, в частности, конвекция Релея-Бенара хорошо исследована. Однако, исследования, как правило, выполнены в предположении несжимаемой
жидкости в приближении Буссинеска, для которого построена математическая модель, разработаны соответствующие численные методы, выполнены двумерные и трехмерные расчеты конвективных течений [1-3]. В основе приближения Буссинеска лежит предположение о малой вариации давления (гидростатического), благодаря чему зависимостью плотности среды от давления можно пренебречь и считать плотность функцией только температуры. Из сказанного следует, что при конвекции газа приближение Буссинеска может быть использовано только в лабораторных условиях, в областях достаточно малой высоты с относительно небольшим изменением гидростатического давления, а конвекция в областях с характерным размером несколько десятков сантиметров и более (17,3 см для воздуха при нормальных условиях) из-за относительно большого изменения гидростатического давления по высоте области требует учета сжимаемости газа на основе полных уравнений газовой динамики [4].
Необходимость рассматривать конвективные течения газовых смесей в областях большой высоты связана с наличием многочисленных практических приложений, возникающих, например, при рассмотрении вопросов взрывобезопас-ности в химической технологии при окислении углеводородов в реакторах, транспортировке их по трубопроводам и хранении в резервуарах [5].
В этой связи отметим, что конвекция газа на основе полных нелинейных уравнений газовой динамики из-за технических трудностей практически не исследована [6], что обусловлено большой жесткостью системы уравнений, что вынуждает проводить вычисления с очень малыми шагами по времени и крайне низким относительным изменением давления [6, 7]. Жесткость системы уравнений связана с наличием двух типов движения с качественно разными характерными временами, распространяющихся со звуковой скоростью быстрых термоакустических волн (аналог волн давления) [8] и развивающегося на их фоне относительно медленного конвективного движения.
Некоторым компромиссом являются упрощенные модели конвекции газа. С математической точки зрения структура таких систем аналогична структуре уравнений вязкой несжимаемой жидкости [7]. Полученные из уравнений газовой динамики в предположении малости числа Маха и параметра гидростатической сжимаемости соответствующие системы уравнений описывают среду, в которой звуковые возмущения распространяются с бесконечно большой скоростью. Такая система с успехом используется для моделировании конвекции в областях малой высоты при больших вариациях температуры и плотности [9], однако, расчеты конвективных течений в областях большой высоты, где в полной мере проявляется сжимаемость среды, необходимо проводить по полным нелинейным уравнениям газовой динамики [4].
В работе [10] расширена область применения уравнений Навье-Стокса на случай сжимаемых течений путем учета в уравнениях переноса тепла работы сил сжатия. Условие возникновения конвекции при этом принимает вид:
Яа ■ (1 - Бс) > Яаег
где Ra и Racr есть число Рэлея для сжимаемой среды и его критическое значение для несжимаемой в приближении Буссинеска, а Sc есть число Шварцшильда, равное отношению адиабатического градиента температуры к заданному. Поскольку предложенное соотношение не учитывает специфики конкретной задачи, а именно, геометрии области, граничных условий для температуры и скорости на боковых и горизонтальных границах, то точность предложенного соотношения может быть только качественной.
В работе [11] в линейном приближении аналитически и в нелинейном численно рассматривается конвекция газа в горизонтальном слое со свободными от касательных напряжений и недеформируемыми горизонтальными границами. Показано, что конвекция может развиваться при числе Шварцшильда строго меньше 1. В работах [4, 11] утверждается, что конвекция невозможна в адиабатическом режиме, и конвективное движение может развиваться только в сверхадиабатическом режиме при превышении температурным градиентом адиабатического.
Качественно возможность адиабатического подавления развития конвекции в сжимаемом газе показали многочисленные работы [11, 12] по численному моделированию конвективных течений, однако, из-за технических трудностей задаваемые значения определяющих безразмерных параметров в этих работах далеки от реальных, а именно, значение критерия гидростатической сжимаемости, характеризующее относительное изменение давления и обусловленное им относительное изменение плотности, завышено на три порядка.
На необходимость учитывать адиабатический градиент температуры c соответствующими поправками результатов измерений указывают многочисленные работы французской группы ученых, где проводилось исследование температурной конвекции в газообразном гелии при криогенной температуре 5 K. Показано, что, несмотря на относительно небольшую высоту области 0,2 м, поправка на сжимаемость газа оказывается очень существенной и изменяет теплообмен в несколько раз [13].
Роль сжимаемости газа становится определяющей в областях большой высоты и ее необходимо в полной мере учитывать даже в лабораторных масштабах, например, при исследовании асимптотических режимов конвекции при сверхбольших значениях числа Рэлея [14].
Но в настоящее время появились новые технические возможности для теоретического исследования конвекции в сжимаемом газе, связанные с возможностью проводить параллельные расчеты на многоядерных центральных процессорах и видеокартах.
При этом безусловная справедливость замечания о жесткости системы уравнений и необходимости проводить вычисления с неоправданно малыми шагами по времени частично нивелируется возможностью проводить вычисления по явной схеме с массивно-параллельной обработкой данных на многоядерных CPU или/и GPU с использованием технологий OpenMP или/и CUDA [15-17]. Преимущества массивно-параллельной обработки данных проявляются тем сильнее, чем больше объем информации обрабатывается, поэтому
такие вычисления особенно перспективны для проведения трехмерных расчетов и подробных двумерных.
В работе впервые на основе модели конвекции, основанной на системе уравнений газовой динамики, на плоскости Яа - АТ рассчитана нейтральная кривая, и получена диаграмма режимов конвекции. Показано, что при переходе высотой области критического значения наблюдается смена режима конвекции, а именно, изобарный режим конвекции (частным случаем которого является конвекция в приближении Буссинеска) сменяется адиабатическим либо сверхадиабатическим при достаточно большой разности температур.
Причем в непосредственно примыкающей к нейтральной кривой сверху узкой ножеобразной области (где реализуется адиабатический режим), благодаря адиабатическим процессам, конвекция может развиваться даже при устойчивой стратификации, что принципиально отличает конвекцию в сжимаемой среде от конвекции в несжимаемой [1-3, 18]. Полученная численно нейтральная кривая описывает зависимость критического числа Рэлея от высоты области и показывает адиабатическое подавление конвекции при большом ее значении.
Конвективное движение представляет собой совокупность поднимающихся и опускающихся плюмов (термиков), и, поскольку многие качественные свойства конвективного течения можно понять путем исследования характеристик движущегося плюма, это и обуславливает интерес к решению такой задачи.
Заметим также, что при специальной геометрии области и не очень высокой надкритичности конвекция развивается как валиковая квазидвумерная [18], что дает основание для рассмотрения сжимаемых и несжимаемых конвективных течений в двумерной постановке.
Данная работа состоит из двух частей.
В первой путем численного моделирования по полным уравнениям газовой динамики исследуются характеристики плюма, поднимающегося в горизонтальном слое сжимаемого вязкого и теплопроводного газа при докритической и за-критической высотах области. Цель данного этапа работы состоит в исследовании зависимости характеристик поднимающегося плюма от высоты области.
Во второй части при докритической и закритической высоте области в горизонтальном слое при подогреве снизу рассчитывается конвективное течение сжимаемого вязкого и теплопроводного газа. Цель данного этапа работы состоит в исследовании зависимости характеристик конвективного течения от высоты области.
Методика исследования
Конвективное движение в сжимаемом, вязком и теплопроводном газе, находящемся в поле силы тяжести, может быть описано следующей системой уравнений:
Яр, + Яр и + и • Ярх + V • Яру = М АЯр,
р, +р и + и •р х + v •р у = М А(р-рй),
1 4 1
и, + и • их + v • иу =--(рТ)х + М (- ихх + иуу +- vxy X
ур 3 ^ 3 '
1 4 1
V + и • ^ + v • ^ = —~ (рТ ) у + М (vxx + 3 ^у + 3 иху ) - СЕ ,
М у — 1 -
Т + и • Тх + v • Ту =—АТ — ^-Т и.
{ х у Рг у
Здесь м, ^ Р, р и Т - безразмерные компоненты скорости, давление, плотность и температура, М = у/((уТоЯ)0,5Я) = 4,608-10"8-Я-1 - число Маха, где вычисленная по кинематической вязкости скорость отнесена к адиабатической скорости звука, в качестве характерных значений для температуры и давления приняты То = 300 К и р = 1 атм., выбранные значения Я = 287 Дж/(кг-К), показателя адиабаты у = 1,4, кинематической вязкости V = 16• 10-6 м2/сек и числа Прандтля Рг = 0,71 соответствуют воздуху, а Се = gH/(yЯTо) = 8,130-10"5-Я - гидростатическая сжимаемость. В качестве масштаба длины выбрана высота области Н, температуры и плотности - их значения на нижней горизонтальной границе, скорости -адиабатическая скорость звука (уЯТз)0,5, давления - ЯР0Т0 и времени - Н/(уЯТ0)0,5. Зависимостью коэффициентов вязкости и теплопроводности от температуры прене-брегается. Высота области конвекции в расчетах варьировалась от 0,004 м до 0,5 м, отнесенный к вертикальному горизонтальный размер области равен п.
Поясним, что для упрощения моделирования, в духе приближения Бус-синеска в уравнение для плотности введен член, описывающий диффузию массы с учетом того, что в совершенном газе коэффициенты кинематической вязкости и диффузии массы по величине совпадают. Отметим, что подобное уравнение для плотности можно получить и в приближении Буссинеска [1], если учесть уравнение для температуры и линейную зависимость от нее плотности.
Для идентификации движущегося плюма в систему уравнений добавлено уравнение пассивной примеси Яр.
Все стенки области считались жесткими с условием прилипания для скорости и изотермическими, температура принимала значения Т0 и Т—ДТ для нижней и верхней горизонтальных границ, соответственно.
Температура на вертикальных границах, ее начальное и равновесное стационарное распределения считались линейными.
На всех границах области плотность, ее начальное и равновесное стационарное распределения принималась равной своим значениям в состоянии гидростатического равновесия, которое определялось согласно уравнению
(рТ ) у =—■ур Се , Т = 1 — АТ • у. (1)
Инициирование движущегося плюма проводилось в малой локальной области вблизи нижней границы, в которой задавалось Яр = 1 и малое значение вертикальной скорости, и, таким образом, инициированный плюм двигается по стационарному фону.
Расчеты проводились по явной схеме по времени [3, 19-21] с использованием массивно-параллельных вычислений, и, поскольку появление ударных волн в решении невозможно, использовалась недивергентная [20] запись системы уравнений с аппроксимацией конвективных нелинейных и диффузионных членов по монотонной схеме А. А. Самарского [3, 21], и, как следствие, используемый численный метод был первого порядка аппроксимации по времени и второго по пространству.
Все расчеты проводились на сетке (240x80) узлов с шагом по безразмерному времени 0,01. Проведенные тестовые расчеты на более подробных сетках по пространству и времени показали достаточность такого выбора.
Численное исследование было организовано по следующей схеме.
Сначала изучается движение плюма. Вычисленное распределение пассивной примеси (краски) позволяет определить координаты плюма и определить с помощью линейной интерполяции значения температуры, плотности и давления в его центре. Полученные данные, представленные как функции высоты области, затем сравниваются со стационарным распределением и с распределением, соответствующим адиабате Пуассона.
Затем рассчитывается конвекция в сжимаемом, вязком и теплопроводном газе. Цель данной серии расчетов установить зависимость основных интегральных характеристик конвекции от высоты области.
Результаты
Сначала опишем результаты исследования характеристик поднимающегося нагретого плюма.
Учитывая, что значение критической высоты для конвекции в рассматриваемой ситуации (в воздухе при нормальных условиях) равно 17,3 см, рассмотрим движение нагретого плюма, для определенности, при докритической высоте области Н = 0,1 м и закритической Н = 0,5 м.
Решение модельной задачи о подъеме нагретого плюма (конвективное движение можно представить как движение поднимающихся и опускающихся плюмов) в сжимаемой среде показывает качественно разные механизмы развития конвективного течения при докритическом и закритическом значениях высот.
На рис. 1 и 2 синей линией (1) показаны плотность и температура подымающегося плюма как функции безразмерной высоты у при Н = 0,5 м и надкритич-ности (отношение числа Рэлея к его критическому значению) г = 2, при этом красная линия (2) соответствует стационарному (начальному) распределению, а зеленая (3) - адиабате Пуассона.
О 0.2 0.4 V о.б
Рис. 1. Плотность подымающегося плюма при Н = 0,5 м
О 0.2 0.4 >' 0.6
Рис. 2. Температура подымающегося плюма при Н = 0,5 м
Видно, что при Н = 0,5 м (закритическая высота области) плотность и температура поднимающегося плюма приближенно следуют адиабате Пуассона, что указывает на квазиадиабатический характер течения. В этой связи заметим, что движение плюма рассчитывается с учетом вязкости и теплопроводности газа, и эти диффузионные процессы обуславливают лишь приближенное соответствие рассчитываемых величин адиабате Пуассона. Описанная картина качественно не меняется при увеличении надкритичности, и, таким образом, при закритической высоте области конвекции роль адиабатических процессов можно охарактеризовать как определяющую.
Принципиально другая ситуация реализуется при докритической высоте области Н = 0,1 м. На рис. 3 и 4 показаны плотность и температура подымающегося плюма при Н = 0,1 м и том же значении надкритичности г = 2, обозначения те же.
1.00010 [
р
1.00005
1
0 0.2 у 0.4
Рис. 3. Плотность подымающегося плюма при Н = 0,1 м
1
т
0.99995
0.99990 -1-1-*
0 0.2 ^ 0.4
Рис. 4. Температура подымающегося плюма при Н = 0,1 м
Сравнение плотности и температуры плюма при закритической высоте области (Н = 0,5 м, рис. 1 и 2) и докритической (Н = 0,1 м, рис. 3 и 4) показывает качественно разный механизм развития течения, а именно, при закритической высоте области плюм поднимается квазиадиабатически, а при докритической плотность и температура скорее близки к стационарным распределениям (красные линии), которые обусловлены чистой теплопроводностью без учета конвекции.
Не менее показательным является сравнение картинок поднимающихся плюмов при Н = 0,5 м (рис. 5) и Н = 0.1 м (рис. 6) при надкритичности г = 2. На рис. 5 и 6 показаны изолинии пассивной примеси, с помощью которой идентифицировались плюмы.
Изолинии пассивной примеси на рис. 5 и 6 соответствуют одному и тому же безразмерному времени ? = 104. Можно видеть, что при докритической высоте области Н = 0,1 м при движении плюма превалируют диффузионные процессы, это согласуется с данными на рис. 3 и 4, так как стационарные распределения обусловлены чистой теплопроводностью без учета конвекции.
Рис. 5. Поднимающийся плюм при Н = 0,5 м
Рис. 6. Поднимающийся плюм при Н = 0,1 м
Теперь несколько слов о зависимости интегральных характеристик конвекции от высоты области. Очевидно, что смена механизма развития конвекции при достижении высотой области критического значения обуславливает также и изменение степенных законов для основных интегральных величин. Физически это связано с тем, что при адиабатическом движении у частицы газа не происходит обмена теплотой с окружающей средой, и это обуславливает большую энергию частицы, увеличение кинетической энергии течения и теплообмена. То, что при закритической высоте области поднимающийся плюм медленнее теряет температуру и имеет большую кинетическую энергию, отчетливо видно при сравнении рис. 5 и 6.
Для примера рассмотрим число Нуссельта как функцию надкритичности г (рис. 7) при Н = 0,1 м (нижняя красная кривая (1)) и Н = 0,5 м (верхняя синяя кривая (2)), для определенности здесь рассмотрен двухвихревой стационарный режим конвекции. Знаком (3) показано значение, вычисленное с целью верификации приведенных данных на сетке (480x160).
Соответствующие аппроксимационные кривые (показаны пунктирами на рис. 7) принимают при Н = 0,5 м и Н = 0,1 м вид:
^0,5 = 2,25 • г0,390, МыоЛ = 1,38 • г0,274.
1 2 3 5 7 10 20 30 50 г 80
Рис. 7. Сравнение чисел Нуссельта как функций надкритичности
при разной высоте области
Отметим, что при одной и той же величине надкритичности г = 10 значения числа Нуссельта при Н = 0,1 м и Н = 0,5 м отличаются более чем в два раза при существенном отличии показателей степенных законов. Аналогичная ситуация наблюдается и для числа Рейнольдса, вычисленного по среднеквадратичной скорости, кинематической вязкости и высоте области.
Обсуждение
Конечно, сама по себе сжимаемость газа при высоте области Н = 0,5 м является ничтожно малой, и это создает впечатление, что ею можно пренебречь и считать среду несжимаемой. Однако это впечатление обманчиво, так как именно эта ничтожно малая сжимаемость при закритической высоте области и обуславливает появление новых режимов течения, где определяющую роль играет адиабатичность. В свою очередь, при адиабатическом движении у частицы газа не происходит обмена теплотой с окружающей средой, и это обуславливает ее большую энергию, увеличение кинетической энергии течения и теплообмена. Очевидно, что смена механизма развития конвекции при достижении высотой области критического значения обуславливает также и изменение степенных законов для основных интегральных величин. Причем отмеченное увеличение оказывается очень значительным, в несколько раз.
Заметим, что степень изменения гидростатического давления по высоте области определяется значением критерия весовой (гидростатической) сжимаемости Се = gH/(yR7o) = 8,13-10-5-Н, которое есть деленное на у относительное изменение гидростатического давления по высоте области, равное относительному изменению плотности из-за вариации давления. Подчеркнем, что смена режимов конвекции наблюдается при ничтожно малой относительной вариации гидростатического давления (и обусловленной им относительной вариации плотности) порядка двух тысячных долей процента от его характерного значения.
Заключение
Таким образом, при конвекции газа в областях высотой больше критической адиабатические процессы ( сжимаемость среды) играют определяющую роль, что находит свое отражение в степенной квазиадиабатической зависимости плотности всплывающего плюма от давления. Отсюда следует, что использовать модель конвекции в приближении Буссинеска при закритической высоте области неправомерно даже для получения качественной картины течения, и следует использовать более полную модель конвекции, основанную на уравнениях газовой динамики и учитывающую сжимаемость среды.
Благодарности
Выражаем благодарность РФФИ, при поддержке которого в виде гранта №20-08-00903-а была выполнена эта работа.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Гершуни Г.З., Жуховицкий Е.М. Конвективная устойчивость несжимаемой жидкости: монография. - М.: Наука, 1972. - 392 с.
2. Палымский И.Б. Турбулентная конвекция Рэлея-Бенара. Численный метод и результаты расчетов: монография. - Germany: LAP, 2011. - 232 с.
3. Пасконов В.М., Полежаев В.И., Чудов Л.А. Численное моделирование процессов тепло- и массообмена: монография. - М.: Наука, 1984. - 285 с.
4. Кригель А.М. Вопросы термодинамики турбулентной конвекции // Журнал технической физики. - 2016. - Т.86. - В.11. - С.136-139.
5. Взрывы газопровода в Гаосюне (2014) [Электронный ресурс]: доступ из ru.wikipedia.org.
6. Лапин Ю.В., Стрелец М.Х. Внутренние течения газовых смесей: монография. - М.: Наука, 1989. - 368 с.
7. Махвиладзе Г.М., Мелихов В.И. Численный метод исследования процессов медленного горения газов // Математическое моделирование. - 1989. - Т.1. - №6. - С. 146-157.
8. Зубков П.Т., Коваленко М.А. Симметрия и несимметрия в слое сжимаемого газа // Теплофизика высоких температур. - 2010. - Т.48. - №3. - С. 438-443.
9. Wan Zh.-H., Wang Q., Wang B., Xia Sh.-N., Zhou Q., Sun D.-J. On non - Oberbeck -Boussinesq effects in Rayleigh-Benard convection of air for large temperature differences // J. of Fluid Mechanics. - 2020. - V.889. - A.10. - P.1-21.
10. H. Jeffreys. The stability of a compressible fluid heated from below. Proc // Cambr. Phil. Soc. - 1930. - V.26. - P. 170-172.
11. Горбунов А. А., Полежаев В.И. Метод возмущений и численное моделирование конвекции для задачи Рэлея в жидкостях с произвольным уравнением состояния: препринт № 897. - М.: Институт проблем механики им. А.Ю. Ишлинского РАН, 2008. - 50 с.
12. Полежаев В.И. Численное решение уравнений Навье-Стокса для течения и теплообмена в замкнутой двумерной области: диссертация на соискание ученой степени к.т.н., 1967. - 196 с.
13. Chavanne X., Chilla F., Chabaud B., Castaing B., Chaussy J., Hebral B. High Rayleigh Number Convection with Gaseous Helium at Low Temperature // Journal of Low Temperature. -1996. - Vol. 104. - Is. 1/2. - P. 109-129.
