CC BY
42
11. Hussein H.M.S. Transient investigation of a two phase closed thermosyphon flat plate solar water heater // Energy Conversion and Management. - 2002. - № 43. - P. 2479-2492.
12. Desrayaud G., Fichera A., Marcoux M. Numerical investigation of natural circulation in a 2D-annular closed loop thermosyphon // International Journal of Heat and Fluid Flow. - 2006. - V. 27. -№ 1. - P. 154-166.
13. Liaqat A., Baytas A.C. Numerical comparison of conjugate and nonconjugate natural convection for internally heated semi-circular pools // International Journal of Heat and Fluid Flow. - 2001. - V.22. - № 12. - P. 650-656.
14. Кузнецов Г.В., Аль-Ани М.А., Шеремет М.А. Численный анализ влияния температурного перепада на режимы переноса энергии в замкнутом двухфазном цилиндрическом термос-
ифоне // Известия Томского политехнического университета. - 2010. - Т 317. - № 4. - С. 13-19.
15. Kuznetsov G.V., SitnikovA.E. Numerical analysis of basic regularities of heat and mass transfer in a high-temperature heat pipe // High temperature. - 2002. - V. 40. - № 6. - P. 898-904.
16. Роуч П. Вычислительная гидродинамика. - М.: Мир, 1980. -616 с.
17. Самарский А.А. Теория разностных схем. - М.: Наука, 1977. -656 с.
18. Пасконов В.М., Полежаев В.И., Чудов Л.А. Численное моделирование процессов тепло- и массообмена. - М.: Наука, 1984. - 288 с.
Поступила 02.02.2011 г.
УДК 53.082.2:550.3
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ РАБОТЫ КОНВЕКТИВНОГО ДАТЧИКА ПРИ ДЕЙСТВИИ ЦЕНТРОБЕЖНОЙ СИЛЫ
И.А. Бабушкин, В.А. Демин, Д.В. Пепеляев
Пермский государственный университет E-mail: demin@psu.ru
Проведено численное моделирование отклика конвективного датчика инерционных ускорений на действие центробежной силы. Описано основное течение в виде пульсирующего конвективного факела и его кризис при определенной частоте вращения. Показано, что непрерывное изменение величины центробежной силы дает возможность постепенно наклонять конвективный факел в плоскости широких граней полости. С помощью воздействия центробежной силы выяснено оптимальное расположение термопар в рабочей камере датчика сучетом наличия осложняющих факторов, связанных с неизбежным несовершенством экспериментальной установки.
Ключевые слова:
Тепловая конвекция, ячейка Хеле~Шоу, центробежная сила, конвективный факел, численное моделирование.
Key words:
Thermal convection, HeleShaw cell, centrifugal force, convective tail, numerical modeling.
Введение
Вращение в большинстве случаев значительно усложняет поведение конвективных систем [1-3]. Даже самое простое вращательное движение с постоянной угловой скоростью отражается натече-нии вследствие добавления к действующим массовым силам дополнительных сил инерции. Действие этих сил (кориолисовой и центробежной) приводит к тому, что даже в типичных условиях тепловая конвекция по причине присущих ей объемных неоднородностей полей скорости и температуры имеет, как правило, трехмерный характер. Если конвективная система находится вблизи или на оси вращения, и величина угловой скорости невелика, то центробежной силой чаще всего можно пренебречь в силу ее малости по сравнению с силой Кориолиса [4]. Однако бывает так, что центробежная сила является неотъемлемым дополнительным фактором, влияющим на конвекцию, или ее действие технологически оказывается востребованным. Например, когда некоторый процесс для своей реализации требует постепенного монотон-
ного изменения подъемной силы в условиях фиксированного нагрева. Такая необходимость возникает при калибровке приборов, рабочим телом которых является неоднородно нагретая жидкость [5, 6].
Активное управление величиной результирующей подъемной силы можно организовать за счет постепенного «включения» центробежной силы. В [6, 7] представлены результаты экспериментов и проведено теоретическое обоснование идеи создания датчика на основе ячейки Хеле-Шоу, позволяющего регистрировать ограниченные по времени вибрационные сигналы и восстанавливать их исходные характеристики. В экспериментах особое внимание уделялось тестированию прибора на возможность возникновения различных нежелательных режимов конвекции, влияющих на показания. В частности, исследования выявили наличие дополнительных функциональных возможностей, позволяющих датчику регистрировать продолжительные по времени монотонные инерционные воздействия. Для определения новых возмож-
ных областей применения датчика было решено провести эксперименты по влиянию вращения на конвективные течения, что, в свою очередь, повлекло за собой необходимость детального теоретического исследования устойчивости конвективных режимов в ячейке Хеле-Шоу при действии центробежных сил.
Постановка задачи
Рабочая полость представляет собой прямоугольный параллелепипед, верхняя и нижняя узкие грани которого играют роль нагревателя и холодильника. По причине конструктивных особенностей датчика подогрев осуществляется снизу точечным источником тепла. Температура нагревателя периодически меняется в зависимости от времени, что создает течение в полости в виде пульсирующего конвективного факела чрезвычайно чувствительного по отношению к внешним инерционным воздействиям. В опытах датчик устанавливался в корзине вращательной установки на достаточно большом расстоянии от оси вращения, так что вклад центробежной силы в результирующее ускорение был существенным. Конвективная камера жестко прикреплялась к корзине, поэтому при повороте вращательной установки на угол 2 ^полость совершала один полный оборот вокруг своей оси.
С целью минимизации действия силы Кориолиса и выполнения приближения плоских траекторий, датчик всегда крепился так, чтобы широкие грани ячейки Хеле-Шоу были ориентированы параллельно штанге. Выберем систему координат так, чтобы ось г была перпендикулярна широким граням. В этой системе координат /(¡апрсоза, созр, -трта) -единичный вектор, направленный вертикально вверх (рис. 1). Полость целиком заполнена жидкостью и вращается с угловой скоростью Пеу в плоскости, которая по причине несовершенства экспериментальной установки слегка наклонена под углом р к горизонту. В опытах величина угла наклона плоскости вращения контролировалась и не превышала 1°. Наклон плоскости вращения приводил в экспериментах к возникновению нежелательных эффектов, однако, как оказалось, полностью исключить этот фактор невозможно.
г = ±1:
= -а(Т -Т,),
(1)
Относительно конвекции в ячейке Хеле-Шоу отмечается высокая чувствительность характеристик течения к тепловым условиям на широких гранях [8], поэтому соответствующее краевое условие должно быть максимально реалистичным. Необходимость учета теплопроводности массива, прилегающего к широким граням, приводит к следующему граничному условию для отклонения температуры жидкости Т от среднего по времени линейного профиля [9]:
дТ_
"&
где Т] - температура окружающей среды; а=Хп/Хр - коэффициент теплоотдачи.
В экспериментах конвективная камера имела высоту к=32 мм, ширину /=24 мм и толщину 2й=4 мм. В качестве рабочей жидкости был выбран гептан. Ранее в работах [6, 7] было показано, что физические свойства гептана заметно более предпочтительны при сравнении его с другими жидкостями. Широкие грани полости были окружены массивом из полиметилметакрилата (ПММА). Выберем за единицу длины полутолщину слоя. Тогда, для фактической толщины массива Б=6 мм, коэффициентов теплопроводности гептана и ПММА 1=0,12 Вт/(м.К), Ят=0,20 Вт/(м.К) коэффициент теплоотдачи в безразмерных единицах равен а=0,54. Оставшиеся условия для Тна узких гранях не так сильно влияют на течение и могут быть приближенно взяты, соответственно, как для теплоизолированных и идеально теплопроводных стенок:
х = 0.
т: дТ
: дх
= 0; у = 0, Н: Т = 0.
(2)
В случае строго вертикального расположения оси вращения на элемент жидкости в рабочей камере действует центробежная сила с ускорением, приблизительно одинаковым для всей полости, т. к. г>>/. В неинерциальной системе отчета, связанной с кюветой, вектор результирующего ускорения у-0[0хг ] подъемной силы, будучи постоянным по направлению и величине, лежит в плоскости широких граней. Результирующая массовая сила ориентирована к горизонту под некоторым углом, значение которого зависит от угловой скорости вращения. Меняя угловую скорость, можно плавно регулировать величину и направление подъемной силы. Если плоскость вращения слегка наклонена, то вектор у в общем случае уже не параллелен широким граням. В системе отсчета, связанной с полостью, вектор у прецессирует вокруг единичного орта у. Строго говоря, это нарушает приближение плоских траекторий и должно приводить к появлению поперечной компоненты скорости. Однако в опытах угол наклона р был мал, а следовательно малы у и у. Прямой расчет течения в поперечном сечении, вызванного покачиванием вектора у в плоскости (у, г), показал, что при амплитуде качаний 1° в полости
возникает длинноволновое колебательное движение, амплитуда скорости которого на два порядка меньше, чем те характерные скорости, которые имеют место в плоскости широких граней при пульсационном нагреве снизу. В результате эффектом прецессии вектора у в дальнейшем будем пренебрегать, сохранив тем самым приближение плоских траекторий.
В ходе экспериментов обрабатывался сигнал, представляющий собой разность температур между симметричными спаями термопары вблизи нагревателя в зависимости от угловой скорости вращения. Определялась амплитуда сигнала и его среднее значение, отнесенные к разности температур между теплообменниками. Оказалось, что существует такое значение угловой скорости вращения (0~1 1/с), при котором как амплитуда, так и среднее значение разности температур имеют экстремум. Результатов термопарных измерений для объяснения этих ярко выраженных пиков на графиках оказалось недостаточно. Для получения полной картины было проведено численное моделирование конвективных течений в плоскости широких граней полости в математической постановке, наиболее близко отражающей натурный эксперимент.
Основные уравнения и методика расчета
Для теоретического описания конвективных течений будем использовать систему уравнений тепловой конвекции в приближении Буссинеска [4]. Особенности движения полости в совокупности с приближением Хеле-Шоу (к,/>>й) позволяют пренебречь влиянием силы Кориолиса на течение в плоскости широких граней. Для характерной скорости в плоскости широких граней у~1 мм/с, плеча г~500 мм и угловой скорости 0~1 1/с отношение кориолисовой силы к центробежной и силе тяжести равно /к/^ц=^/0г~2.10-3; Fк/Fт=vQ/,g~10-4. Оценка показывает, что действующая на течение сила Кориолиса пренебрежимо мала по сравнению с центробежной силой и практически не нарушает приближение плоских траекторий. Уравнения тепловой конвекции с учетом только центробежной силы в безразмерной форме имеют вид:
= -Чр + ЛV+Ra Ту + Та0х[0хг]Т, (3)
Pr — + (VУ)Т = ЛТ, divV = 0, (4)
дt
где V, Т, р - поля скорости, температуры и давления. При обезразмеривании выбирались следующие единицы: полутолщина полости й (длина), й/у (время), (скорость), 0 (температура), рУ/й2 (давление). Здесь 0 - средний по времени перепад температур на верхней и нижней гранях; р - плотность жидкости; у, % - коэффициенты кинематической вязкости и температуропроводности. Си-
стема уравнений (3), (4) явно содержит три безразмерных параметра
Pr =у/%, Ra = gрt0d3 /у%, Ta = Pt0O 2й 4 ¡у% ,
соответственно, числа Прандтля, Рэлея и Тэйлора (у - ускорение силы тяжести, Д - коэффициент теплового расширения). Дополнительными безразмерными параметрами в задаче являются длина полости Ь, высота Н, угол Д, характеризующий наклон плоскости вращения и коэффициент теплоотдачи а. Для скорости на границах полости будем использовать условие прилипания. На нижней границе полости поддерживается распределение температуры, меняющееся с течением времени по закону
Т|у=0 = 2в~(х-т2)%т2(©0.
Эта зависимость моделирует пульсирующий с течением времени локальный подогрев, который создается точечным нагревателем в середине нижнего теплообменника. Такой нестационарный нагрев повышает чувствительность конвективной системы. Считается, что частота пульсаций температуры на нижней границе полости достаточно велика по сравнению с характерным временем рассасывания тепловых возмущений на широких боковых гранях. Далее в решении для температуры выделим осредненный линейный профиль, соответствующий подогреву снизу Т=1-у/Н. Эта модель, скорее всего, оправдывается в реальных условиях и неоднородности профиля температуры на широких гранях не будут сильно отражаться на форме и динамике конвективного факела.
Ограничения на толщину ячейки позволяют использовать приближение плоских траекторий, поэтому дальнейшее рассмотрение будет проводиться на основе уравнений, записанных в терминах функции тока и температуры. Функция тока связана с компонентами скорости соотношениями ух=5Т/5у, уу=-5Т/5х. Зависимость функции тока от координаты г моделировалась тригонометрической функцией, чтобы удовлетворялось условие прилипания на широких гранях. Поле температуры описывалось суперпозицией базисных функций, характерных для двух предельных случаев: идеально теплопроводных и теплоизолированных широких граней (ниже Т - отклонение температуры от осредненного линейного профиля):
^ (х, у, г, 0 = у( х, у, 0 соэ(жг/2),
Т(х, у, г, 0 = в(х, у, 0 + &(х, у, 0соэ(жг/2), (5)
где у, в, & - зависящие от времени амплитуды, характеризующие распределение полей функции тока и температуры в плоскости (х, у). Амплитуды в и & связаны между собой, т. к. выражение для температуры в (5) должно удовлетворять граничному условию (1). Разложения (5) подставлялись в исходную систему, после чего уравнения усреднялись по г в соответствии процедурой Галеркина. В терминах амплитуд вихря скорости р, функции тока у и температуры в эти уравнения имеют вид:
др | 8 (ду др ду др ^ =
дt ЗпРп
-Ra
ду дх 4этД
дх ду
\~Л1р-^р-
пН
х| ^^соэД-■двsinрcos0t
дх ду
+г Та
4
пН
2(2+а) дв_
п ду
(6)
„ дв 2 (3п2+3п2а+8а2)
Рг----1-- --------------— х
дt Зп (п2+8а+2а2)
'дудв дудв")_Л в п(2+а) ;
ду дх дх ду \ 1 (п2+8а+2а2) ”
ду а (2+а)п2
в,
дх 2(п2+8а+2а2) ( )
где г0 - безразмерное плечо, р=Л1у, Л1 - плоский оператор Лапласа. Уравнения (6), (7) совместно с граничными условиями (2) образуют краевую задачу. В предельных случаях а^да и а^0 имеем переходы, соответственно, к идеально теплопроводным и теплоизолированным широким граням.
Расчеты проводились методом конечных разностей по явной схеме. Компьютерный модуль был написан на языке программирования Рог1гап-90. Задача решалась в переменных у и р, т. е. использовался двухполевой метод [10]. Основные расчеты выполнялись на сетке 51:41. При составлении конечно-разностного аналога уравнений тепловой конвекции пространственные производные аппроксимировались центральными разностями, а производные по времени - односторонними. Значения вихря на границах полости находились по формулам Тома, которые получались разложением функции тока в ряд Тейлора в приграничной точке с точностью до квадратичных членов. Уравнение Пуассона для поля функции тока решалось методом простых итераций. В ходе расчетов использовался метод установления; определялись мгновенные значения полей у и Т, находились
максимальное и минимальное значения функций тока. Для моделирования показаний термопар в трех узлах с координатами (26,32), (18,16) и (34,16) определялась температура, и анализировалась амплитудно-частотная характеристика временной зависимости 5(0=Т(34,16)-Т(18,16).
Обсуждение результатов
Для сопоставления результатов численного моделирования и экспериментальных данных безразмерные параметры принимали значения, наиболее близкие к действительности: Рг=6,9 (гептан), Ь=12, Н=16, г0=500. Угол р, определяющий наклон плоскости вращения к горизонту и характеризующий несовершенство экспериментальной установки, в соответствии с оценками лежал в интервале 0...2°. Частота пульсаций нагревателя ограничивалась сверху естественным временем кондуктивно-конвективного теплоотвода, и при пересчете в безразмерные единицы для нее была установлена величина ш=0,17. Средняя по времени разность температур между верхним и нижним теплообменниками была равна 0=0,9 °С, что соответствовало числу Рэлея Ёа«190. Угловая скорость вращения менялась в диапазоне 0=0,3...1,3, обуславливая верхнюю границу значений для числа Тэйлора Та<0,05.
При строго вертикальном положении полости в отсутствие вращения конвективный факел над нагревателем характеризуется лево-правой симметрией. Приведение системы во вращение нарушает симметрию полей скорости и температуры: за счет действия центробежной силы теплый конвективный факел наклоняется в плоскости широких граней в сторону оси вращения. В дополнение небольшой наклон плоскости вращения создает колыхания конвективного факела вблизи среднего отклонения. На рис. 2 представлены поля функции тока (а, в) и температуры (б, г) в фиксированный момент времени для угла наклона плоскости вращения р=1°. В случае кратных частот возможен параметрический резонанс, однако в эксперименте подобные ситуации специально не рассматривались. Две частоты со и О в общем случае считались несоизмеримыми. С другой стороны, в случае
Рис. 2. Поля функции тока и температуры при 0=0,4 (а, б) и 0=1,2 (в, г) для а=0,54
а=0 в расчетах для кратных частот при 0=0,51 и 0=0,68 (0=3о, О=4о) была зафиксирована смена режима колебаний, сопровождавшаяся резким изменением характеристик течения (пики на рис. 3, б). Если значения частот сопоставимы, но несоизмеримы, то в системе наблюдается сложное колебательное движение. Разность температур в двух симметричных точках 5(0 совершает колебания с усредненной амплитудой А вблизи некоторой среднего значения (5(0). Фурье-анализ временной зависимости 5(0 показывает, что наибольшей амплитудой обладает гармоника с частотой 2о.
Сначала при увеличении угловой скорости вращения амплитуда А и среднее значение (5(0) начинают быстро расти по модулю (рис. 3). Отсутствие гладкости у графиков, наблюдающееся в большей степени в теории, и менее заметное в эксперименте, имеет несколько причин. Во-первых, начальная фаза состояния нагревателя по отношению к положению кюветы не контролировалась в эксперименте и поэтому случайным образом выбиралась в расчетах, т. е. отслеживались в определенном смысле осредненные эффекты, не зависящие от момента начального приготовления системы. Как следствие при моделировании этой ситуации колебания, складываясь в разных фазах, дают немонотонное изменение амплитуды с ростом 0. Во-вторьа, применительно к сложному колебательному процессу имела место значительная ошибка при вычислении А и (5(0) за счет грубости алгоритма определения этих величин. Также, вследствие больших характерных тепловых и гидродинамических времен релаксации в ячейке Хе-ле-Шоу могло наблюдаться некоторое неустано-вление колебательного режима.
Несмотря на наличие трудноустранимых погрешностей, из рис. 3 отчетливо видно, что при достижении некоторой частоты А и (5(0) графики для а=0,54 и 0 имеют ярко выраженные максиму-
мы. Визуализация полей показывает, что при этой частоте происходит смена конвективного течения. Центробежная сила становится настолько большой, что конвективный факел опрокидывается, и в полости устанавливается трехвихревое колебательное течение (рис. 2, в) с одним большим вихрем, занимающим всю центральную часть полости и простирающимся по диагонали из левого нижнего в правый верхний угол. Наряду с главным вихрем в левом верхнем и правом нижнем углах полости всегда существуют два сателлитных вихря меньших размеров с противоположной закруткой. Условия смены конвективного режима определяются многими факторами, в том числе углом наклона плоскости вращения. Однако в первую очередь опрокидывание конвективного факела зависит от угловой скорости вращения. В предельном случае ¡3=0° циклическая частота О исчезает из краевой задачи (6), (7), (2), в результате чего колебания в значительной степени регуляризуются. Тем не менее, качественно зависимости А и (5(0) от 0 остаются прежними, т. е. при определенной частоте эти графики имеют ярко выраженный максимум, который свидетельствует о смене конвективного режима.
Сравнение опытных данных и результатов численного моделирования проводилось по следующим признакам: наличие максимумов у зависимостей А и (5(0), совпадение А и (5(0) с экспериментальными значениями по величине и близость расчетной частоты к экспериментальному критическому значению, при котором происходит смена режима. Зависимости А и (5(0) для а=0,54 и 0 (ПММА и теплоизолированные широкие грани) демонстрируют хорошее согласие с опытом. По выше перечисленным характеристикам они приблизительно одинаковы и значительно дистанцируются от другой группы графиков для «=110 (сталь) и а^да. В случае идеально и высоко тепло-
Рис. 3. Амплитуда А и среднее значение (5(1)) в зависимости от частоты О; 3=1°: О, О - а=0,54 (ПММА), 110 (сталь); □, х - а=0,да; 3=0°: А, + - а=0,54, 110
проводных широких граней, характерный максимум на графике (5(0) отсутствует. Натеплопро-водных широких гранях все возмущения быстро рассасываются. По этой причине амплитуда колебаний разности температур в жидкости лежит значительно ниже, чем для а=0 или 0,54.
Иными словами, широкие грани ПММА, в роли регулятора конвективного движения, по своим свойствам находятся значительно ближе к теплоизолированным, чем к идеально теплопроводным стенкам. С другой стороны, для теплоизолированных широких граней на графике (5(0), в отличие от реалистичного случая а=0,54 на резонансных частотах, отчетливо видна смена конвективного режима (пики на рис. 3, б). Конвективная система в случае теплоизолированных широких граней более чувствительна к различного рода тепловым возмущениям, которые не рассасываются на стенках, а, напротив, развиваются на фоне периодического воздействия во вторичный колеба-
СПИСОК ЛИТЕРАТУРЫ
1. Бубнов Б.М., Голицын Г.С. Конвективные режимы во вращающейся жидкости // Доклады АН СССР. - 1985. - Т. 281. -№3.- С. 552-555.
2. Бердников В.С., Захаров В.П., Марков В.А. Тепловая гравитационно-центробежная конвекция в подогреваемом снизу слое жидкости // Инженерно-физический журнал. - 2001. - Т. 74. -№4. - С. 111-115.
3. Ингель Л.Х. Конвекция во вращающейся среде над термически неоднородной горизонтальной поверхностью // Известия РАН. Механика жидкости и газа. - 2008. - № 6. - С. 25-35.
4. Гершуни Г.З., Жуховицкий Е.М. Конвективная устойчивость несжимаемой жидкости. - М.: Наука, 1972. - 392 с.
5. Бабушкин И.А., Богатырев Г.П., Глухов А.Ф., Путин Г.Ф., Авдеев С.В., Бударин Н.М., Иванов А.И., Максимова М.М. Изучение тепловой конвекции и низкочастотных микроускорений на Орбитальном комплексе «Мир» с помощью датчика «Дакон» // Космические исследования. - 2001. - Т. 32. -№2.- С. 150-158.
тельный режим. При а=0,54 метод продолжения по параметру дает возможность избежать развития вторичного колебательного течения даже при точном попадании на резонансную частоту, что хорошо согласуется с экспериментом.
Выводы
По совокупности признаков результаты расчетов, учитывающие специфику материала, из которого изготовлены широкие грани, наиболее адекватно отражают картину конвективных течений и перестройку режимов в ячейке Хеле-Шоу. Выбранная теоретическая модель может использоваться для описания конвекции и теплопереноса в рабочей камере датчика. Показано, что геометрия полости позволяет использовать физические возможности центробежной силы для управления тепломассопереносом при имитации действия продолжительных по времени и постоянных по величине инерционных сигналов.
6. Бабушкин И.А., Глухов А.Ф., Демин В.А., Дягилев РА., Мало-вичко Д.А. Сейсмоприемник на основе ячейки Хеле-Шоу // Прикладная физика. - 2008. - № 3. - С. 134-140.
7. Бабушкин И.А., Демин В.А., Пепеляев Д.В. Принципы регистрации инерционных сигналов с помощью конвективных датчиков // Известия Томского политехнического университета. - 2010. - Т. 317. - № 4. - С. 38-43.
8. Бабушкин И.А., Демин В.А. Вибрационная конвекция в ячейке Хеле-Шоу. Теория и эксперимент // Прикладная механика и техническая физика. - 2006. - № 2. - С. 40-48.
9. Мызникова Б.И. Численное исследование надкритических режимов тепловой конвекции в ячейке Хеле-Шоу при подогреве снизу // Исследование тепловой конвекции и теплопередачи. - Свердловск: УНЦ АН СССР. - 1981. - С. 23-31.
10. Тарунин Е.Л. Вычислительный эксперимент в задачах свободной конвекции. - Иркутск: Изд-во Иркут. ун-та, 1990. - 225 с.
Поступила 02.02.2011 г.
