СЛАБОНЕЛИНЕЙНЫЙ АНАЛИЗ УСТОЙЧИВОСТИ ТЕРМОАКУСТИЧЕСКОГО АДВЕКТИВНОГО ТЕЧЕНИЯ Текст научной статьи по специальности «Физика»

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

2019

• ФИЗИКА •

Вып. 3

УДК 532.5 PACS 43.25.Nm

Слабонелинейный анализ устойчивости термоакустического адвективного течения

А. О. Иванцов^

a Пермский государственный национальный исследовательский университет 614990, Пермь, ул. Букирева, 15 email: ivants2@psu.ru

b Институт механики сплошных сред УрО РАН, 614013, Пермь, ул. Ак. Королева, 1

Изучается влияние бегущей звуковой волны на устойчивость адвективного течения по отношению к спиральным возмущениям. Проведены слабонелинейный анализ устойчивости и прямое численное моделирование. Из условия разрешимости задачи третьего порядка получено дифференциальное уравнение, описывающее эволюцию амплитудных функций. При этом для вычисления постоянных, входящих в данное уравнение, необходимо определить собственные функции линейной задачи устойчивости и решения некоторых неоднородных линейных задач. Использовался метод стрельбы с применением процедуры ортогонализа-ции. Показано, что надкритические режимы возникают мягко и устойчивы вблизи порога относительно модуляции. Прямое численное моделирование задачи осуществлено с помощью метода конечных разностей. Проведено исследование конвективного течения при значительном удалении от порога устойчивости. Установлено, что при определенном значении параметров в слое возникают колебания сложной структуры. Ключевые слова: адвективное течение; слабонелинейный анализ; метод амплитудных функций

Поступила в редакцию 03.06.2019; принята к опубликованию 16.09.2019

Weakly non-linear analysis of thermoacoustic advective flow

A. O. Ivantsova,b

a Perm State University, Bukireva St. 15, 614990, Perm email: ivants2@psu.ru

b Institute of Continuous Media Mechanics, 614013, Perm

The influence of a traveling sound wave on the stability of an advective flow with respect to spiral perturbations is studied. A weak nonlinear stability analysis and direct numerical simulation were performed. From the condition of solvability of a third-order problem, a differential equation is obtained, which describes the evolution of amplitude functions. At the same time, to calculate the constants in this equation, it is necessary to determine the eigenfunctions of the linear problem of stability and the solution of some inhomogeneous linear problems. The shooting method was used using the orthogonalization procedure. It is shown that supercritical regimes appear soft and stable near the threshold with respect to modulation. Direct numerical simulation of the problem was carried out using the finite difference method. A study of convective flow with a significant distance from the stability threshold was carried out. It has been established that for a certain value of parameters in the layer, oscillations of a complex structure occur. Keywords: advetive flow; weak nonlinear analysis; amplitude function method

Received 03.06.2019; accepted 16.09.2019

doi: 10.17072/1994-3598-2019-3-28-44

© Иванцов А. О., 2019

распространяется на условиях лицензии

Creative Commons Attribution 4.0 International (CC BY 4.0).

1. Введение

Интерес к изучению адвективных течений связан с рядом геофизических и технологических приложений. К ним относятся, в частности, атмосферная циркуляция Хэдли, некоторые типы движения в океане, коре и мантии Земли, движение расплава в установках для получения кристаллов. Для многих технических приложений особенно важна задача управления устойчивостью адвективных течений, поскольку потеря устойчивости обычно приводит к резкому изменению характера течения.

Линейная задача устойчивости адвективного течения рассматривалась многими авторами, подробный анализ работ можно найти в [1, 2]. В частности, известно, что при больших и умеренных числах Прандтля важную роль играют монотонные трехмерные возмущения типа продольных валов (спиральная мода). Конвективные валы локализованы в пристеночных областях, где жидкость стратифицирована неустойчиво.

Влияние бегущей звуковой волны на устойчивость адвективного течения изучалось в работе [3]. Получен профиль комбинированного течения и изучена его устойчивость; показано, что в зависимости от параметров задачи наиболее опасными являются бегущие плоские возмущения, стационарные продольные валы и наклонные волны. Проведен также слабонелинейный анализ задачи, получено дифференциальное уравнение, описывающее эволюцию амплитудных функций. Позднее, при численном исследовании этого амплитудного уравнения, было обнаружено, что для спиральной монотонной моды возможно жесткое возбуждение конвекции [4]. Этим фактом и вызван интерес к более подробному исследованию вторичного течения в виде продольных валов.

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

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

Рассмотрим бесконечный горизонтальный слой жидкости толщиной 2Ъ, заключенный между твердыми границами, на которых поддерживается линейное распределение температуры: Т (х, у, ±к ) = Ах (плоскость х, у расположена на

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

и = асо/сое ®(/-х/сзв),

(1)

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

жидкости, / - орт оси х. Длина звуковой волны

предполагается сравнимой с толщиной слоя.

Для описания возникающих в слое конвективных течений применимы осредненные уравнения термоакустической конвекции, полученные в работе [5]. Вследствие однородности звукового поля уравнения термоакустической конвекции в данном случае совпадают с обычными уравнениями тепловой конвекции в приближении Буссинеска:

— + (v • V) v = —- Vр + + я ВТ у, д1 V ' р

дТ_

~дг

+ (у.у)Т = хАТ,

(2)

С1гуу = 0.

Здесь у, Т - средняя скорость и температура жидкости; р - конвективная добавка к гидростатическому распределению давления; р - средняя плотность; V , % - коэффициенты кинематической вязкости и температуропроводности, предполагаемые постоянными; g - ускорение свободного падения; р -коэффициент теплового расширения; у - единичный вектор, направленный вертикально вверх.

При распространении звуковой волны вдоль твердой границы в динамическом пограничном слое генерируется осредненное течение [1, 5], которое можно учесть с помощью эффективных граничных условий. В данном случае вследствие однородности акустического поля граничное условие запишется в виде

г = ±к:

3 а2а>2

4 с

(3)

Остальные компоненты скорости обращаются в нуль на твердых стенках, температура границ линейно зависит от продольной координаты х :

г =±Н : у = V = 0, Т = Ах.

(4)

В ходе решения задачи предполагается, что в среднем расход жидкости через плоскость, перпендикулярную оси у, г, отсутствует:

Q = | | V,¿гау = 0.

(5)

Перейдем к безразмерным переменным, выбрав в качестве единиц измерения длины, времени, скорости, давления, температуры следующие величины: к, к2¡V , Vк, рv2/к2, Ак . Тогда задача (2)-(4) принимает вид

т-1

^ + (■v • v) v = -Ур + ау + вгТу, дТ 1

— + (у-У)Г = — АТ, (6)

д1 V ' Рг

сЦуу = 0,

2 = ±н: V = Яе, V = V = 0, Т = х, (7)

где Ог = gPAh'í/V2 - число Грасгофа, Рг = у/% -число Прандтля и Яе = 3а2а>2к/^сзв - число Рей-нольдса, определенное по скорости акустического течения. Оценим возможные значения числа Рей-нольдса. Пусть частота звуковой волны, распространяющейся в слое, равна 106 Гц, амплитуда вибраций - 10-5 см, толщина слоя - 10 см, а скорость звука в жидкости примем равной 1.5 -105 см/с. Тогда Яе будет порядка 50 при вязкости жидкости 0.01 Ст.

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

Система уравнений и граничных условий (6)-(7) допускает решение, описывающее плоскопараллельное конвективное течение. Поля скорости

и температуры для такого течения имеют вид

Ог з Яе,,, 2 1ч

щ =— (г3 -2) + — (322 -1), 6 2

ОГРГ

Т = X + $п = х + -^(325 - 1023 + 72) +

360

(8)

+Яе*Т( 2 2 - 1)2. 8

Основное течение представляет собой суперпозицию двух течений - классического адвективного [6] и течения с квадратичным профилем. Профили скорости и добавки к линейному распределению температуры (2) для разных чисел Рейнольдса представлены на рис. 1. Как видно, акустическое воздействие вызывает существенное изменение структуры плоскопараллельного течения. Устойчивость решения (8) в линейной постановке изучена в работе [3]. В частности, показано, что в случае достаточно интенсивного звукового воздействия наиболее опасными становятся спиральные возмущения (продольные валы). Настоящая работа посвящена исследованию нелинейных конвективных течений этого типа. Будем предполагать, что возникающее вторичное течение имеет вид про-

Рис. 1. Профили скорости (а) и добавки к температуре (Ь) для плоскопараллельного течения при Gr = 650: 1) Re = 0, 2) Re = 50, 3) Re = 100

СГ&Л

Рис. 2. Возмущения типа спиральных валов

x

дольных конвективных валов (рис. 2).

В этом случае поле скорости зависит лишь от двух поперечных по отношению к основному течению координат (у, г):

dv*

dx

dvv dv„

_v

dx

dx

= 0.

(9)

В распределении температуры можно выделить линейное по x слагаемое:

T = x + &( у, z).

(10)

Покажем, что в данном случае давление также является линейной функцией х . Для этого спроектируем уравнение движения (6) на оси координат:

dvx dvx dvx dp

— + vv—- + V— ~^vx =——, dt dy dz dx

dvv dvv dvv Л dp

—- + v —- + v— -Av =——,

dt 8y dz v 8y

dv dv dvz —- + v—- + vz—- -Avz - Gr5 = dt dv dz

= -dP + Grx. dz

(11) (12)

(13)

f = X (v, z),

dx

-8t = ¥ (V, z),

dv

(14)

-8p + Grx = Z (у, z).

dz

Тогда, интегрируя по z последнее соотношение, получаем

p = J Z(у, z)dz + Grxz + f (x, у): = p(y, z) + Grxz + f(x, y).

(15)

Подставляя (15) во второе соотношение (14), будем иметь

р(у, г) + = -У (у, г).

оу ду

Из полученного уравнения видно, что второе слагаемое является функцией, не зависящей от координаты х :

d

f (x,y) = f(y).

ду

Следовательно,

f (*, у) = J Ry)dy + G(x) = F(y) + G(x).

Тогда (15) можно переписать так: р = р{у, z) + F(y) + Grxz + G(x) = = p (у, z) + Grxz + G(x).

(16)

Наконец, подставляя (16) в первое из условий (14), получим

G'(x) = X(у, z) - Grz .

Левая часть уравнения зависит только от x, в то время как правая часть зависит от координат y, z. Следовательно, обе части уравнения должны быть равны одной и той же постоянной

G'( x) = const = c.

Таким образом, искомая производная имеет

вид

dp „ — = Grz + c, dx

(17)

где с - некоторая константа (среднее по слою значение продольного градиента давления). Поле давления запишется в виде

р = Grxz + сх + р{у, z).

(18)

Левые части уравнений (11)—(13) являются функциями, не зависящими от координаты х :

Константа с определяется из условия замкнутости потока.

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

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

3.1. Метод амплитудных функций

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

. С этим обстоятельством связана относительно простая форма вторичных течений, которые удается описать с помощью медленно меняющихся в пространстве и во времени функций, называемых амплитудными функциями [2].

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

и давления основного состояния и разложим производные по горизонтальным координатам и времени в ряды по формальному малому параметру е :

V = м0 + £v(1) + £2 v<2) + ...,

3 = 30+£3(1) +£23(2) +..., р=р<0) +£р(1) +£2р<2) + ...

д д д

— =-+ £— + ...,

дх дх0

д д д

-=-+ £-+ ...

ду ду0 дт]

д д д — =-+ £-+ ...

dt ôt0 дт

(19)

(20)

(21)

Число Грасгофа также представим в виде разложения

Gr = Gr0 + ¿'Cîг, +£ Gr2 +...

(22)

Подставляя разложения (19)-(21) в исходную краевую задачу, получим в первом порядке линейную задачу устойчивости плоскопараллельного течения (8) относительно пространственных возмущений (интересуясь слабонелинейным поведением стационарных возмущений в виде продольных валов, положим д/дх0 =д/д^ = 0):

(1) „.(1) du0 _

Au ' - w

= 0,

dz

Av(i)-^ = 0,

fyo

Aw«-^-Gr 3(1) = 0,

dz

(23)

A3«- Pr ^ u(1)- w«^ J = 0,

dv(1) dw(1)

-+-= 0.

dy0 dz

Здесь u , v, w - компоненты вектора скорости v. Возмущения скорости и температуры на твердых границах равны нулю, следовательно,

z = +1: ы(1) = v(1) = w(1)=3(1) = 0. (24)

Решение данной однородной задачи может быть представлено в виде

v'1' = as j (z) exp (iky0 ) + к.с., (25)

3(1) = a3( z) exp (iky0 ) + к.с., (26)

p(1) = ap (z) exp (iky0 ) + к.с. + b, (27)

где символом к.с. обозначены комплексно-сопряженные слагаемые, у1,31,р1 - собственные функции линейной задачи устойчивости. Амплитудная функция а = д (¿\ //. г...) характеризует интенсивность вторичного течения; Ь(^,т],т...) -

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

Для рх получаем линейную однородную

краевую задачу, изученную в [3]:

Du - ЩЫg = 0,

Dv - ikpx = 0, Dw - p[- Gr3 = 0,

Pr (u - w3g ) = 0, ikv + w[ = 0.

z = ±1 : ы = v = w = 3 = 0.

(28)

(29)

Здесь и далее штрихом обозначено дифференцирование по 2 . В системе (28) введен линейный дифференциальный оператор

D = 4 - k2. dz2

(30)

Собственные числа задачи (28)-(29) определяют порог устойчивости плоскопараллельного течения. Амплитуды движения а и Ь в первом порядке не определяются.

В ходе решения задач во втором и третьем порядке по е получим неоднородные системы дифференциальных уравнений, левые части которых совпадают с левыми частями уравнений задачи (28):

Дм(2) - ^^ = /, д2

дУ0

д (2)

Дм^^- -Г,З(2) = /з, (31)

д2

ДЗ2 -Рг[и<2) -м^3^ = /,

ду(2) дм(2) , дУ0 д2

Условие разрешимости таких неоднородных систем имеет вид

j и + Vf2 + Wf + ©f + Pf ) dz.

(32)

Здесь и, V, Ж, 0, Р являются собственными функциями задачи, сопряженной с (28)-(29):

Би - Ог0© = 0, DV + 1кР = 0,

БЖ - Р'- Ог0т'0- ыои = 0, Б0 + Рг Ж = 0, ikv - Ж' = 0.

г = ±1: и = V = Ж = 0 = 0.

(33)

(34)

Второе приближение по е приводит к задаче для нахождения у<2), 3'1', р'1', неоднородные части которой имеют вид

/1 =

ди(1) др(1) ди(1)

+ ип

дт д' 0 д' {(1]

+ у->— + щ(1) —--2

(1) ди(1) (1) ди(1) „ д2и(1)

Л >--и -и/ > -

ду0

дг ду0 дц'

/2 =

дУ(1) др(1) дУ(1)

+и

2 дт дц 0 д'

(1) ду(1) (1) ду(1) „ д2у(1) + У(1)-+ щ(>--2-

/з =

ду0

дщ(1) дщ(1)

дг ду0 дц

дт

- + ип

д'

+ V

,(1)

сы(1) (1) дм>() „ 52щ(1)

(35)

- + щ

дг

- 2-

/4 = РГ

дЗ(1) дЗ(1)

--+ и„--+

дт

д'

+ V

(1)

дЗ

(1)

д^0

- + щ

(1)

дЗ

(1)|

дг

,2о(1)

- 2

д 2З

5у0 дц'

/ =

ди(1) -+-

д' дц

Из условия разрешимости системы второго порядка следует равенство нулю линейной поправки вг . Таким образом, с точностью до е2 разложение (22) имеет вид

вг = вг + е вг2.

(36)

и(2) ={а\2еъку° + ( да да

{д' ' дц

+ 1 — и. +—и |е+ к.с. К

(37)

I |2 1 дЬ 2 1Ч + а и20 +---(г -1)

2 д'

={а\2е2Лу0 + (да да

{д' дц

+ 1 — V. + — V |ек° + к.сЛ +

(38)

1 дЬ

+ а ^0 + (г -1) 2 дц

щ(2)={а2 щгеъку0 +

(да да | кку \ , ,2 + 1— щ + — Щ Iе + к.с.\ + а щп,

[д' ' дц ц 1 ' 11 20

(39)

+(даз+ьЗ |¿у + к.си

З(2)= |а|2 З20+{а2З22е2ку0 3

{д'Г дц 11 Рг дЬ 24 д'

(40)

Рг дЬ 4 2 +---(г - 6г + 5),

р{2) = \а\2 р20+р(2) + {а2р22е2^ +

(да да

+{1Г'р'+Тцрц\е""0 +к£.

а

(41)

Задача второго порядка содержит как линейные, так и нелинейные неоднородности. Линейные неоднородности приводят к появлению в решении слагаемых, пропорциональных ек0, нелинейные -к слагаемым, не зависящим от у0, и к слагаемым,

пропорциональным е2йу°. Решение неоднородной задачи второго порядка может быть записано в виде

Таким образом, каждая из компонент поля скорости и поле температуры зависят от четырех функций координаты г . Подставляя решение (37)-(41) в систему уравнений второго порядка, получим четыре неоднородные задачи для нахождения данных функций.

Для функций и, , щ ,З , р^ будем иметь

Би,, - щи'0 = (Л + и)и + р, -ikp( = (Л + и0)V,

Бщ - р'(- вг0З = (Л + и0) щ, (42)

БЗ(-Рг-щЗ') = Рг(Л + и)З,

ikví + = -и,

г = ±1: и = V = Щ =З= 0. (43)

Здесь коэффициент Л определяется из условия разрешимости этой задачи и вычисляется по формуле

А = -=±-

Ди^ + ри - щР) ¿2

1

I Бй2

(44)

5 = щи + уу + му + -Г, З©.

Функции ит, ут, мт,Зт ,рт находятся из следующей задачи:

Би - м и'0 = -21кщ, Бу^ -1кр = -Икух + р, БМ]- р]- Ог0Зч=-2гкм1, БЗ]- Рг (и]- М]З'0) = -21кЗ1, 'кут + wT=-уl,

2 = ±1: и = у = м = З = 0.

— т т т т

(45)

(46)

Условие разрешимости задачи третьего порядка приводит к дифференциальным уравнениям для амплитудных функций:

да

= а%+ gaтт+ 1$ат +

+ а (1 - К\а\2-кЬ%- 1хь]), (52)

д а

■ = Ь„ + Ь .

й тт

Здесь индексами обозначаются производные по соответствующим координатам, g, 8, к, х и K - вещественные постоянные, определяемые численно по следующим формулам:

а а Ш

Условие разрешимости данной задачи представляет собой условие минимума нейтральной кривой Gr(k) и монотонности возмущений.

Для функций и22, у22 , щ2,З22 , Ри получим следующую неоднородную задачу:

Би22 -м>22и'0 = Иш1у1 Бу22 - Икр22 = 1к\>\ + у[ м>1, Бм>22-р'22-Сх0322 =0, (47)

Б322 - РГ(и22 -м>223'0) = РГ(//ЦЗ + 1^3'), Ику22 + м22 = 0,

2 = ±1: и22 = у22 = -м22 = З22 = 0, (48)

где

~ У

Б=--4к2. (49)

аг

Функции и20, у20, Мо, З , Рм находятся из следующей системы дифференциальных уравнений:

и"о - Мои2 = 2Ш (мМЩ ) ,

у2о = 0,

М20 - Р2о - -г0З20 = 2 (|^ ) ,

З20 - Рг (и20 - М2оЗ0) = 2Рг м>З )',

(50)

м2О = 0,

2 = ±1 : и20 = У20 = М20 = З20 = 0 ,

(51)

¿О, хО,

\к\ \к\

где

1 'г

5 =(иеи + у у + му + -г з©),

рг

§ = и1и+уу + м>11Г + —?- 3©, Рг

I'

-Б + Ру-УлР)с1г, Ят=(ии + уу + МтЖ),

(53)

(54)

(55)

(56)

. 1

$ = ^1{(Л + и0)(8ч+Сг0Зч®) + рчи-ичР-

-т(^ +З©-г/Рг)+ру-ур)¿2, (57)

г-Ц

1 (2 -2)ми -ЗУ +

ОгоРг 360

(325 -1023 + 72)м©

¿2, (58)

где ЭТ - операция взятия реальной части комплексного числа, звездочкой обозначена операция комплексного сопряжения.

к

+ 0^23 - 32)© ¿2, (59)

-1

2

-1

•4

I

11 ,/ог

(а) (Ь)

Рис. 3. Бифуркации решений: а) - прямая; б) - обратная. Пунктирной линией показаны неустойчивые ветви решений

К=-

^ 1

-1((и'22 +и'20)и) +у'22г + м>'22)¥ +

ьОг0 ((^ )©)+

+ щ22 (и[и - у[¥ + + вг03©) - (60)

- 2гку1 (и22П + у22У + + Сг0 322©) -(щи - уу + щЖ + вг03©)}Ж,

б = | и 20

(61)

I |2 1 а =—. 1 1 К

(62)

лов. Показано, что вблизи минимума нейтральной кривой периодическое вторичное течение устойчиво, если

+ +1 + я + кя + дх> 0 .

(64)

Характеризующий нелинейность коэффициент К называется константой Ландау. Характер ветвления решений вблизи порога зависит от знака К. Если знак константы Ландау положителен, то в критической точке происходит прямая вилочная бифуркация (рис. 3, а). Следовательно, возбуждение неустойчивости будет мягким. В случае отрицательного знака коэффициента К бифуркация решений будет обратной (рис. 3, б). Значит, возможно подкритическое развитие неустойчивости.

Для стационарных режимов квадрат модуля амплитуды определяется по формуле

Задача (52) аналитически изучена в работе [3]. Показано, что продольные валы возникают мягко при выполнении условия

1 + к> 0. (63)

Коэффициент я в уравнении (52) положителен, так как разложение (19)—(21) проводится вблизи минимума нейтральной кривой Ог0 (к).

Кроме того, в [3] исследована устойчивость регулярных вторичных течений в виде продольных ва-

Отметим, что в [3] формула (63) приводится дважды: как условие мягкого возбуждения (второе из условий (2.15) в цитируемой работе) и как одно из условий устойчивости вторичного течения (формула, следующая за формулой (2.18), там же). При этом в расчетах нарушение условия (63) ошибочно интерпретировалось как условие неустойчивости регулярных валов, а не как условие жесткого возбуждения. Эта неправильная трактовка условия (63) была отмечена в [4].

Чтобы вычислить коэффициенты, входящие в амплитудные уравнения (52), необходимо найти собственные значения и собственные функции линейной однородной задачи (28), а также собственные функции сопряженной задачи (33) и решения неоднородных задач второго порядка (42), (45), (47), (50). Для этого использовался метод стрельбы с применением процедуры ортогонализации.

3.2. Описание численных методов решения и результаты вычислений

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

(квг) 1 и (квгРг) 1 при старших производных. В

ходе интегрирования среди решений появляются

а

-1

осциллирующие и быстро растущие. При этом теряется линейная независимость частных решений [7].

Обозначая амплитуды возмущений функции тока, проекции скорости и температуры через ам можно переписать систему (28) в виде

dai ft \ az

г = 1... 8,

(65)

В рассматриваемой задаче Ю1 = щ, Ю2 = иу\,

Юэ = ^1, Ю4 = #1, © = щ', © = V', ®7 = ^1, © = 3 . Граничные условия на левой границе слоя определяют начальные значения для четырех из восьми функций ю^ Остальные четыре функции могут задаваться произвольно, что позволяет построить четыре линейно независимых решений системы (65) ©к) (здесь к - номер частного решения). Начальные значения произвольно задаваемых функции

а

,(*)

(г = 5...8) в разных решениях выбираются

г = 5... 8,

(66)

(и)

(67)

г = 1...8, А: = 1...4.

Потребуем, чтобы преобразованные функции

(г = 5... 8) были такими же, как и в начале интегрирования, т.е. удовлетворяли уравнениям (66). Получаем систему из 16 линейных уравнений для величин аЛ . С помощью вычисленных коэффициентов по формулам (67) можно определить остальные компоненты частных решений. Полученные значения функций а>)(г = 1... 8) принимаются в качестве начальных условий для следующего шага интегрирования.

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

требует проведения восстановительных операций для построения собственных функций [8].

Собственные функции рассматриваемой краевой задачи находятся как линейная комбинация частных решений:

щ = с^ + с2а>\2) + с343) + с4(44) = сщ . (68)

Для построения собственных функций по данной формуле необходимо провести над ортогона-лизованными решениями обратные преобразования. На некотором шаге интегрирования т выражение (68) приобретает вид

(69)

как компоненты четырех взаимно ортогональных векторов:

где 8Ч - символ Кронекера. В процессе интегрирования начальные направления этих векторов восстанавливаются дополнительным линейным преобразованием решений [7].

Пусть после некоторого шага интегрирования получены частные решения

аа/)3\ со^ (г =1,...,8). Построим четыре линейных комбинации этих решений:

где Ат-...-А2А1 - матрицы, составленные из коэффициентов % .

Для вычисления собственных значений и собственных функций написана программа на языке Fortran 90. Алгоритм программы построен согласно рекомендациям, сделанным в [8]. Сначала определяется собственные значения краевой задачи. Затем система дифференциальных уравнений интегрируется повторно. После каждой процедуры ортогонализации вычисляется матрица

Bm = AmBm_l, (B0 )j=Sj ,

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

3i(zm) = Bm1cdi(zm).

(70)

Коэффициенты ^ в формуле (68) находились из условий на правой границе:

z = 1: J(k)=а%)=а%)=$) = 0.

(71)

Получаем систему однородных алгебраических уравнений для с4:

1 (1) , (2) , (3) , (3) А

z = 1: c1со)) + е2со) ) + е3 со) ) + с3со)) = 0, i = 1.. .4.

(72)

Определитель данной системы всегда равен нулю, следовательно, значение одного из коэффициентов произвольно. Этим объясняется невозможность определения амплитуды течения в первом порядке. После вычисления остальных коэффициентов в точках ортогонализации определяются собственные функции Ж.

Необходимое для сходимости метода количество точек ортогонализации определялось опытным путем: оказалось, что необходимо порядка 300 точек ортогонализации. Для повышения точности определения собственных функций расчеты проведены при количестве точек ортогонализации, равном 500.

n=1

(а) (б)

Рис. 4. Кривые устойчивости при: а - Pr = 1; б - Pr = 0.1

Для нахождения обратных матриц, умножения матриц, решения систем линейных уравнений и вычисления определителя матрицы использованы функции библиотеки IMSL Fortran 90 MP версии 4.0. Интегрирование осуществляется методом Рун-ге-Кутта 4-5 порядка точности. При нахождении собственных значений критические числа Грасго-фа определялись с помощью метода одномерных секущих, а соответствующие волновые числа наиболее опасных возмущений вычислялись поиском минимума функции Gr(k), используя метод парабол.

В результате расчетов определены собственные числа линейной задачи устойчивости и собственные функции линейной и сопряженной задач. Собственные числа с хорошей степенью точности совпали с критическими числами, определенными в работе [3]. Построены кривые устойчивости, отражающие зависимость критического числа Грасгофа от числа Рейнольдса (рис. 4). Из графиков видно, что акустическое воздействие значительно снижает границу устойчивости плоскопараллельного течения.

Собственные функции линейной задачи также хорошо согласуются с видом вторичного течения, полученным в [3]. Акустическое воздействие приводит к существенному изменению надкритических режимов. В частности, подавляется вторичное течение в нижней части слоя. Более подробное исследование влияния акустического воздействия на вторичное течение проведено в ходе прямого численного моделирования.

Для нахождения собственных функций неоднородных задач второго порядка (42), (45), (47) удобно представить ее решение в виде суперпозиции общего решения однородной задачи и частного решения неоднородной:

W = а>] + c ®,(1) + c ®,(2) + С ®,(3) + С ®,(4), (73)

(1) (2) (3) (4)

где Щ) , Щ ), Щ ), Щ ) - частные решения однородной задачи, Щ - частное решение неоднородной задачи. Метод нахождения общего решения однородной задачи рассмотрен ранее.

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

é4. = + + Ь2со\2) + й3®,(3) + Й4®,(4) . (74)

Коэффициенты b подбираются таким образом, чтобы в преобразованном частном решении функции Щ, Щ, Щ, Щ были такими же, как и в начале интегрирования, т.е. нулевыми.

Для восстановления частного решения проводится преобразование, обратное (74). При этом необходимо учитывать, что частные решения однородной задачи также ортогонализуются, следовательно, перед восстановлением частных решений неоднородной задачи необходимо восстановить частные решения однородной задачи. Коэффициенты ci находились из граничных условий на правой границе.

Таким образом, используя результаты вычислений собственных функций линейной задачи, можно получить решение данной неоднородной задачи. При решении линейной задачи собственные функции определялись только в некотором

наборе точек, поэтому для интегрирования частного решения неоднородной задачи необходимо интерполировать собственные функции линейной задачи. Интерполяция осуществлялась квадратичным сплайном с использованием процедуры, входящей в математическую библиотеку IMSL Fortran.

Неоднородная система (50) состоит из дифференциальных уравнений, которые могут быть проинтегрированы аналитически. Для собственной функции скорости получаем

и20 = 2 ^

J

Л

+ Cz + B,

(75)

где коэффициенты С и В определяются из условия отсутствия возмущений на границах:

B = C = -Я| Jwud^

\

(76)

Собственная функция температуры определяется по следующей формуле:

^20 = Pr

J J u20d^d£ _-i-i

+ 2«(j w^d^ где K, L определяются из граничных условий:

+ Kz + L

(77)

K = L =-

1 i c ( i — J Ju20d^-^l J w^d^ 2 1 1 I 1

(78)

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

4. Прямое численное моделирование

4.1. Описание метода

Численное моделирование влияния акустического воздействия на надкритические режимы адвективного течения в слое проводилось методом конечных разностей. Получим нелинейные уравнения, описывающие конвективные течения рассматриваемого типа. Применив к уравнениям (12) и (13) операцию rot и спроектировав результат на ось х, исключим добавку к распределению давления p{y,z):

8Q 8Q SQ ^ 83 — + v — + v — = АП + Gr—,

8t 8y 8z 8y

(79)

где £ 2 = г<Л, V = су_ /су - су , - вихрь скорости.

Уравнение непрерывности позволяет ввести функцию тока, связанную с компонентами скорости следующим образом:

Vy 8z

V =-

8у

Z'

8у

8y

(80)

В результате, трехмерное поле скорости определяется двумя скалярными функциями: проекцией скорости на ось х и функцией тока:

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

С помощью рассмотренных численных методов могут быть найдены все функции, необходимые для определения коэффициентов, входящих в дифференциальные уравнения (52). По формулам (54)-(61) проведено вычисление этих коэффициентов в широком диапазоне значений параметров задачи. Установлено, что константа Ландау К всегда положительна, следовательно, характер ветвления решений вблизи порога мягкий.

Коэффициенты системы (52) ранее были получены в [3]. С использованием этих коэффициентов в [4] проведено численное моделирование уравнений (52); в частности, сделан вывод о возможности жесткого возникновения неустойчивости.

Однако вычисления, проведенные в настоящей работе, показали, что в действительности оба условия (63), (64) справедливы во всём диапазоне параметров, где существует спиральная мода неустойчивости. Это означает, что в задаче о термоакустическом адвективном течении в слое про-

v = Vj ех л--еу--ег = v^ ех + rot

8z 8y

(we,)-

Система нелинейных дифференциальных уравнений, описывающих эволюцию возмущений типа продольных валов, принимает вид

дО дш дО дш дО ^ ^ дЗ

— + —-------= ДО + Gr—,

дt дz ду ду дz ду

дЗ дш дЗ дш дЗ 1 „

— + V + —-------= —ДЗ,

дt х дz ду ду дz Рг (81)

дvv дш дvv дш дvr , ^

—- + —--------- = Д^ - + с,

дt дz ду ду дz

д ш = -О.

Граничные условия для этой системы запишутся так:

z = ±1: v = 0, д-^ = 0, 8z

v = Re, 3 = 0.

(82)

Вследствие сделанных предположений о форме вторичного течения решение данной системы бу-

z

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

2 = Цу^ йуй2 = 0,

& 2

2й ■

М = -

к2

4 + А\ш \'

I т тах I

на сетках 30x30, 40x40, 50x50, 60x60, 70x70. Течения, рассчитанные с использованием сеток 60x60 и более, имеют практически одинаковую интенсивность (см. таблицу).

Зависимость результатов вычислений от размера сетки. Расчеты проведены при Кг = 10, От = 700, Рт = 1

где - Ь горизонтальный размер расчетной области. К граничным условиям (82) следует добавить периодические условия на боковых границах расчетной области:

уу(0,2) = уу(Ь, 2), 3(0, 2) = 3(Ь, 2), 0(0, 2) = П(Ь, 2), уД0, Ь) = ух(Ь, 2).

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

При составлении конечно-разностного аналога уравнений (81) пространственные производные аппроксимировались центральными разностями, а производные по времени - по явной схеме второго порядка. Значения вихря на твердых границах находились по формулам Пирсона, полученным путем разложения функции тока в ряд Тейлора в приграничной точке с точностью до кубических членов:

' = к = '

Размер сетки Максимальное значение функции тока

30x30 1.136

40x40 1.097

50x50 1.079

60x60 1.070

70x70 1.069

где к - шаг сетки по оси 2 . Уравнение Пуассона для функции тока решалось методом установления.

Шаг по времени в соответствии с требованиями устойчивости счета выбирался по формуле [9]:

где - максимальное значение функции тока,

А - эмпирический параметр (в расчетах выбирался равным 3/2). В начальный момент времени поля скорости и температуры задавались в соответствии с формулами (8), а для функции тока и вихря скорости вводились малые случайные возмущения. При нахождении зависимости характеристик течения от какого-либо параметра задачи вычисления проводились с помощью метода продолжения по параметру, что позволяет значительно сократить время вычислений.

Вычисления проводились на сетке размером 60x60 узлов. Для определения необходимого размера сетки сравнивались результаты, полученные

4.2.Полученные результаты

Задача содержит четыре безразмерных параметра: число Рейнольдса (Re), определяющее интенсивность акустического воздействия, число Грасгофа ^г), характеризующее величину конвекции, число Прандтля (Рг) и горизонтальный размер расчетной области, определяющий пространственный период конвективной структуры.

Для классического адвективного течения реализуется трехвихревой режим течения с интенсивными пристеночными вихрями [6]. При увеличении числа Рейнольдса нижний вихрь подавляется, а верхний усиливается (рис. 5). При небольших надкритичностях изолинии функции тока у схожи с собственными функциями линейной задачи, обсуждавшимися в [3].

Построенные по результатам расчетов кривые подтверждают корневой закон зависимости интенсивности вторичного течения от числа Грасгофа и хорошо согласуются с результатами слабонелинейного анализа. На рис. 6 представлена амплитудная кривая вторичного течения - зависимость максимального значения функции тока ц/т от числа Грасгофа. Точками показаны результаты нелинейных расчетов, сплошная кривая получена в слабонелинейном приближении. Как видно, результаты, полученные двумя методами, неплохо согласуются друг с другом даже достаточно далеко от порога: при относительной надкритичности 1.5 отличие результатов чуть более 3%. Подобные амплитудные кривые позволяют определить критическое число Грасгофа с помощью экстраполяции по корневому закону.

Определенные по амплитудным кривым критические числа с хорошей точностью совпали со значениями, полученными в рамках линейной теории устойчивости (рис. 7). Расчеты проведены при Pr = 1, 0.5, 0.1. Расхождение значений критических параметров, определенных в ходе линейного исследования и прямого численного моделирования, не превышало 2%.

-1.00- —1.00— 1.00 1 .оо

3.50-3.5С

'1.00-1.00-1.00-1.00-

— -1.50--1.5С

0.00 0.20 0.40 0.60 0.80 1.00 1.20 1.4

1.60 0.00 0.20 0.40 0.60 0.80 1.00 1.20 1.40 1.60

Рис. 5. Линии тока и изотермы при Gr = 640, Re = 10, Pr = 1, к = 3.78

шп

0.6 1 1.4 10'3-°г

Рис. 6. Максимальное значение функции тока в зависимости от числа Грасгофа, Рг = 1, Re = 10, k = 3.78

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

Проследим за эволюцией вторичного течения при увеличении числа Грасгофа. При числах Грасгофа, близких к критическому, значение завихренности для верхних конвективных валов значительно больше, чем для нижних. Однако при возрастании Gr разница в интенсивности верхних и нижних валов уменьшается, и при достаточно больших значениях числа Грасгофа интенсивность верхних и нижних валов практически одинакова. Это объясняется тем, что, увеличивая число Грасгофа, мы усиливаем интенсивность конвективного течения, оставляя при этом неизменным

Рис. 7. Кривая устойчивости. Точками показаны критические значения числа Gr, полученные в ходе прямого численного моделирования. Сплошная линия - линейная теория устойчивости

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

При увеличении надкритичности стационарное решение становится неустойчивом, в слое возникают колебания сложной структуры. На рис. 10 приводится эволюция максимального значения функции тока за период колебаний. С ростом числа вг частота колебаний возрастает (рис. 11).

Колебания развиваются следующим образом. Сначала появляются два дополнительных конвективных вала в нижней части области (рис. 12, а), при этом максимальное значение функции тока снижается. Затем дополнительные конвективные валы формируются в верхней части области (рис. 12, б). Далее интенсивность дополнительных конвективных валов, находящихся вблизи верхней границы, уменьшается (рис. 12, в) и они исчезают.

0.80

0.80

0.60

0.60

0.40

0.40

0.20

0.20

0.20

-0.20

0.40

-0.40

0.60

-0.60

-0.80

4

2

(а) (б)

Рис. 8. Изолинии функции тока при Gr = 1000, Рг = 1: а - Кг = 10, к = 3.78; б - Кг = 20, к = 3.16

(а)

(б)

Рис. 9. Изолинии функции тока при Рг = 1: а - Кг = 30, От = 150, к = 2.16; б - Кг = 40, От = 80, к = 1.89

1. 2

/

0 . 8

0 . 4

0 1 2 3 4 t 5

°Н—1—I—1—I—1—I—1—I—1—I 1500 1600 1700 1800 1900 2000

Gr

Рис. 10. Изменения максимального значения функции тока за период колебаний при Gr = 1580, Re = 10, Pr = 1, к = 3.78

Рис. 11. Зависимость частоты колебаний от числа Грасгофа при Pr = 1, Re = 10, к = 3.78

0 20 0 40 0 60 0 80 1 00 1 20 1 40 1 6

0 20 0 40 0 60 0 80 1 00 1 20 1 40 1

(а)

(б)

(в)

00 0 20 0 40 0 60 0 80 1 00 1 20 1 40 1 6

(г)

(д)

(е)

7

6

5

4

3

2

0 40

0 20

-0 20

-0 40

-и

-и

-0

-0

Рис. 12. Эволюция изолиний функции тока за период колебаний при Gr = 1580, Re = 10, Рг = 1, к = 3.78 .• a - t = 0.13 ; б - t = 2.26; в - t = 3.80; г - t=3.89; д - t = 4.09; е - t = 4.72

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

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

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

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

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

1. Ниборг В. Акустические течения // Физическая акустика. Т. 2. М.: Мир. 1969. С. 302-377.

2. Гершуни Г. З., Жуховицкий Е. М., Непомнящий А. А. Устойчивость конвективных течений. М.: Наука, 1989. 320 с.

3. Любимов Д. В., Шкляев С. В. Об устойчивости адвективного термоакустического течения // Известия РАН. Механика жидкости и газа. 2000. № 3. С.10-21.

4. Седельников Г. А. Об устойчивости спирального вторичного течения в задаче термоакустической конвекции // Мат. X Всерос. студ. конф. «Математическое моделирование в естественных науках». Пермь, 2001. С. 38-39.

5. Любимов Д. В. О тепловой конвекции в акустическом поле // Известия РАН. Механика жидкости и газа. 2000. № 2. С. 28-36.

6. Мызников В. М. Конечно-амплитудные конвективные движения жидкости в горизонтальном слое с продольным градиентом температуры //. Мат. модели течений жидкости: труды VI Всес. семинара по численным методам механики вязкой жидкости. Новосибирск: Изд-во ин-та теор. и прикл. мех. СО АН СССР. 1978. С. 176-186.

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

8. Бирих Р. В., Рудаков Р. Н. Применение метода ортогонализации в пошаговом интегрировании при исследовании устойчивости конвективных течений. Часть I // Гидродинамика. 1974. Вып. V. C.43-46.

9. Гершуни Г. З., Жуховицкий Е. М., Юрков Ю. С. Конечно-амплитудные конвективные движения в прямоугольных полостях с внутренними источниками тепла // Гидродинамика. 1974. Вып. V. С. 3-23.

References

1. Nyborg W. Acoustic streaming. In: W. Mason (Ed.) Physical acoustics, vol. 2B. New-York: Academic Press. 1965, pp. 290-295.

2. Gershuni G. Z., Zhukhovitskij E. M., Nepom-njaschij A. A. Ustojchivost' konvektivnyh techenij (Stability of convective flows). Moscow: Nauka, 1989, 320 pp. (In Russian).

3. Lyubimov D. V., Shklyaev S. V. Thermal convection in an acoustic field. Fluid Dynamics, 2000, vol. 35 (3), pp. 321-330.

4. Sedelnikov G. A. Ob ustoichivosti spiral'nogo vtorichnogo techeniia v zadache termoakustich-eskoi konvektsii (On stability of spiral secondary flow in thermoacoustic convection problem). Proc. of 10th All-Russian student conf. "Mathematical modeling in science", Perm, 2001, pp. 38-39. (In Russian).

5. Lyubimov D. V. Thermal convection in an acoustic field. Fluid Dynamics, 2000, vol. 35 (2), pp. 177184.

6. Myznikov V. M. Konechno-amplitudnye konvek-tivnye dvizheniia zhidkosti v gorizontal'nom sloe s prodol'nym gradientom temperatury (Finite-amplitude convective motions of liquid in horizontal layer with longitudinal thermal gradient). Mathematical model of fluid flows: Proc. of 6th AllUnion Seminar on numerical methods in viscous fluid mechanics. Novosibirsk: Institute of theoretical and applied mechanics SB USSR, 1978, pp. 176-186. (In Russian).

7. Birikh R. V., Rudakov R. N., Semakin I. G. Primenenie metoda ortogonalizatsii v poshagovom

integrirovanii pri issledovanii ustoichivosti konvektivnykh techenii. Chast' I (Application of orthogonalization method in step-by-step integration for convective flows stability problems. Part I) Convective flows. Perm: Perm State Pedagogical Institute, 1979, pp. 58-60. (In Russian).

8. Birikh R. V., Rudakov R. N. Primenenie metoda ortogonalizatsii v poshagovom integrirovanii pri issledovanii ustoichivosti konvektivnykh techenii. Chast' I (Application of orthogonalization method

in step-by-step integration for convective flows stability problems. Part I) Hydrodynamics, 1974, vol. 5, pp. 43-46 (In Russian).

9. Gershuni G. Z., Zhukhovitskij E. M., Yur-kov Yu. S. Konechno-amplitudnye konvektivnye dvizheniia v priamougol'nykh polostiakh s vnu-trennimi istochnikami tepla (Finite-amplitude con-vective motions in rectangular cavities with volume heat sources). Hydrodynamics, 1974, vol. 5, pp. 3-23 (In Russian).

Просьба ссылаться на эту статью в русскоязычных источниках следующим образом:

Иванцов А. О. Слабонелинейный анализ устойчивости термоакустического адвективного течения // Вестник Пермского университета. Физика. 2019. № 3. С. 28-44. doi: 10.17072/1994-3598-2019-3-28-44

Please cite this article in English as:

Ivantsov A. O. Weakly non-linear analysis of thermoacoustic advective flow. Bulletin of Perm University. Physics, 2019, no. 3, pp. 28-44. doi: 10.17072/1994-3598-2019-3-28-44