Устойчивость сверхзвукового течения Куэтта колебательно-возбужденного двухатомного газа Текст научной статьи по специальности «Физика»

2015 Математика и механика № 1(33)

УДК 532.5:532.517.4 Б01 10.17223/19988621/33/5

И.В. Ершов

УСТОЙЧИВОСТЬ СВЕРХЗВУКОВОГО ТЕЧЕНИЯ КУЭТТА КОЛЕБАТЕЛЬНО-ВОЗБУЖДЕННОГО ДВУХАТОМНОГО ГАЗА1

В рамках линейной теории исследована устойчивость течения Куэтта колебательно-возбужденного двухатомного газа с параболическим профилем статической температуры. Исходной математической моделью течения газа служила система уравнений двухтемпературной аэродинамики. В результате было показано, что при определенном сочетании значений параметров исследуемого течения (чисел Рейнольдса Яе, Маха М, объемной вязкости аь степени колебательной неравновесности уу1Ь и времени колебательной релаксации т) оно может быть как устойчиво, так и неустойчиво по отношению к малым возмущениям. Для вязких возмущений рассчитаны спектры собственных значений, инкременты нарастания и кривые нейтральной устойчивости в плоскости (Яе, а) для первой и второй растущих мод в диапазоне чисел М = 2-6 и Яе = 104-107. Найден диапазон изменения критических чисел Рейнольдса Яесг и (2-5)-104. Показано, что при всех уровнях возбуждения наиболее неустойчивой является вторая мода. Возбуждение практически не меняет форму области неустойчивости, но ее границы с ростом возбуждения смещаются в сторону больших волновых чисел. Можно констатировать, что, в общем, возбуждение внутренних степеней свободы молекул газа снижает инкременты нарастания возмущений и оказывает стабилизирующее воздействие на течение.

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

В работах [1, 2] устойчивость плоского дозвукового течения Куэтта термически неравновесного молекулярного газа рассматривалась на основе нелинейной энергетической теории. Проведенное в этих работах обобщение теории на случай сжимаемых течений позволило получить значения критических чисел Рейнольдса Яесг, в том числе и для слабо неравновесного газа. Найденные значения Яесг по порядку величины совпадают с критическими числами Рейнольдса, полученными в аналогичной постановке для несжимаемого течения [3]. Этот результат подтверждает известное представление о том, что дозвуковое течение Куэтта можно считать практически несжимаемым. Вместе с тем экспериментальные данные по Яесг в обоих случаях превосходят расчетные значения на несколько порядков. При этом для несжимаемой жидкости в настоящее время отсутствуют подходы, которые позволили бы сблизить данные теории и эксперимента. До последнего времени единственной альтернативой энергетической теории была классическая линейная теория устойчивости. В ее рамках плоское течение Куэтта несжимаемой жидкости изучалось многими авторами. Краткий обзор работ в этом направлении приведен в [4]. Центральным здесь является строгий математический результат

1 Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (проект № 14-01-00274).

[5] о его абсолютной устойчивости при всех числах Рейнольдса и произвольных длинах волн возмущений.

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

Основные численные результаты во всех трех работах получены методом кол-локаций с использованием QZ-алгоритма для нахождения спектра фазовых скоростей возмущений. Тем не менее результаты работ [6-8] противоречат более ранним результатам [4]. Авторы [4] констатировали сильное стабилизирующее влияние вязкости и отсутствие растущих вязких мод вплоть до чисел Яе = 5-106 при числах Маха М < 5. Отсутствие растущих вязких мод было также зафиксировано на основе асимптотических поправок к результатам в невязком пределе. В то же время в численных расчетах [6-8] были найдены неустойчивые вязкие моды при близких числах Рейнольдса. Более того, авторами [6, 7] было обнаружено, что в некотором диапазоне длин волн, чисел Рейнольдса и Маха вязкость оказывает дестабилизирующее воздействие. В частности, возникает неустойчивость выделенной моды, устойчивой в невязком пределе.

Возможной причиной такого расхождения является несовершенство реализации численного метода в [4], где использовалась авторская разработка, в отличие от работ [6-8], применявших профессиональное математическое обеспечение, которое в [6, 7] дополнительно тестировалось на основе альтернативного конечно-разностного метода.

Общие характеристики линейной устойчивости плоскопараллельных течений колебательно-возбужденного газа рассматривалось в [10, 11], где было показано значительное стабилизирующее воздействие релаксационного процесса. Линейная устойчивость течения Куэтта в условиях сильного отклонения от термодинамического равновесия до последнего времени не исследовалась. Следует отметить, что в цитированных работах влияние объемной вязкости, отражающей слабую неравновесность внутренних степеней свободы молекул газа, исключалось с помощью соотношения Стокса. Поэтому обращение к линейной теории с целью исследования влияния термической неравновесности на характеристики устойчивости классического течения представляет самостоятельный интерес. Результаты, полученные на основе линейной теории для невозбужденного газа, в принципе позволяют надеяться на сближение по порядку величины расчетных и экспериментальных значений критических чисел Рейнольдса по сравнению с результатами [1, 2], по крайней мере, для сверхзвуковых чисел Маха.

Постановка задачи и основные уравнения

Рассматривается линейная устойчивость плоского течения Куэтта термически неравновесного двухатомного газа. В координатной плоскости (х, у) поток ограничен двумя бесконечными параллельными плоскостями, находящимися на расстоянии И друг от друга. Считается, что плоскость у = 0 покоится, а граница у = И движется равномерно в собственной плоскости со скоростью и0. Исходной математической моделью течения газа служит система уравнений двухтемпературной аэродинамики. В соответствии с физическими представлениями [12, 13] модель двухтемпературной аэродинамики является общепринятой физико-математической моделью течений колебательно-возбужденного молекулярного газа, когда диссоциацией молекул, возбуждением верхних колебательных уровней и поправками на ангармонизм колебаний можно пренебречь.

В качестве характерных величин для обезразмеривания были выбраны ширина канала И, скорость границы и0, плотность р0 и температура Т0 основного течения на движущейся границе канала и образованные из них время т0 = Ь/П0 и давление р0 = р0и02. В безразмерных переменных система уравнений двухтемпературной аэродинамики имеет вид

др + д р щ

д / д X:

(

= 0, р

Л

д щ д и —- + и :—:

, д г 1 д х,

V 1

д р 1

д 2щ

1

д х Яе д х

2 + Яе

1

а, +-

д2и,

3 ) дх,дх

\д Т + 3 Т ^ 5 щ

р1 ^ + и> ^1 +(у - 1)р Т =

д г д х: ) д х:

у д2Т + ууР (Ту - Т) + у(у -1)

ЯеРг д х2

2Яе

^ +ди± ^ д хд х,

V 1 1 )

2

+* * - Ж

5 Ту 5 Ту

20 ууу д2Ту УуР (Ту - Т)

ЗЗЯеРг дх*

у М р = р Т, уу =

У уь

1 - У УЬ

1,1 = 1, 2,

(1)

где х! = х, х2 = у, а по повторяющимся индексам подразумевается суммирование. Параметры, входящие в уравнения системы (1), определяются следующим образом. Коэффициент а! = пь/п есть отношение объемной и сдвиговой вязкостей. Коэффициент у = Ср/Су - показатель адиабаты, су = су су г, ср = су+Я - соответственно удельные теплоемкости при постоянных объеме и давлении, где выделены составляющие, обусловленные поступательным су г и вращательным су г движением молекул газа, Я - газовая постоянная. Коэффициент уУ1Ь = су у/(су 4+ су г+ су у) характеризует степень неравновесности колебательной моды, су у - удельная теплоемкость при постоянном объеме, связанная с колебательным движением молекул газа, т - характерное время колебательной релаксации. Параметры Яе = р0Ии0/п и М = и0/(уЯТ0)12 есть соответственно числа Рейнольдса и Маха несущего потока. Рг = псу/Х - число Прандтля, где коэффициент теплопроводности X = X (+Х г определяется поступательными и вращательными степенями свободы молекул газа.

Нижний предел уУ1Ь = 0 соответствует случаю невозбуждения колебательной моды молекул. С другой стороны, равнораспределение энергии по степеням свободы молекул не является здесь верхним пределом для параметра уУ1Ь. Поскольку закон равнораспределения энергии неприменим в неравновесной ситуации, описываемой системой уравнений (1), когда разрыв между статической температурой потока Т и колебательной температурой ТУ1Ь может быть достаточно велик. В [13] показано, что при Т = 300 К неравновесная теплоемкость сУ1Ь и 1,8^. Используя равнораспределение энергии в состоянии термодинамического квазиравновесия по поступательным и вращательным модам молекул, получаем, что параметр уУ1Ь и 0,42. С ростом разрыва между температурами ТУ1Ь и Т значение уУ1ь увеличивается, приближаясь в пределе к единице, когда энергия колебательной моды молекул существенно превышает температуру квазиравновесного термостата, определяемого поступательными и вращательными степенями свободы молекул. В расчетах максимальное значение параметра уУ1Ь было выбрано равным 0,4 для того, чтобы остаться в рамках используемой модели, избежав возбуждения высоких колебательных уровней энергии.

Предполагается, что в невозмущенном стационарном потоке все параметры зависят только от поперечной координаты у, статическая Т5 и колебательная ТУ1Ь,5 температуры равны: Т5 (у) = ТУ1Ь 5 (у). Для исходного течения ставятся следующие

граничные условия:

и5(0) = 0, и5(1) = 1,

йТ?

йу

= 0, Т5 (1) = 1.

у=0

В результате получаем, что точное решение системы (1) имеет вид

(у - 1)РгМ2 2

и5 (у) = у, Т8 (у) = Туь, ? (у) = 1 + ■-(1 - у2),

Р?(у) =

Т? (у)

Р5 ( у) =

у М2

(2)

Подстановка мгновенных значений гидродинамических переменных возмущенного потока в виде

й1 = и5 + йх , и2 =й у , Р = Р5 + Р, Т = Т5 + T, ТУ1Ь = ТУ

у1ь, 5 + ту1ь ,

Р = Р5 + Р

(где йх (х, у, /), й у (х, у, /), Р( х, у, ^ , р( х, у, /), Т( х, у, /), Т^ х, у, /) - возмущения

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

др д /

+и

др ( дй

д й у

5 д х Р? х д у

+ й

= 0,

Р? I

дйх д /

-+и5

дйх

д х

-+и,.

диА

д у

=-дР+—

д х Яе

^2 й х

д 2й ^

д х д у

дР 5 у 5 у

+—I а, +-Яе1 1 3

2 -д и х

д 2й у

дх2 дхду

р5&+^^1=-дР+±

^ дг д х ) д у Яе

( д2 й у д 2й

2- Л

дх2 ду2

+—I а, +-Яе1 1 3

д2д и

2

дудх ду2

(3)

дТ_ д г

д_

д

р5 I-+ и8— + й у |ч| , ^

д х д у ) ( д х д у

+ у(у - 1)М2 р8\ +

д йх д г<г

Т

2

ЯеРг

д2 _ д2 _

2 ^Л

д х2 д у д_

Ту Р5

д _

у1Ъ

Ту Р5 (_*Ъ - Ъ , 2 у (у - 1)М2 (д йх , д 1 йП,

Яе

ду дх ) йу

д г

+ие

у1Ъ

д_е1

д х

+ и у

д у

20 уу у

ЗЗЯеРг

д2 _уъ , д 2 1 ту р5 (_уъ - г)

уЛ

д х2

д у 2

у М2 р = р5_ + р _5.

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

Представляем периодические по х возмущения в виде бегущих плоских волн:

Ч (х, у, г) = qo( у) • ег

а(х-сг)

Ч(хy,г) = (йх, йу, р _, _у1ъ,р) =

Чо(у) = (й, аV, р, 0, еу, р),

(4)

где а - волновое число вдоль периодической переменной х, с = Оу+к^ - комплексная фазовая скорость, / - мнимая единица. Подставляя (4) в уравнения системы (3), получим, что амплитуды возмущений будут удовлетворять следующим уравнениям:

1 а

Бр + ар^у + р5с = 0, —Д й -р5Бй - ар5уи'5 -/ае = 0, —Ду - ар5Бу - е'= 0, Яе Яе

7 Д0 - р5Б0 - ар5у_' - (у -1)с + 2у(Т 1)М

ЯеРг

Яе

( + /а2V) + Ь^(0у -0) = 0, (5)

зЗЯерГ Д 0у - YvPsD0v - атурзУ_' -^(0у - 0) = 0, ТМ2р = р5 0 + р_5

й\ п = й , = V п = V , 1у=0 1у=1 1у=0 1у=1

где

0 = 0| 1 = 0 I 0 = 0у| 1 = р| 0 = р| 1 = 0,(6)

1у=0 1у=1 у 1у=0 у 1у=1 г 1у=0 г 1у=1 > 4 '

Б = /а(из - с), с = а(у' + ш),

е = р-

Яе

1

а1 +-

Ту =

Т уъ

1 - Т у1Ъ

Д = ■

й2

йу 2

--а

а штрих у функций здесь и далее обозначает дифференцирование по переменной у. Система (5) вместе с однородными граничными условиями (6) определяет спектральную задачу, в которой собственными значениями являются комплексные фазовые скорости возмущений с = сг+1с,, а числа Маха М, Рейнольдса Яе и волновое число а служат параметрами.

Методы решения спектральной задачи

Для расчета собственных значений c = cr+ici неустойчивых мод система уравнений (5) сводилась к матричному виду и далее решалась численно в среде пакета Matlab. Использовался метод коллокаций, основанный на полиномиальной интерполяции собственных функций полиномами Чебышева [14, 15]. В качестве узлов коллокации (интерполяции) выбирались точки Гаусса - Лобатто

1 L ( п п 1 + cos I

Уп =-

N

п = 0,1,

, N,

в которых полином Чебышева Ж-й степени имеет экстремумы на отрезке у = [0, 1]. Дифференциальные операторы первого порядка, входящие в спектральную задачу, аппроксимируются на данном шаблоне матрицей коллокационных производных размером (Ж+1)х(Ж+1), элементы которой определяются по формулам [14, 15]

DN,ij =

(-1)

/+1

si( у- - У]) _ У]

2(1-У 2 )

Si =

2 N2 +1 6

2 N2 +1

{2

j = 0, N,

j = 1,2, ..., N -1.

* * ], 1 < I = ] < N -1,

I = ] = 0, I = ] = N,

При этом элементы 1-й строки матрицы являются коэффициентами разностной аппроксимации первой производной в 1-м узле коллокации на шаблоне {ук}, к = 0, 1, ..., N. Дифференциальные операторы второго порядка аппроксимируются суперпозицией ^ = [14, 15].

В терминах введенных аппроксимаций задача (3) сводится к обобщенной задаче на собственные значения (линейному матричному пучку) относительно спектрального параметра с = ег+1е{.

X (Gkj -cFkj)Г] = 0, k = 0,1,2, ...,5N + 4.

(7)

j=0

В (7) вектор неизвестных г размером 5(Ж+1) состоит из значений собственных функций в узлах коллокации:

r (x2) = ^ Pl,., Р N , u1

, uN, v0, vb

vv,0'vv,1'

),

а матрицы О, Е размером 5(Ж+1)х5(Ж+1) вычисляются с использованием специальной процедуры МаИаЪ по формулам

О = А1 ® ОЖ + А2 ® ОЖ + А3 ® , Е = А4 ® ^ .

Здесь знак «®» обозначает прямое (тензорное) произведение матриц [16]; -единичная матрица размером (Ж+1)х(Ж+1); А4 (] = 1, 2, 3, 4) - матрицы размером 5х5:

(0 0 1

А =

0

Яе

а

о о —| а +-

Яе I 1 3

0 0 о

0 0

0

ЯеРг 0

0 0

0

0

20 УУ. 33ЯеРг

А2 =

0 0

Т*

у М2

0 0

/а ( 1

—I а! +-

Яе I 1 3

2у(у - 1)М2^;

Яе 0

-Р5 /а2 ( 1

—I а, +-

Яе I 1 3

0

-а (у -1) 0

0 0^ 0 0 Р5

0

у М2

00 0 0

( -iUs ' а Т5

А3 =

у М2 Т'

у М2 0

0

-'Р5 -а а1

0

-/а (у -1) 0

-Р5 -аР8и'8

-а2 а2 а а5

-аУ. Р5Т5

0

_'аР 5 у М2

у М2

-а3

У.Р5

0 А 0

0

У .Р,5 т

- Ту а4

А4 =

(-' 0 0 0 0

0 -'Р5 0 0 0

0 0 -'а2Р 5 0 0

0 0 0 -'аР5 0

1 0 0 0 0 -'аУ.

Л

а ( 4А . тт а . тт уа у.Р5

а1 = 1 а1 + т| + фЛ . а2 — + 'РЛ . а3 = +-+ /аР5и5 =

Яе V 3) Яе Яе Рг т

а4 = + £«. + /а^ , а5 = 2""('- 1)М2^ -рд. Т. =■ Тл

33ЯеРг т

Яе

1 - Уу1Ъ

Однородные граничные условия (6) для уравнения (7) учитываются неявно через оператор Б1Ы и на дискретном уровне реализуются заменой матриц БN (к = 1, 2) на окаймленные матрицы размером (Ж—1)х(Ж-1) [14, 15], которые получаются при выполнении условий

Б - Б - 0 Б - Б - 0

^0, ] ~ иы, ] - и> \,0 ~ ,,N - и>

I - 0,1,...,N, у - 0,1,...,N, £ -1,2.

Для нахождения всех собственных значений и функций обобщенной спектральной задачи (7) использовалась процедура МаНаЪ, реализующая р2-алгоритм, который позволяет одновременным ортогональным преобразованием привести пару матриц О, Е к обобщенной верхней треугольной форме [17]. В результате применения данной процедуры для фиксированных значений чисел Рейнольдса Яе и Маха М, объемной вязкости аь степени неравновесности колебательной энергии уу1Ь, времени колебательной релаксации т и волнового числа а получается набор (N+1)^ собственных значений с = ег+1е

Для проверки точности вычислений параллельно были проведены расчеты собственных значений с = сг+1с, с помощью метода «стрельбы». Для этого уравнения (5) заменялись фундаментальной системой уравнений и граничными условиями для вещественных и мнимых частей функций р, и, V, 6 и 6у1Ь. Полученная система при фиксированных наборах параметров Яе, М, уу1Ь, т и а интегрировалась численно с помощью процедуры Рунге-Кутты четвертого порядка на интервалах >>6 [0; 0,5] и >>6 [0,5; 1] с шагом Ду = 10—3. Шаг по волновому числу Да = 10—3. Точкой «прицеливания» служила середина канала у = 0,5. Значения сг и с, подбирались таким образом, чтобы вычисленные «слева» и «справа» в точке у = 0,5 значения функций р,, иг, V,., 6,, 6у1Ь, г и р,, и, V, 6,, 6у1Ь, , совпадали с точностью до 10—8. Соответствующее такому совпадению значение с принималось в качестве собственного значения при заданном наборе параметров Яе, М, а!, уу1Ь, т, а. Сравнение результатов, полученных с помощью методов коллокаций и «стрельбы», показало, что различия в значениях с = с,+,с,, наблюдаются лишь в шестом-седьмом десятичных знаках после запятой. Таким образом, была обеспечена необходимая точность вычисления инкрементов (декрементов) возмущений.

Расчеты велись при следующих значениях параметров: уУ1Ь = 0-0,4; т = 10—2—10; а! = 0—2; М = 0,5—15; Рг = 3/4; у = 7/5. Значение волнового числа менялось в диапазоне а = 0—10 с шагом Да = 10—3. Число узлов коллокаций в интервале [0, 1] варьировалось в диапазоне от N+1 = 100 до N+1 = 500 и в большинстве расчетов принималось равным N+1 = 300.

Результаты расчетов

Классификация вязких мод возмущений на четные и нечетные осуществляется аналогично классификации этих мод в невязком приближении [4] и сохраняется для случая колебательно-возбужденного газа [11]. Параметрические расчеты спектральной задачи показали, что изменение значений времени колебательной релаксации в диапазоне 10—2 < т < 10 слабо влияет на поведение кривых сг(а,аьуУ1ь, М, Яе) и с,(а,аьуУ1ь, М, Яе). Поэтому ниже расчетные данные приведены для одного значения характерного времени т = 1.

Расчет нейтральных кривых для п-й вязкой моды возмущений проводится следующим образом. Для фиксированных значений параметров аь уУ1Ь и числа Маха М вычисляются двухмерные массивы инкрементов (декрементов) п-й вязкой моды юПд (а-, Яек) = а- сП (а-, Яек), где одномерные массивы а,-, Яек рассчитываются по формулам а- = а0+/ Да (- = 0, 1, ..., I) и Яек = Яе0+к ДЯе (к = 0, 1, ..., К), Да и ДЯе - шаги по волновому числу и числу Рейнольдса соответственно. Массив Ю— (а-, Яек) определяет поверхность ю,(а, Яе) для п-й вязкой моды возмущений. Координаты точек, определяющих геометрическое место данной изолинии поверхности ю,(а, Яе) на плоскости (а, Яе), вычисляется по формуле | юП -к (а-, Яек) - С | < 10-8, где С представляет собой некоторое заданное числовое

значение. При С = 0 получаем нейтральную кривую ю,(а, Яе) = 0 для п-й вязкой моды возмущений. Используя подход, описанный выше для случая, когда фиксируются значения параметров а1, уУ1Ь и Яе, получим поверхность ю,(М, а) и нейтральную кривую ю,(а, М) = 0 на плоскости (а, М) для п-й вязкой моды возмущений.

Графики зависимостей фазовых скоростей сг(а) для семейств четных и нечетных мод возмущений приведены на рис. 1, где сплошной и штриховой линиями показаны зависимости сг(а) для четных и нечетных невязких мод возмущений соответственно при уУ1Ь = 0 и уУ1Ь = 0,4. Из данного рисунка, следует, что варьирование числа Рейнольдса Яе при фиксированных значениях параметров М, а1, уУ1Ь не оказывает влияния на сг(а). Из поведения кривых на рис. 1 следует, что рост значений параметра уУ1ь приводит к уменьшению для сг < 0 и сг > 1 и к возрастанию для сге(0; 1) по модулю значений фазовых скоростей сг, но при этом их значения не выходят за границы сг = 0 и сг = 1. Для нечетных мод при а^-0 сг > 1 и для всех мод, кроме моды I, сг ^ +да. В то же время для четных мод при а^-0 сг < 0 и для всех мод, за исключением моды II, сг ^ -да. Выделенные моды I и II при а = 0 имеют конечные пределы.

2

-2

4X1; I 7'

Ы\г а

2

4 а

Рис. 1. Фазовые скорости с,(а) при а! = 0 (а - М = 2, б - М = 5; 1,1' - мода I, 2, 2' - мода II, 3, 3' - мода III, 4, 4' - мода IV, 5, 5' - мода V, 6, 6' - мода VI, 7, 7' - мода VII, 8, 8' -мода VIII; о - Яе = 5-105, х - Яе = 5-106; сплошная линия - невязкие моды при уу1ь = 0, пунктирная - невязкие моды при уу1ь = 0,4)

1

0

Расчеты зависимостей сг(а) для четных и нечетных вязких мод возмущений показали, что первыми пересекают границы интервала сге [0; 1], в котором в соответствии с первым условием Рэлея [9-11] возможно развитие неустойчивости, вязкие моды I и II (см. рис. 1). При этом, как мы видим, моды I и II пересекают линию сг = 1 для М = 2 при а и 3,4, а для М = 5 при а и 1,9. Значения сг для старших четных (IV, VI и т.д.) и нечетных (III, V и т.д.) мод возмущений в диапазоне волновых чисел 0 < а < 3,4 для М = 2 и при 0 < а < 1,9 для М = 5 лежат вне интервала неустойчивости сге [0; 1]. Из рис. 1 следует, что рост числа Рейнольдса Яе не оказывает влияния на изменение фазовых скоростей сг как четных, так и нечетных мод возмущений.

а) Влияние вязкости на инкременты наиболее неустойчивых вязких мод возмущения в совершенном газе. На рис. 2 и 3 демонстрируется влияние варьирования числа Рейнольдса Яе на поведение зависимостей инкрементов ю,- = /ас,(а) для наиболее неустойчивых вязких мод I и II в совершенном газе (а! = уУ1Ь = 0) при фиксированных числах Маха М = 2 и М = 5 соответственно. Из графиков рис. 2 следует, что при М = 2 для широкого интервала изменения числа Рейнольд-са Яе значения ю,- возрастают с ростом Яе, но остаются в отрицательной полуплоскости ю,- во все диапазоне изменения волнового числа а.

0

-4

-12

1

" -____ -----

——--- -Л \ \

а

2

з -2["

Л __У \ ' - - -

б

3

4

5

6

3

4

5

6

Рис. 2. Зависимости <в,(а) для моды I (а) и моды II (б) в совершенном газе при М = 2 (линия в виде точек - Яе = 107, штрихпунктирная - Яе = 106, штриховая - Яе = 105, сплошная -Яе = да)

0

4

а

а

Из рис. 3 следует, что графики зависимостей ю (а) для наиболее неустойчивой моды II при М = 5 и различных значениях Яе имеют два четко выраженных «пика», ширина первого из которых существенно меньше ширины второго. С ростом значений Яе, как мы видим, ширина второго «пика» увеличивается, а ширина первого «пика» уменьшается. При этом высота первого «пика» падает, приближаясь в пределе Яе^-да к нулю, а высота второго «пика» наоборот растет. Причем рост значений ю в диапазоне волновых чисел а, принадлежащих области второго «пика», ограничен сверху кривой, полученной в невязком приближении Яе^-да (сплошная линия на рис. 3). В полуплоскости ю,- < 0 рост значений Яе сопровождается ростом значений ю .

а

Рис. 3. Зависимости <в,(а) для моды II в совершенном газе при М = 5 (линия в виде точек - Яе = 5-105, штрихпунктирная - Яе = 5-106, штриховая - Яе = 5-107, сплошная - Яе = да)

На основании сказанного выше можно констатировать, что, в общем, рост значений числа Рейнольдса Яе приводит к росту значений инкрементов ю,- наиболее неустойчивых вязких мод возмущений и тем самым оказывает дестабилизирующее воздействие на течение. Отметим, что ранее аналогичный результат был также получен в работах [6-8].

б) Влияние термической неравновесности газа на инкременты наиболее неустойчивых вязких мод возмущений. На рис. 4 представлены изолинии инкрементов роста ю,(М, а) вязких мод I и II для числа Рейнольдса Яе = 2-105 и различных значений объемной вязкости а! и степени колебательной неравновесности уУ1Ь. Сплошные линии соответствуют совершенному газу, а штриховые -а! = 2, уУ1Ь = 0,4. Замкнутые линии ю,(М, Яе, а, а!, уУ1Ь) = 0 представляют собой нейтральные кривые устойчивости первой и второй вязких мод возмущений на плоскости (М, а) для заданных значений параметров Яе, аь и уУ1Ь. Внутри областей, ограниченных линиями ю,(М, Яе, а, аь уУ1Ь) = 0, значения ю,- > 0 (течение неустойчиво), а вне - ю,- < 0 (течение устойчиво). Точки Л\, Л\, Л2, Л'2 и Б\, В\, В2, В'2 на нейтральных кривых ю,(М, Яе, а, аь уУ1Ь) = 0 представляют собой соответственно левые и правые критические точки на плоскости (М, а).

Из рис. 4 следует, что рост параметров а!, и уУ1Ь не меняет форм областей неустойчивости мод I и II, но их границы с ростом возбуждения внутренних степеней свободы молекул газа смещаются в сторону больших волновых чисел а. При этом абсциссы точек Ль Л2 на нейтральных кривых ю,(М, Яе, а, аь уУ1Ь) = 0 смещаются в сторону больших значений числа Маха М (абсциссы точек Л'ь Л'2), а абсциссы точек Вь В2 - в сторону меньших значений М (абсциссы точек В'ь В'2). Следовательно, с ростом значений параметров а!, уУ1Ь размеры областей неустойчивости вязких мод I и II уменьшаются.

3,5

а 3

2,5

Л % э а

XV чХ^З

4,8

а 3,1

1,4

2,4

3

М

3,6

Рис. 4. Изолинии инкрементов роста ю,(М, а) при Яе = 2-105 (а - мода I, б - мода II; сплошная линия - совершенный газ, штриховая - а1 = 2, у„ь = 0,4)

Влияние варьирования параметров аь уУ1Ь на поведение инкрементов роста ю,(а) представлено на рис. 5. Из графиков на рис. 5 следует, что при Яе = 5-105 и М = 3 кривые ю,(а) как для моды I, так и для моды II достигают своих максимальных значений, которые лежат в полуплоскости ю,- > 0. С ростом а1, уУ1Ь значения максимумов падают, а значения волновых чисел а, при которых они достигаются, смещаются в область больших а. При этом мы наблюдаем, что на достаточно узких интервалах изменения волнового числа 2,68 < а < 2,78 для моды I и 2,29 < а < 2,40 для моды II рост значений параметров а1, уУ1Ь приводит к дестабилизации первой и второй мод возмущений, а вне указанных интервалов к их стабилизации.

Рис. 5. Зависимости инкрементов роста ю,(а) при Яе = 5-105 и М = 3 (а - мода I, б - мода II; сплошная линия - совершенный газ, штриховая - а1 = 2, уУ1Ь = 0,4)

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

в) Критические параметры течения. Влияние сжимаемости и неравновесности внутренних мод молекул газа. На рис. 6 представлены изолинии инкрементов роста ro,(Re, а) вязких мод I и II для различных значений числа Маха M и параметров аь yvlb. Сплошные линии соответствуют случаю ai = yvlb = 0, а штриховые - а1 = 2, yvlb = 0,4. Точки K1, K\, K2, K'2 и K3, K'3 на нейтральных кривых устойчивости ю,(а, Re) = 0 представляют собой критические точки для вязких мод I и II, а их координаты определяют критические значения числа Рейнольдса Recr и соответствующие им значения волновых чисел acr.

Ö 2,5

к\ а I

6,3

Ö 3,8

104

5105 Re

1,3

106 104

5-105 Re

Рис. 6. Изолинии инкрементов роста ю,(Ке, а) при Яе = 5105, М = 3 (а) и М = 5 (б) (сплошная линия - совершенный газ, штриховая - а! = 2, уУ1Ь = 0,4)

3

2

Из рис. 6 следует, что с ростом возбуждения внутренних степеней свободы молекул газа (параметров аь уУ1Ь) нейтральные кривые = 0 не меняют своей формы, но смещаются в сторону больших значений чисел Рейнольдса Яе и волновых чисел а. В результате рост параметров а! и уУ1Ь приводит к росту значений Яесг и асг. Так, например, критические значения числа Рейнольдса и соответствующие им значения волновых чисел вязкой моды I при числе Маха М = 3 равны Яесг: = 123900, асг1 = 2,837 для а! = уУ1Ь = 0 и Яесг1 = 125600, асг1 = 2,859 для а! = 2, уУ1Ь = 0,4, а вязкой моды II - Яесг11 = 50060, асгп = 2,546 для а1 = уУ1Ь = 0 и Яесг11 = 50060, асгп = 2,546 для а1 = 2, уУ1Ь = 0,4. В качестве критического числа Рейнольдса течения Яесг при фиксированных значениях параметров а1, уУ1Ь принимается минимальное числовое значение из Яесг1 и Яесг11: Яесг = ш1п(Кесг1, Яесг11).

На рис. 7 представлены графики зависимостей Яесг(М) и асг(М) для совершенного (а1 = уУ1Ь = 0) и термически неравновесного (а1 = 2, уУ1Ь = 0,4) газов. Из вида графиков на этом рисунке следует, что на интервале 3 < М < 6,2 кривые Яесг(М) падают, а при 6,2 < М < 15 наоборот возрастают. При М = 6,2 кривые Яесг(М) достигают своих минимальных значений, которые равны Яесг, ш1п = 20169 для а1 = уу1Ь = 0 и Яесг, ш1п = 20473 для а! = 2, уУ1Ь = 0,4. В таблице приведены значения Яесг и асг для различных чисел Маха М и параметров а!, уУ1Ь. Из данных таблицы следует, что критические значения числа Рейнольдса Яесг в среднем по числу Маха М увеличиваются примерно на 12,6 % при возрастании параметров а! и уУ1Ь до их максимальных значений, принятых в расчетах.

Рис. 7. Зависимости Recr(M) (а) и acr(M) (б) (сплошная линия - совершенный газ, штриховая - а! = 2, y„b = 0,4)

Критические числа Рейнольдса Recr и соответствующие им значения волновых чисел acr

Совершенный газ Колебательно-возбужденный газ

M (а1 = Yvib = 0) (а1 = 2, Yvib = 0,4)

Recr acr Recr acr

3 50060 2,5460 50801 2,5570

5 23830 2,1310 24186 2,1830

7 21640 1,9301 21967 1,9934

9 35083 1,8706 35613 1,9275

11 55753 1,8794 56594 1,9246

13 75244 1,8839 76379 1,9240

15 85150 1,8110 86436 1,8650

Заключение

Проведен анализ линейной устойчивости сжимаемого течения Куэтта термически неравновесного молекулярного газа. В результате получено, что при определенном сочетании значений параметров исследуемого течения (чисел Рейнольдса Яе, Маха М, объемной вязкости аь степени колебательной неравновесности уУ1Ь и времени колебательной релаксации т) оно может быть как устойчиво, так и неустойчиво по отношению к малым возмущениям.

Для вязких возмущений рассчитаны спектры собственных значений, инкременты нарастания и кривые нейтральной устойчивости в плоскости (Яе, а) для первой и второй растущих мод в диапазоне чисел М = 2-6 и Яе = 104-107. Найден диапазон изменения критических чисел Рейнольдса Яесг « (2-5)-104. Показано, что при всех уровнях возбуждения наиболее неустойчивой является вторая мода. Возбуждение практически не меняет форму области неустойчивости, но ее границы с ростом возбуждения смещаются в сторону больших волновых чисел. Можно констатировать, что, в общем, возбуждение внутренних степеней свободы молекул газа снижает инкременты нарастания возмущений и оказывает стабилизирующее воздействие на течение.

ЛИТЕРАТУРА

1. Григорьев Ю.Н., Ершов И.В. Энергетическая оценка критических чисел Рейнольдса в сжимаемом течении Куэтта. Влияние объемной вязкости // ПМТФ. 2010. Т. 51. № 5. С. 59-67.

2. Григорьев Ю.Н., Ершов И.В. Критические числа Рейнольдса в течении Куэтта колебательно возбужденного двухатомного газа. Энергетический подход // ПМТФ. 2012. Т. 53. № 4. С. 57-73.

3. Гольдштик М.А., Штерн В.Н. Гидродинамическая устойчивость и турбулентность. Новосибирск: Наука, 1977. 367 с.

4. Duck P.W., Erlebacher G., Hussaini M.Y. On the linear stability of compressible plane Couette flow // J. of Fluid Mech. 1994. V. 258. P 131-165.

5. Романов В.А. Устойчивость плоскопараллельного течения Куэтта // ДАН СССР. 1971. Т. 196. № 5. С. 1049-1051.

6. Hu S., Zhong X. Linear stability of viscous supersonic plane Couette flow // Phys. of Fluids. 1998. V. 10. No. 3. P. 709-729.

7. Ebring A.A. High speed-viscous plane Couette-Poiseuille flow stability. Ph.D. Thesis. Department of Mechanical Engineering. Middle East Technical University. Ankara, Turkey, 2004. 125 p.

8. MalikM., Dey J., Alam M. Linear stability, transient energy growth, and the role of viscosity stratification in compressible plane Couette flow // Physical Rev. E. 2008. V. 77. Issue 3. P. 036322(1)-036322(15).

9. ЛиньЦзя-цзяо. Теория гидродинамической устойчивости. М.: ИЛ, 1958. 195 с.

10. Григорьев Ю.Н., Ершов И.В. Линейная устойчивость невязкого сдвигового течения колебательно возбужденного двухатомного газа // ПММ. 2011. Т. 45. Вып. 4. С. 581-593.

11. Григорьев Ю.Н., Ершов И.В. Устойчивость течений релаксирующих молекулярных газов. Новосибирск: Изд-во СО РАН, 2012. 230 с.

12. Жданов В. М., Алиевский М.Е. Процессы переноса и релаксации в молекулярных газах. М.: Наука, 1989. 336 с.

13. Нагнибеда Е.А., Кустова Е.В. Кинетическая теория процессов переноса и релаксации в потоках неравновесных реагирующих газов. СПб.: Изд-во С.-Петербурского ун-та, 2003. 272 с.

14. Canuto C., Hussaini M.Y., Quarteroni A., Zang T.A. Spectral methods in fluid dynamics: Springer series in Computational Physics. Berlin: Springer-Verlag, 1988. 564 p.

15. Trefethen L.N. Spectral methods in Matlab. Philadelphia: Soc. for Industr. and Appl. Math., 2000. 160 p.

16. Корн Г., Корн Т. Справочник по математике. М.: Наука, 1970. 720 с.

17. Moler C.B., Stewart G.W. An algorithm for generalized matrix eigenvalue problems // SIAM J. Numer. Anal. 1973. V. 10. No. 2. P. 241-256.

Статья поступила 08.10.2014 г.

Ershov I.V. STABILITY OF A SUPERSONIC COUETTE FLOW OF VIBRATIONALLY EXCITED DIATOMIC GAS. DOI 10.17223/19988621/33/5

Within the linear theory, we study stability of the Couette flow of vibrationally excited diatomic gas with a parabolic profile of static temperature. The original mathematical model of the gas flow is the system of equations of two-temperature aerodynamics. As a result, it has been shown that when a certain combination of values of the parameters of the test flow (Reynolds number Re, Mach number M, bulk viscosity a!, the degree of vibrational nonequilibrium yvib, and vibrational relaxation time т), it can be both stable and unstable with respect to small perturbations. For viscous perturbations, the spectra of eigenvalues, the growth increments, and neutral stability curves in the plane (Re, a) were calculated for the first and second growing modes in the range of numbers M = 2-6 and Re = 104-107. The range of variation of the critical Reynolds number Recr « (2-5)-104 was found. It is shown that the second mode is most unstable for all lev-

els of excitation. The excitation does not actually change the shape of the region of instability, but its boundaries shift to higher wave numbers with increasing excitation. It can be stated that, in general, the excitation of internal degrees of freedom of the gas molecules reduces the disturbance growth increments and has a stabilizing effect on the flow.

Keywords: hydrodynamic stability, vibrational relaxation, equations of two-temperature aerodynamics, unstable viscous excitation modes, critical Reynolds number.

ERSHOVIgor Valer'evich (Doctor of Physics and Mathematics, Assoc. Prof., Novosibirsk State University of Architecture and Civil Engineering, Novosibirsk, Russian Federation)

E-mail: i_ershov@ngs.ru

REFERENCES

1. Grigor'ev Yu.N., Ershov I.V. Energeticheskaya otsenka kriticheskikh chisel Reynol'dsa v szhimaemom techenii Kuetta. Vliyanie ob"emnoy vyazkosti. Prikladnaya mekhanika i tekhnicheskaya fizika, 2010, vol. 51, no. 5, pp. 59-67. (in Russian)

2. Grigor'ev Yu.N., Ershov I.V. Kriticheskie chisla Reynol'dsa v techenii Kuetta kolebatel'no vozbuzhdennogo dvukhatomnogo gaza. Energeticheskiy podkhod. Prikladnaya mekhanika i tekhnicheskaya fizika, 2012, vol. 53, no. 4, pp. 57-73. (in Russian)

3. Gol'dshtik M.A., Shtern V.N. Gidrodinamicheskaya ustoychivost' i turbulentnost'. Novosibirsk, Nauka, 1977, 367 p. (in Russian)

4. Duck P.W., Erlebacher G., Hussaini M.Y. On the linear stability of compressible plane Couette flow. J. of FluidMech., 1994, vol. 258, pp. 131-165.

5. Romanov V.A. Ustoychivost' ploskoparallel'nogo techeniya Kuetta. Doklady AN SSSR, 1971, vol. 196, no. 5, pp. 1049-1051. (in Russian)

6. Hu S., Zhong X. Linear stability of viscous supersonic plane Couette flow. Phys. of Fluids, 1998, vol. 10, no. 3, pp. 709-729.

7. Ebrinf A.A. High speed-viscous plane Couette-Poiseuille flow stability. Ph.D. Thesis. Department of Mechanical Engineering. Middle East Technical University. Ankara, Turkey, 2004, 125 p.

8. Malik M., Dey J., Alam M. Linear stability, transient energy growth, and the role of viscosity stratification in compressible plane Couette flow. Physical Rev. E, 2008, vol. 77, issue 3, pp. 036322(1)-036322(15).

9. Lin' Tszya-tszyao. Teoriya gidrodinamicheskoy ustoychivosti. Moskow, IL Publ., 1958, 195 p. (in Russian)

10. Grigor'ev Yu.N., Ershov I.V. Lineynaya ustoychivost' nevyazkogo sdvigovogo techeniya kolebatel'no vozbuzhdennogo dvukhatomnogo gaza. Prikladnaya matematika i mekhanika, 2011, vol. 45. Vyp. 4, pp. 581-593. (in Russian)

11. Grigor'ev Yu.N., Ershov I.V. Ustoychivost' techeniy relaksiruyushchikh molekulyarnykh gazov. Novosibirsk, SO RAN Publ., 2012, 230 p. (in Russian)

12. Zhdanov V.M., Alievskiy M.E. Protsessy perenosa i relaksatsii v molekulyarnykh gazakh. Moskow, Nauka Publ., 1989, 336 p. (in Russian)

13. Nagnibeda E.A., Kustova E.V. Kineticheskaya teoriya protsessov perenosa i relaksatsii v potokakh neravnovesnykh reagiruyushchikh gazov. St. Petersburg, St. Petersburg Univ. Publ., 2003, 272 p. (in Russian)

14. Canuto C., Hussaini M.Y., Quarteroni A., Zang T.A. Spectral methods in fluid dynamics: Springer series in Computational Physics. Berlin, Springer-Verlag, 1988, 564 p.

15. Trefethen L.N. Spectral methods in Matlab. Philadelphia, Soc. for Industr. and Appl. Math., 2000, 160 p.

16. Korn G., Korn T. Spravochnik po matematike. Moskow, Nauka Publ., 1970, 720 p. (in Russian)

17. Moler C.B., Stewart G.W. An algorithm for generalized matrix eigenvalue problems. SIAM J. Numer. Anal., 1973, vol. 10, no. 2, pp. 241-256.