2011
ВЕСТНИК ПЕРМСКОГО УНИВЕРСИТЕТА
Математика. Механика. Информатика Вып. 4(8)
УДК 532.5.013.3
Эффекты асимметричных вращательных колебаний в задаче о тепловой конвекции в полости квадратного сечения с подогревом снизу
А. Б. Мелентьев
Пермский государственный национальный исследовательский университет Россия, 614990, Пермь, ул. Букирева, 15
a.b.melentyev@mail.ru; +7 906 889-89-49
Методами математического моделирования исследованы эффекты, возникающие при вращательных колебаниях полости квадратного сечения в случае тепловой конвекции жидкости, подогреваемой снизу. Вращательные колебания были вызваны поворотами полости относительно внутренней оси. Были рассмотрены как симметричные, так и асимметричные колебания. Выявлены пороговые значения параметров, при которых происходит смена направления течения для случая колебательных и стационарных режимов. Выяснено, что отклонение от симметрии колебаний приводит к смещению пороговых значений параметров.
Ключевые слова: конвекция; вращательные колебания; асимметрия.
Введение
В конвективных задачах важным является учет вибраций и других модуляций параметров, так как эти воздействия ведут к различным эффектам, особенно в случае подогрева снизу [1-3]. Вибрации неоднократно исследовались теоретически и экспериментально. Сравнительно часто рассматривался высокочастотный предел, позволяющий использовать метод осреднения. Модуляции низкой частоты исследовались значительно реже. Кроме того, как правило, исследовалась модуляция по гармоничному закону. В то же время исследования [4] показали, что асимметрия модуляций приводит к ожидаемым эффектам.
В работе численно исследуются конвективные режимы в бесконечном цилиндре квадратного сечения под действием симметричных и асимметричных (кусочно-синусои-
© А. Б. Мелентьев, 2011
дальных) низкочастотных изменений угла поворота цилиндра вокруг собственной оси.
1. Постановка задачи
Область квадратного сечения, для которой решалась задача, изображена на рис. 1.
Здесь 0Х^о - инерциальная система отсчета; 0ху - неинерциальная система отсчета, связанная с полостью; а - угол отклонения 0ху от 0Х^о.
Для описания движения жидкости использовались уравнения свободной конвекции в приближении Буссинеска [1]. Особенностью
постановки является добавление слагаемых в правой части уравнения для скорости, соответствующих силам, возникающим при поворотах системы: силе Эйлера, центробежной силе и силе Кориолиса [5-7] при допущении независимости плотности от температуры в этих слагаемых (отсутствие значительного влияния эффекта теплового расширения на указанные силы) вследствие малых частот модуляций угла. ду
— + (у У)у = —Ур + Рг Ау +
дt
+ Gr Т(1 sinа(t) + j cosа(t)) -—(ё^) х г + ё(^ х (ё^) х г) + 2ё^) х У), (1.1)
— + (у У)Т = АТ, дt
div у = 0.
Здесь и = (н,и,0) - скорость жидкости, а - угол поворота полости, ё = (0,0, а) -
угловая скорость вращения полости, ё - про-и зводная от угловой скорости по времени, г = (х, у ,0) - радиус-вектор компоненты
жидкости, Т - поле температуры, р - давление, 1, j - единичные векторы, направленные вдоль осей х и у соответственно.
Уравнения записаны в безразмерных переменных. Масштабам расстояния, времени, скорости, давления и температуры соответствовали: Ь, Ь / %, %/ Ь, р%2 / Ь,
^/Г2/у2 (X - коэффициент температуропроводности, р - плотность, V - кинематическая вязкость, в - разность температур между нижней и верхней границей полости).
На границах полости задано условие прилипания скорости и значения температуры:
= 0,Т\г = 0.5 - у.
(1.2)
Задача зависит от числа Прандтля Рг = V / % , числа Грасхофа Gr = gP0^3 / V2, а также от параметров, определяющих поведение зависимости угла от времени а (і):
а(і) =
,лі
а0 cos(—), і є [0, і1]; і1
,п(і - і,).
-а0 cos(---------—), і є [і1, і0].
і0 - І1
(1.3)
Здесь а0 - амплитуда, t0 - полный период модуляций, tl- время первой доли периода, при котором а^) изменяется от а0 до —а0. Параметр асимметрии [4] равен
(1.4)
2. Метод решения
Задача решалась численно двухполевым методом [1, 2, 8] в переменных функции тока / и вихря скорости р :
дш дш
и = V =--------, V = V =-----------.
1 ду у дх
р = (гоіи) = -Аш.
(2.1)
Соответствующие уравнения в неинер-циальной системе отсчета 0ху (см. рис. 1) имеют вид
др дш др дш др
— +-------------------= Рг Ар +
ді ду дх дх ду
дТ дТ
+ Gг(—cosa(і)--------sinа(і)) - 2а(і),
дх ду (2.2)
Аш + р = 0,
дТ дш дТ дш дТ
---+ —----------------= АТ.
ді ду дх дх ду
Здесь а (і) - вторая производная по
времени от а (і). В расчетах в качестве точки
начала координат использовалась левая нижняя точка полости, а не центр полости, как на рис. 1. Для решения системы (2.2) использовалась явная схема с центральными разно-
стями с порядком аппроксимации 0(т + h2).
Для большинства расчетов шаг сетки равнялся 1/20. При таком шаге погрешность интегральных характеристик была менее 3%. Алгоритм решения соответствовал [2], в цикле по времени:
1) на первом этапе находились значения температуры Т и вихря скорости р для внутренних узлов по явной схеме;
2) на втором этапе методом ПВР [2, 9] находились значения функции тока у/ для внутренних узлов;
3) на третьем этапе по формуле Тома [2] находились значения для вихря на границе.
V
г
Шаг по времени определялся по формуле, с запасом, обеспечивающим устойчивость в исследуемом интервале параметров [2]:
h2
Т = -------:--
12 Рг + h /
| • т
В расчетах типичное значение т < 2 • 10-1. Решение задачи Дирихле для функции тока методом ПВР на каждом временном шаге проводилось с точностью 10—5.
Счет по времени продолжался до установления стационарного максимального значения функции тока с точностью 10—6. Колебательный режим считался установившимся при относительном изменении амплитуды на соседних периодах максимума функции тока менее 10—6.
3. Результаты вычислительных экспериментов
В расчетах менялись значения числа Грасхофа Gr и параметры модуляции: амплитуда а 0, полный период t0 и параметр асимметрии % . Число Прандтля было фиксированным Рг = 1.
3.1. Стационарные значения угла наклона
Сначала были рассмотрены стационарные конвективные режимы при постоянном угле наклона для определения параметров бифуркации. Как будет показано, эти параметры бифуркации позволяют предсказать эффект модуляции. Аналогичный прием применения параметров бифуркации был использован в задаче управления равновесием [10].
Начальные условия соответствовали равновесному значению температуры Т0 = 1 — у для а = 0 и возмущению вихря
скорости: р0 = Аху(1 — х)(1 — у). (3.1)
Типичное значение А в формуле (3.1) при расчётах составляло ±100. В дальнейшем использовался метод продолжения по параметру, в котором при смене параметров использовалось решение, полученное при предыдущих значениях параметров.
Для угла наклона, равного а = 5°, были получены стационарные решения для различных значений числа Грасхофа Ог < 20 • 103. Результаты этих решений представлены зависимостью экстремальных значений функции тока от числа Грасхофа на рис. 2.
Рис. 2. Зависимость экстремума функции тока от числа Грасхофа (стационарное решение)
Штриховая линия (линия 3) соответствует тривиальной бифуркации при угле наклона = 0 . Критическое значение числа Грасхофа, соответствующего значению = 0 для квадратной области, Огх ~ 5.1 • 103 [1, 2].
Кружками отмечены значения функции тока, соответствующие "благоприятной" ветви, на которой вращение жидкости происходит против часовой стрелки. Треугольниками обозначены отрицательные экстремальные значения функции тока на "неблагоприятной" ветви, на которой основная масса жидкости вращается по часовой стрелке. Выход на неблагоприятную ветвь осуществляется при определенных условиях. Для выхода на эту ветвь можно в начальном условии для вихря скорости (3.1) задать закрутку по часовой стрелке выбором знака А , а затем задать значение , уменьшающее интенсивность этого течения.
Отметим, что при смене знака угла наклона "благоприятная" и "неблагоприятная" ветви меняют свои знаки. При этом картина ветвления является зеркально симметричной относительно горизонтальной оси /т = 0 для случая а с другим знаком.
Подобные бифуркационные зависимости, описываемые складкой Уитни, рассматривались, например, в [12, 13]. В работе [12] причиной, нарушающей равновесие, был угол поворота нагрева относительно вертикали для полости кругового сечения. В работе [13] причиной, нарушающей равновесие, было движение верхней границы. При малых значениях параметра, нарушающего равновесие, теория [13] подсказывает описание этой зависимости кубической параболой:
/т — С (Ог — Ог )/т — с2Ог = 0. (3.2)
Обработка результатов расчета методом наименьших квадратов (далее - МНК) позволила получить коэффициенты этой
зависимости c1 « 0.004, c2 « 0.0023 . Ветви
этой зависимости на рис. 2 изображены сплошными линиями как для "благоприятной" (линия 3), так и для "неблагоприятной" (линия 5) ветвей. Согласно теории при Gr > Gr * существует 3 решения, соответсву-ющих корням кубического уравнения. Два из них устойчивы относительно бесконечно малых возмущений, а третье нустойчиво.
Неустойчивая ветвь (линия 6) отмечена пунктирной линией.
Важной характеристикой ветвления является точка (у*, Gr*). Для коэффициентов
с1 и е2, найденных методом МНК, эти величины равны у* « -2.142 и Gr* « 8.587 • 103.
Обсудим точность описания ветвления кубическим уравнением. Как известно [13], для кубического уравнения при Gr = Gr * отношение модулей амплитуд "благоприятной" ветви к "неблагоприятной" равно 2. Соответствующее отношение амплитуд при Ог = 8.65 • 103 превышает ожидаемое значение 2 на ~ 18%. Это свидетельствует о значительной погрешности описания амплитудной зависимости кубическим уравнением при а = 5°.
Поэтому для получения характеристик (у*, Ог*) ветвления в зависимости от угла а были выполнены следующие вычислительные эксперименты (далее - ВЭ). В этих ВЭ после получения решения на "неблагоприятной" ветви уменьшалось значение Ог с шагом АОг = 50 . Последнее значение, при котором реализуется "неблагоприятная" ветвь, давало приближенные значения характеристик ветвления: у* « -2.23, Ог* « 8.65 • 103. Отличие этих значений от полученных из зависимости (3.2) менее 5%.
На рис. 3 представлены изолинии функции тока и температуры для "благоприятной" и "неблагоприятной" ветвей при одном и том же значении числа Грасхофа Ог = 8.65 • 103, близким к Ог * .
1)
2)
-0.2 0 1 2 3 4
Рис. 3. Изолинии 1) функции тока и 2) температуры для а) "благоприятной" и б) "неблагоприятной ветви" (Ог = 8.65 • 10 , а = 5°)
Как видно, искривление изотерм на "благоприятной" ветви значительно больше, чем на "неблагоприятной"; соответствующие тепловые потоки примерно равны 2.174 и 1.519.
Обсудим особенности перестройки течения при переходе с "неблагоприятной" ветви на "благоприятную" при значениях числа Грасхофа, близких Ог * . Для таких значений числа Грасхофа время перестройки велико, при удалении от Ог * время перестройки уменьшается. Иллюстрацией к отмеченной особенности является рис. 4.
6 4 2
о -2 -4 -6
-У'т і
-
-
\ : * - і
\( -- 2 .... з
- ^2 І
0 С 5
Рис. 4. Зависимость максимума функции тока от времени
На этом графике показана зависимость экстремума функции тока от времени для четырёх значений числа Грасхофа: Gr = 8.65 • 103, 8.6 • 103, 8.55 • 103, 8.45 • 103 и а = 5°. До момента і = і1 « 0.467 значение числа Грасхофа было фиксированным
Gr = 8.65 • 10 и а = 0 . При і > і1 было установлено значение угла а = 5° и к моменту і = і2 « 2.072 было получено стационарное решение на "неблагоприятной" ветви. При
і > і2 представлены решения для трех значений числа Грасхофа, уменьшенных на 50 (линия 1), 100 (линия 2), 200 (линия 3). Интервалы времени с момента і2 до момента смены
знака экстремума функции тока соответственно равны 1.635, 0.926, 0.591.
В дальнейших ВЭ определялись значения а *, приводящие к перестройке течения с "неблагоприятной" ветви на "благоприятную", для различных значений числа Грасхофа в интервале от 6 -103 до 11.3 -103. Эти ВЭ позволили определить зависимость а * (Ог). Обработка этой зависимости методом МНК позволила представить ее в виде линейной зависимости:
а* = 1.27S • 103(Gr - Gr0).
(3.3)
Отличие расчетных значений от зависимости (3.3) на указанном интервале чисел Грасхофа не превосходит 10%.
3.2. Модуляция угла наклона
Рассмотрим конвективные режимы при колебании угла наклона (1.3). В этом случае при поиске установившихся режимов использовались нулевые начальные условия, так как наклон полости сам вырабатывал возмущения определенного знака. Колебательные режимы жидкости рассматривались без учета Шлих-тинговского механизма [3, 14], так как использовались низкочастотные колебания, при которых толщина пограничного слоя
(3.4)
сравнима с характерным размером полости.
Обсудим вначале колебательные режимы при наличии симметричной (£ = 1) модуляции угла наклона с амплитудой а0 = 5° и периодом колебаний ^0 = 1.
Соответствующие результаты представлены на рис. 5. Пунктирная линия (корневой закон Ландау) соответствует стационарным решениям без модуляции (а = 0). При модуляциях периодическая смена направления вращения отчётливо реализуется до значений
числа Грасхофа Gr < 7 • 10 . Для чисел Грасхофа Gr > 8 • 103 в результате установления реализуется однонаправленный режим колебаний (без смены направления вращения жидкости). Соответствующий диапазон значений экстремума функции тока указан штриховкой. В области вблизи Ог « 7500 велико время выхода на установившиеся колебания. Эта область на рис. 5 отмечена отсутствием штриховки.
Рис. 5. Зависимость амплитуды функции тока от числа Грасхофа Отметим, что при значении периода колебаний t0 > 1 и значениях числа Грасхофа
Gr > Gr* = 8.6s • 103 max(ym (t)) близок к зна-
t
чениям функции тока стационарного решения на "благоприятной" ветви, а min(ym (t)) - к
t
значениям на "неблагоприятной" ветви. Этот факт свидетельствует о большом значении периода модуляции.
Перейдем к рассмотрению результатов для больших значений амплитуды модуляции а0. Пусть а0* - критическое значение амплитуды модуляции, разделяющее два различных режима установившихся колебаний - без смены направления вращения жидкости и со сменой. Эта граница зависит от числа Грас-хофа. На рис. б показано, что при а0 = бо < а0 да б.бо реализуется режим с однонаправленным вращением жидкости (ли*
ния 1), а при а0 = 70 > а0 - режим со сменой направления (линия 2).
Л —
ч о _ - .... 2 \
с. П — \
и о _ г :
-4- £ ■У
и 1
и ' 1 1 0.5 1 1 1 1.5 1 I 2.
Рис. 6. Зависимость амплитуды функции тока от времени
Модуляция включается после получе-
ния
стационарного
решения
для
Gr = 8.65 • 10 и а = 0 в момент ^ * 0.467 с фазой, обеспечивающей а(^) = 0.
Как видно, до времени t * 1 обе линии примерно совпадают. При t > 1 эти зависимости существенно различны: при а0 = 6° экстремальное значение функции тока остается отрицательным, а при а0 = 7° на интервале
времени 1.3 < t < 1.8 существуют положительные экстремальные значения функции тока, что свидетельствует о смене направления вращения жидкости.
Следущая серия ВЭ была направлена на выяснение зависимости а0 ^г). Обработка
результатов этих многочисленных ВЭ позволила получить приближенную формулу
(Ч = 1): ________
а0 * 0.112^/Gr — Gr1. (3 5)
На исследуемом интервале чисел Грасхофа от 6 • 103 до 15 -103 погрешность этой формулы менее 13%. В следующих экспериментах выяснялась зависимость а0* от периода модуляций.
Для периода модуляций t0 = 0.1, например, соответсвующая формула имеет вид а0* * 1.127 у!0г—0г[ . (3.6)
Формула дает значения угла с отклонением не выше 8% для чисел Грасхофа от 6 • 103 до 20 -103. Отметим, что отношение коэффициентов в формулах (3.5) и (3.6) примерно совпадает с отношением частот, равным 10.
Выполненные ВЭ показывают, что для реализации колебаний со сменой направления вращения жидкости требуется выполнение неравенства а0 > а**^г, t0). Логично предположить, что для реализации течений со сменой направления вращения жидкости требуется выполнение двух условий: *
Дат = (а0 — а ) > 0 и достаточного интервала
*
времени ЫС1, при котором а > а . Оба этих
условия могут быть выполнены при достаточных значениях интеграла
*2
|(а^) — а ^. (3 7)
В формуле (3.7) ^, t2 - время начала и
*
окончания превышения а(1) над а на одном периоде модуляций (см. рис. 7).
Рис. 7. Превышение а над а *
В следующих ВЭ определялся минимум
*
величины Ба , обеспечивающий колебательный
режим со сменой направления вращения жидкости. Перебор значений амплитуды модуляции а0 с шагом Да = 0.1° при фиксированном периоде модуляций t0 = 1 и значениях числа
Грасхофа от 8.65 • 103 до 14 • 103 позволил получить методом МНК линейную зависимость
*
5а с погрешностью, не превышающей 5%:
= 7.86 • 103 — 4.34 • 107 Gr.
(3.8)
Далее были исследованы конвективные режимы с параметром асимметрии (1.4), отличным от 1. При большом периоде модуляции t0 > 1 и параметре асимметрии % = 3 и % = 1/3 критическое значение угла а0* меняется незначительно.
Однако асимметрия сказывается на экстремальных значениях функции тока. При % = 3 минимальное значение функции тока -
4.94, максимальное - 2.52. При £ = 1/3 наоборот минимальное значение функции тока -2.52, максимальное - 4.94 (укажем, что при £ = 1 минимальное значение функции тока -4.60, максимальное - 4.60).
При уменьшении периода модуляции значение а0 уменьшается. Так, например, при t0 = 0.1 и £ = 3 коэффициент в формуле (3.6) уменьшается 2.87 раза.
Заключение
Найдены зависимости характеристик ветвления при постоянных значениях угла наклона от угла наклона а . Показано, что эти зависимости являются вспомогательными при определении эффектов медленных модуляций с периодом t0 > 1. Найдено критическое значение величины S* (3.8), при превышении
которой в установившемся решении реализуется режим модуляции со сменой направления. Отмечены особенности асимметричных модуляций.
Автор выражает благодарность
В.Г.Козлову и Т.П.Любимовой за полезные дискуссии и Е.Л.Тарунину за руководство.
Список литературы
1. Гершуни Г.З., Жуховицкий Е.М. Конвективная устойчивость несжимаемой жидкости. М.: Наука, 1972. С.242-255.
2. Тарунин Е.Л. Вычислительный эксперимент в задаче свободной конвекции: учеб. пособие. Иркутск: Изд-во. Иркут. ун-та, 1990. 228 с.
3. Gershuni G.Z., Lyubimov D.V. Thermal Vibrational Convection. Wiley, NY etc. 1998. 358 p.
4. Тарунин Е.Л. Обзор особенностей асимметричных колебаний // Проблемы механики и управления: Нелинейные динами-
ческие системы. Пермь: Изд-во Перм. унта, 2005. №37. С. 169-187.
5. Козлов В.Г. Вибрационная тепловая конвекция в полости, совершающей высокочастотные вращательные качания // Изв. АН СССР. Механика жидкости и газа. 1988. №2. С.138.
6. Детлаф А.А., Яворский Б.М. Курс физики: учеб. пособие для втузов. 4-е изд., испр. М.: Высш. шк., 2002. 718 с.
7. Ландау Л.Д., Лифшиц Е.М. Теоретическая физика: учеб. пособие. В 10 т. Т.1. Механика. 4-е изд., испр. М.: Наука, 1988. 216 с.
8. Ландау Л.Д., Лифшиц Е.М. Теоретическая физика: учеб. пособие. В 10 т. Т.У1. Гидродинамика. 3-е изд., перераб. М.: Наука, 1986. 736 с.
9. Самарский А.А. Теория разностных схем: учеб. пособие. М.: Наука, 1977. 656 с.
10. Келлер И.О., Тарунин Е.Л. Управление устойчивостью конвективного равновесия жидкости, подогреваемой снизу // Изв. АН СССР. Механика жидкости и газа. 1990. №4.
11. Чернатынский В.И., Шлиомис М.И. Конвекция вблизи критических чисел Рэлея при почти вертикальном градиенте температуры // Изв. АН СССР. Механика жидкости и газа. 1973. №1. С.64-70.
12. Тарунин Е.Л. Ветвление решений уравнений конвекции в замкнутой полости с подвижной границей при подогреве снизу // Современные проблемы тепловой гравитационной конвекции. Минск, 1974.
С.51-58.
13. Вайнберг А.А., Треногин В.А. Теория ветвления решений нелинейных уравнений. М.: Наука, 1969. 528 с.
14. Иванова А.А., Козлов В.Г. Вибрационная конвекция при непоступательных колебаниях полости (Изотермический случай) // Изв. РАН. Механика жидкости и газа. 2003. №2. С.25-32.
Asymmetric rotation oscillations effects in free convection task in a square-cross-section cavity with heating from below
A. E. Menenmbee
A. B. Melentyev
Perm State National Research University, Russia, 614990, Perm, Bukireva st., 15 a.b.melentyev@mail.ru; +7 906 889-89-49
Effects of heat convection with heating from below and presence of rotation oscillations were investigated with math modeling method. The task was heat convection in a square-cross-section cavity. Rotation oscillations were generated by rotation of the cavity around inner axis. Symmetric and asymmetric oscillations were considered. Threshold values of parameters were found, these values correspond change of a flow direction in case of stationary and oscillation modes. It was found that deviation from oscillations symmetry leads to threshold parameters values shift.
Key words: convection; rotation oscillations; asymmetry.