О возникновении конвекции в слое жидкости со свободной поверхностью. Слабонелинейный анализ Текст научной статьи по специальности «Физика»

ВЕСТНИК ПЕРМСКОГО УНИВЕРСИТЕТА

2011 Серия: Физика Вып. 1 (16)

УДК 532.5

О возникновении конвекции в слое жидкости со свободной поверхностью. Слабонелинейный анализ

А. Е. Самойлова

Пермский государственный университет, 614990, Пермь, ул. Букирева, 15

В работе проведен слабонелинейный анализ задачи о колебательной устойчивости плоского горизонтального слоя жидкости с деформируемой свободной поверхностью. При этом учет плавучести проводился за рамками приближения Буссинеска. Разложение по малому параметру проведено до третьего порядка включительно. Задача первого порядка совпадает с линейной задачей устойчивости. Её решение ищется в виде двух бегущих в разные стороны волн. В третьем порядке из условия разрешимости получена система уравнений Ландау. Определены стационарные решения этой системы, а также их устойчивость. В пространстве параметров задачи обнаружены области мягкого и жесткого возбуждения конвекции. Показано, что в системе течение возбуждается в виде бегущих волн или в виде стоячей волны.

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

1. Введение

Исследования конвективной неустойчивости горизонтального плоского слоя жидкости имеют давнюю историю. Началом её систематического изучения принято считать эксперименты Бенара и последующий теоретический анализ Релея (см. подробный обзор в работе [1]). Релею удалось аналитически определить порог возникновения монотонной тепловой конвекции в плоском горизонтальном слое жидкости для случая обеих свободных границ при подогреве снизу. Позже было показано существование колебательных возмущений при нагреве сверху. Их генерация обусловлена действием подъемной силы.

Как сказано выше, Релей считал свободные границы рассматриваемого горизонтального слоя плоскими и недеформируемыми. В действительности же свободная поверхность может деформироваться под действием возмущений. Анализ влияния деформируемости свободной границы на развитие конвекции Релея-Бенара приведен в [1]. Заметим, что этот анализ опирался на результаты работы [2], которая теперь представляется ошибочной.

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

верхностного натяжения. Исследования термокапиллярного механизма неустойчивости восходят к работе Пирсона [3].

Конвекцию Релея-Бенара-Марангони в рамках приближения Буссинеска изучали многие ученые. В работе [4] рассматривалась устойчивость относительно монотонных ячеистых возмущений. Задача решалась в рамках приближения Буссинеска, свободная поверхность предполагалась не деформируемой. Учет деформаций свободной поверхности был выполнен в работе [5]. Колебательная неустойчивость была обнаружена в [6], где показано (в рамках приближения Буссинеска), что она существует при отрицательных числах Марангони. Подробный обзор исследований неустойчивости Ма-рангони приведен в книге [7].

Как известно, корректное описание эффектов плавучести при деформируемой свободной поверхности возможно лишь ценой отказа от приближения Буссинеска. При этом переменность плотности должна быть учтена везде, а не только в слагаемом с подъемной силой. Ясно, что количественные характеристики неустойчивости могут оказаться очень чувствительными к виду уравнения состояния. В работе [8] предложена модель учета влияния плавучести на неустойчивость Релея-Бенара-Марангони, при этом использована линейная и экспоненциальная зависимость плотности от температуры. В рамках упомянутой модели рассмотрена монотонная длинноволновая и ячеи-

© Самойлова А. Е., 2011

стая неустойчивость Бенара-Марангони-Релея. Исследования выполнялись в предположении отсутствия теплоотдачи со свободной границы. Показано, что использование приближения Буссинеска в условиях пониженной гравитации приводит к неверным результатам.

Линейная теория устойчивости строится на предположении о том, что возмущения основного состояния малы. Эта теория позволяет определить границу устойчивости. В области неустойчивости малые возмущения со временем нарастают экспоненциально, становятся конечными и перестают описываться линейной теорией. Нелинейная теория устойчивости сопряжена с большими математическими трудностями. Достаточно строгую математическую основу удается подвести под исследования вторичных режимов в припороговой области. Согласно предложенной Л.Д.Ландау [9] идее исследования характера установившегося нестационарного движения в надкритической области, форма вторичных возмущений вблизи порога может описываться всего одной функцией - амплитудной. В слабонелинейном анализе используется также метод малого параметра, предложенный В.С.Сорокиным [10].

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

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

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

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

Рассмотрим поведение неоднородно нагретой жидкости в горизонтальном слое исходной толщины h со свободной деформируемой верхней границей и твёрдой нижней границей. Градиент температуры вертикален, слой находится в постоянном поле тяжести. Нижняя граница слоя идеально теплопроводная, ее температура T. На верхней границе происходит теплоотдача по закону Ньютона, при этом граница будет иметь некоторую температуру T (рис. 1).

Рис. 1. Геометрия задачи

Жидкость считается изотермически несжимаемой. Плавучесть учитывается за рамками приближения Буссинеска. В данной работе рассматриваются двумерные возмущения состояния равновесия, распространяющиеся в плоскости Х01 .

Полная краевая задача, описывающая поведение данной системы, в безразмерном виде такова:

р ( ду Л

—I-------+ V • Уу I = -V» + Ду +

Рг ^ дґ )

+1 У&уу - вару,

дТ

---Ь V •УТ = ДТ ,

дґ

Сілу = еДТ ,

2 = 0: у = 0, Т = 1,

дС .

2 = С : и = ^Т+, дґ

-Т С Т

хЬх г+ , 2 = ВіТ -Т) ,

Ф+с л/ї+ГТ

- p + 2

u z + v x - Cx (u x + v z )

1 + Cx2 '

-—divv =Гх^ (l + Г2 )-32,

x

—divv =

3 Cr

2 (uz - v x )Гх +(ux + vz )(1 -Cx2 ) =

2\V2

= -Ma (Tx + Tzrx )(l + Cx2 )1 .

(1)

Здесь y - орт оси z ; Q - описывает отклонение свободной поверхности от состояния равновесия; v = (v, 0, u); - температура вне слоя на беско-

нечности; нижние индексы означают взятие частной производной.

Безразмерные параметры задачи: число Пран-дтля Pr = v/x (при этом v=^/p(T0) , p(T0) -

плотность на верхней границе в состоянии механического равновесия, ц - коэффициент динамической вязкости), параметр Буссинеска s = р& ,

число Галилея ва = р(Т0)ghъ /щ , параметр капиллярности Сг = г)х1^, число Марангони Ма = ат©h/щх , число Био Ы = mh|к . Здесь © = Т - Т - характерный масштаб температуры, А = -р- (Т0) ёр/ёТ - коэффициент объемного расширения, х - коэффициент температуропроводности, ат = -да/дТ - температурный коэффициент поверхностного натяжения, т - коэффициент теплоотдачи с верхней границы, к -коэффициент теплопроводности.

Уравнение состояния принято в экспоненциальной форме: р = е~еТ .

Уравнения в задаче (1) допускают решение, отвечающее состоянию механического равновесия:

Т0 = 1 - ^ Ро = -1), дРо/дг = -Рова. (2)

3. Слабонелинейный анализ

Для исследования устойчивости решения (2) по отношению к конечным возмущениям с произвольной длиной волны вблизи порога возникновения конвекции применяется метод многих масштабов [12]. Разложим по малому параметру поля скорости, давления, температуры, плотности, отклонение от плоской поверхности, а также производную по времени:

V = дг1 + 5^2 + 5^3 + К ,

Р = Ро + 5Р, + 5Р2 + 5Рз +К ,

Т = Т + 5% + 5% + 5% + К , р = р0 +5р1 +§2р2 + 5р + К, (3)

С=5С! +52С2 +53СзК,

д д д 2 д

— = — + 5 — + 8 — +... д, д,0 д, д,2

Управляющим параметром обнаруженной ранее дополнительной моды колебательной неустойчивости является параметр Буссинеска. В данной работе исследуется поведение системы в припоро-говой области: е=е0+52£2. Малый параметр 5 здесь имеет смысл корня из надкритичности.

При подстановке разложения (3) в задачу (1) в первом порядке по малому параметру получаем линейную краевую задачу. Её решение ищется в виде двух бегущих в разные стороны волн:

т (х, г, г) = А (,, ,2,...) е^-4 т(1) (г) +

*2

+В (,, ,2,...) е-кх-и° Т(2) (г) + к.е.

(4)

Здесь Т обозначены поля скорости, давления и температуры, а также отклонение от плоской поверхности. Сокращение “к.с.” означает комплексно-сопряженные слагаемые.

Из условия разрешимости нелинейной задачи второго порядка следует равенство нулю производной по первому медленному времени: д/ д, = 0 .

В третьем порядке из условия разрешимости нелинейной задачи получаем систему уравнений Ландау:

— = (лС - Ко |А|2 - К, |В|2) А,

2

ёВ ^ - Ко \В\2 - К, |А|2) В.

(5)

Ж,

Коэффициент /л и коэффициенты Ландау Ко и К выражаются через собственные функции первого и второго порядка:

| (% % - кX ) + роОайТ (1 - е%)) ёг

Л= -

I +£ %и,+тт +£о%? т | &

Рг

I (% ао + %Фо + %о + % Ф) ёг

Ко = Т

(6)

+;р° %и1+%Т1 +£о%? Т1 ]

а1 + %рА + % + % ф, )^

к, = т

1(р°%^ + р°%и, + ТТ + С>%Т ]ёг

Рг

где значком : обозначены собственные функции

сопряженной линейной задачи, штрих означает взятие производной по координате г , а

а, А, X,-, Ф - достаточно громоздкие функции от собственных значений и функций задач первого и второго порядка.

При анализе системы (5) обнаружено 3 равновесных решения, помимо тривиального (когда

II? II?

А = В = о):

1) |а|2 = о, |В|2 =

Яе (л) С 2 Яе (Ко)

> о;

2) 'В'2 = ^

Яе (/“)С2

3) |А|2 = |В|2 =

Яе (Ко + К,)

> о;

> о.

(7)

Здесь Яе означает взятие действительной части числа.

Из вида решения линейной задачи (4) следует, что случаи 1) и 2) отвечают бегущим волнам, а случай 3) - стоячей волне.

Анализ устойчивости решений (7) выявил, что в зависимости от значений констант л , К0 и К реализуются различные ситуации.

Рассмотрим случай, когда

Яе (К0) > о, Яе (К0 + К) > о. Если при этом

Яе(л) > о , то на фазовой плоскости |А|2 - |В|2 существуют четыре положения равновесия: вышеперечисленные три и нулевое, т.е. мы находимся в надкритической области. При этом нулевое положение равновесия оказывается равновесием типа “неустойчивый узел”, а три других - типа “седло” или “устойчивый узел” (в зависимости от соотношения между коэффициентами Ландау). Если Яе (л) < о , то система (5) имеет только тривиальное решение, которое устойчиво относительно малых возмущений (устойчивый узел); мы находимся в подкритической области.

В случае Яе (Ко) < о, Яе (Ко + К,) < о при

Яе (л) > о существует только нулевое положение

равновесия (подкритическая область), малые возмущения которого нарастают (неустойчивый узел). В случае Яе ( К ) < о, Яе ( К + К ) < о при Яе (л) < о существуют ещё три, помимо нулевого,

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

Из вышесказанного следует, что случай Яе (Ко )< о, Яе ( К + К ) < о отвечает жесткому возбуждению конвекции, а случай

Яе (Ко )> о, Яе (Ко + к, ) > о - мягкому возбуждению. Для последнего случая имеет смысл исследовать, в какой из форм - бегущей или стоячей волны - возбуждается конвекция. Из анализа устойчивости положений равновесия (7) следует, что при Яе ( К ) < Яе ( К ) одновременно устойчивы равновесия 1) и 2), вторичное возмущение формируется в виде бегущей волны. При Яе (Ко )> Яе (К) устойчивым оказывается равновесие 3), возбуждается стоячая волна.

4. Методы решения

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

Для решения линейной задачи используется метод построения фундаментальной системы решений [13]. Для сохранения линейной независимости частных решений на каждом шаге интегриро-

вания применялась процедура ортогонализации [14-16]. Использовались два варианта процедуры ортогонализации. При нахождении собственных значений линейной задачи устойчивости применялась ортогонализация по Грама-Шмидту. При определении собственных функций спектральной краевой задачи устойчивости использовался следующий вариант ортогонализации. Тройка векторов частных решений поворачивается так, чтобы последние три компоненты этих векторов являлись элементами единичной матрицы. При таких манипуляциях с векторами частных решений собственные значения задачи не меняются, но собственные функции претерпевают существенные изменения. Поэтому перед нахождением общего решения следует восстановить векторы частных решений в каждой точке ортогонализации.

При интегрировании систем обыкновенных дифференциальных уравнений применяется метод Рунге-Кутты-Мерсона. Модификация Мерсона (метод пятого порядка точности) позволяет оценить погрешность приращения, что в свою очередь позволяет гибко менять величину шага интегрирования, получая решение уравнения с требуемой точностью за приемлемое время.

Для решения систем нелинейных уравнений используется метод двумерных секущих.

Основные результаты, представленные ниже, получены при значении числа Био В1 = о .

5. Основные результаты

На рис. 2 представлены графики зависимости Яе (К) и Яе (К + К) от числа Прандтля (сплошная и пунктирная линии соответственно).

Рис. 2. Коэффициенты Ландау в зависимости от числа Прандтля (е = -в.1, ва = о, Сг = 2-Ю-6)

Графики на данном рисунке получены для области параметров, в которой ранее автором была обнаружена дополнительная мода колебательной

неустойчивости. Для наглядности на рис. 3 приведем зависимость критических значений числа Ма-рангони и волнового числа от числа Прандтля.

12

Рис. 3. Зависимость экстремальных волновых чисел (пунктирная линия) и значений числа Марангони (сплошная линия) от числа Прандтля (є = -0.1, ва = 0,

Сг = 2 -10-6)

Как видно из рис. 2, существует область мягкого возбуждения обнаруженной дополнительной моды колебательной неустойчивости (в интервале значений числа Прандтля от 0.0002 до 0.005).

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

Рис. 4. Зависимость действительной части разности коэффициентов Ландау от числа Прандтля (є =-0.1, ва = 0,

Сг = 2 -10-6)

Из данного рисунка видно, что в том интервале значений числа Прандтля, где конвекция возбуждается мягко, возмущение формируется в виде бегущей волны.

6. Заключение

Проведен слабонелинейный анализ задачи о колебательной устойчивости плоского горизонтального слоя жидкости с деформируемой свободной поверхностью. Разложение по малому параметру проведено до третьего порядка включительно. Задача первого порядка совпадает с линейной задачей устойчивости. Её решение ищется в виде двух бегущих в разные стороны волн. В третьем порядке из условия разрешимости получена система уравнений Ландау. Определены стационарные решения этой системы, а также их устойчивость. Определены области мягкого и жесткого возбуждения неустойчивости для значений параметров задачи s=-0.1, Ga = 0, Cr = 2 -10-6 (в этой области параметров существует дополнительная мода колебательной неустойчивости). Показано, что в области мягкого возбуждения конвекция возникает в виде бегущей волны.

Автор выражает глубокую признательность проф. Д.В.Любимову и к.ф.-м.н. С.В. Шкляеву за неоценимую помощь и постоянное внимание к исследованию.

Работа выполнена при частичной финансовой поддержке из средств гранта RUX0-009-PE-06 Американского Фонда Г ражданских Исследований и Развития (АФГИР).

Список литературы

1. Гершуни Г. З., Жуховицкий Е. М. Конвективная неустойчивость несжимаемой жидкости. М.: Наука, 1972. 392 с.

2. Изаксон В. Х., Юдович В. И. О возникновении конвекции в слое жидкости со свободной границей // Прикл. мех. и техн. физ. 1969. №3. С. 89-92.

3. Pearson J. K. A. On convection cells induced by surface tension // J. Fluid Mech. 1958. Vol. 4, N 5. P. 489-495.

4. Nield D. A. Surface tension and buoyancy effects in cellular convection // Ibid. 1964. Vol. 19. P. 341-352.

5. Scriven L. E., Sterling C. V. On cellular convection driven by surface-tension gradients: Effects of mean surface tension and surface viscosity // Ibid. 1964. Vol. 19. P. 321-332.

6. Takashima M. Surface-tension driven instability in a horizontal liquid layer with a deformable free surface. II. Overstability // J. Phys. Soc. Jаp. 1981. Vol. 50. P. 2751-2756

7. BirikhR. V., Briskman V.A., Velarde M. G, Legros J.-C. Liquid Interfacial Systems: Oscillations and Instability. New York-Basel: Marcel Dekker Inc., 2003. 367 р.

8. Lyubimov D. V., Lyubimova T. P., Alexander Iwan J. D. and Lobov N. I. On the Boussinesq approxi-

в

4

mation for fluid systems with deformable interfaces // Adv. Space Res. 1998. Vol. 22, N 8. P. 11591168.

9. Ландау Л. Д., Лифшиц Е. М. Теоретическая физика. Т. 6. Гидродинамика. M.: Наука, 1986. 736 с.

10. Сорокин В. С. О стационарных движениях жидкости, подогреваемой снизу //Прикл. мех. и матем. 1954. Т.18, №2. С. 197.

11. Лобов Н. И., Самойлова А. Е. Колебательная устойчивость плоского слоя жидкости со свободной деформируемой поверхностью // Конвективные течения...: сб. науч. тр. / Перм. гос. пед. ун-т. Пермь, 2009. Вып.4. С. 35-50.

12. Найфе А. Х. Методы возмущений. М.: Мир, 1976. 455 с.

13. Лобов Н. И., Любимов Д. В., Любимова Т. П. Численные методы решения задач теории гид-

родинамической устойчивости: учеб. пособие // Перм. ун-т. Пермь, 2004. 101 с.

14. Бирих Р. В., Рудаков Р. Н. Применение метода ортогонализации в пошаговом интегрировании при исследовании устойчивости конвективных течений. Часть I // Гидродинамика: межвуз. сб. науч. тр. / Перм. ун-т. Пермь, 1974. Вып. V. С. 149-158.

15. Бирих Р. В., Рудаков Р. Н.. Семакин И. Г. Применение метода ортогонализации в пошаговом интегрировании при исследовании устойчивости конвективных течений. Часть II // Конвективные течения: сб. науч. тр. / Перм. гос. пед. ун-т. Пермь, 1979. С. 58-6о.

16. Корн Г., Корн Т. Справочник по математике для научных работников и инженеров. М.: Наука, 1970. 720 с.

On the onset of the convection in the liquid layer with a deformable surface. Weakly nonlinear analysis

A. E. Samoilova

Perm State University, Bukirev st., 15, 614990, Perm

In this paper the weakly nonlinear analysis of oscillatory instability of the flat horizontal liquid layer with a deformable free surface is performed. The problem is solved within the correct account buoyancy model, where the buoyancy is taken into account in every density term. The small-parameter expansion is carried out till the third order. The first order problem coincides with the linear instability problem. We accept two traveling waves-form solution for the first order problem. System of Landau equations is obtained from the solvability condition of third order problem. The stationary states of this system and their stability are determined. Parameter region of soft and hard convection excitation is revealed. It is showed, that the convection onsets in the form of the traveling wave or the standing wave.

Keywords: flat liquid layer, deformable surface, weakly nonlinear analysis.