Известия РАН. Механика жидкости и газа, 2019, № 3, стр. 3-15

ЛИНЕЙНАЯ УСТОЙЧИВОСТЬ СВЕРХЗВУКОВОГО ПОГРАНИЧНОГО СЛОЯ РЕЛАКСИРУЮЩЕГО ГАЗА НА ПЛАСТИНЕ

Ю. Н. Григорьев ab*, И. В. Ершов ac**

a Институт вычислительных технологий СО РАН
Новосибирск, Россия

b Новосибирский государственный университет
Новосибирск, Россия

c Новосибирский государственный аграрный университет
Новосибирск, Россия

* E-mail: grigor@ict.nsc.ru
** E-mail: ivershov1969@gmail.com

Поступила в редакцию 16.07.2018
После доработки 18.10.2018
Принята к публикации 18.10.2018

Полный текст (PDF)

Аннотация

На основе линейной теории устойчивости исследовано развитие невязких и вязких двумерных дозвуковых возмущений в сверхзвуковом пограничном слое колебательно возбужденного газа на плоской пластине. В качестве исходной модели использовалась система двухтемпературной газовой динамики, включающая релаксационное уравнение Ландау–Теллера. Невозмущенный поток описывался автомодельным погранслойным решением для совершенного газа. Показано, что в невязком приближении возбуждение снижает максимальные инкременты нарастания наиболее неустойчивой второй моды на 10–12% по сравнению с идеальным газом. Рассчитаны кривые нейтральной устойчивости для первой и второй наиболее неустойчивых мод для чисел Маха М = 2.2, 4.5, 4.8. Для обеих мод критические числа Рейнольдса при максимальном возбуждении превышают на 12–13% соответствующие значения для совершенного газа.

Ключевые слова: линейная устойчивость, неустойчивые моды, критические числа Рейнольдса

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

В этом направлении наиболее систематические экспериментальные результаты, подтвержденные адекватным численным моделированием, были получены группой профессора Г. Хорнунга [4, 5, 7, 8]. Для обтекания острых конусов гиперскоростными потоками воздуха, азота N2, кислорода O2 и углекислого газа CO2 при полных энтальпиях потока H0 = 3–15 Дж/кг было, в частности, показано, что при одинаковых условиях число Рейнольдса ламинарно-турбулентного перехода в CO2 примерно в четыре раза больше по сравнению с двухатомными N2, O2 и воздухом. Это объясняется тем, что частота второй акустической моды, через которую происходит ламинарно-турбулентный переход в гиперзвуковых пристенных течениях [9], близка к релаксационной частоте легко возбудимой изгибной колебательной моды CO2, в то время как частоты колебательных мод двухатомных молекул существенно выше.

В серии статей авторов, результаты которых представлены в монографии [10], для ряда модельных задач, в том числе для сверхзвукового течения Куэтта и свободного сдвигового слоя на основе линейной и энергетической теорий устойчивости было показано, что релаксация колебательных степеней свободы молекул приводит к уменьшению инкрементов нарастания возмущений и повышению критических чисел Рейнольдса в пределах 10–15%.

Влияние релаксации внутренних степеней свободы молекул на устойчивость пограничного слоя на плоской пластине рассматривалось в [6, 11, 12]. Однако из-за неадекватной постановки задачи полученные результаты вызывают серьезную критику. Так, в [11, 12] и ряде цитированных в них источников используется нефизическое понятие знакопеременной “второй вязкости”. Анализ его противоречивости выходит далеко за рамки этой статьи. Но чтобы подчеркнуть нефизичность этой характеристики, достаточно сказать, что в [11, 12] в исходных уравнениях релаксационной газовой динамики (по сравнению ниже с системой (1.1)–(1.4)) опущено дивергентное слагаемое с коэффициентом объемной вязкости. Последний по современным представлениям кинетической теории газов [13] реально описывает диссипативный процесс, возникающий при слабой термической неравновесности молекулярных газов, положительно определен и имеет прямую связь с введенным Стоксом коэффициентом второй вязкости.

В обеих статьях [11, 12] в релаксационном уравнении Ландау–Теллера отсутствует диффузионный перенос колебательной энергии, который необходимо учитывать в пристенных течениях в силу существенных поперечных градиентов колебательной температуры. В то же время в это уравнение вводится источник энергии (накачку), с тем, чтобы обеспечить стационарный процесс релаксации. Однако при этом в постановке задачи отсутствует процесс отвода энергии, который мог бы описываться эндотермической реакцией в объеме или конвективным теплообменом на поверхности. Тем не менее в уравнении энергии нет соответствующего слагаемого, а на стенке ставится условие адиабатичности.

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

В [6] на основе численного моделирования рассматривалась задача линейной устойчивости пограничного слоя при обтекании плоской пластины термически неравновесным воздухом в полной постановке, которая используется в данной работе. Рассматривались условия реального полета на высоте H = 12 км при числе Маха M = 4.5, когда релаксация состоит в возбуждении колебательных мод, и его моделирование в аэродинамической трубе, где, наоборот, энергия колебательных мод, “замороженных” при температуре торможения, в процессе релаксации повышает статическую температуру потока. В статье получены неоднозначные результаты воздействия релаксации на устойчивость течения. В случае, когда в рассматриваемых условиях учитывалось возбуждение колебательных мод, которое в данных условиях не могло быть значительным, отмечено сильное дестабилизирующее влияние на первую моду с волновым вектором β = π/3 относительно несущего потока. При этом, как следует из расчетов, первая мода с этим волновым вектором оказалась наиболее неустойчивой и для совершенного газа. Это противоречит общепризнанным результатам [9, 16], согласно которым для M > 4 наиболее неустойчивой считается вторая акустическая мода с волновым вектором β = 0. Для этой моды получено лишь слабое стабилизирующее влияние для пограничного слоя на пластине с острой передней кромкой.

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

1. ОСНОВНЫЕ УРАВНЕНИЯ И ПОСТАНОВКА ЗАДАЧИ

В рамках линейной теории устойчивости рассматривается развитие двумерных возмущений при обтекании колебательно возбужденным газом плоской полубесконечной пластины. Начало декартовой системы координат (x, y) совпадает с носиком пластины, координата х ориентирована вдоль пластины по направлению невозмущенного потока, координата у соответственно направлена в поток по нормали к пластине. Течение описывается системой уравнений двухтемпературной газовой динамики [10, 18].

В качестве характерных величин для обезразмеривания выбраны текущее расстояние x = L вдоль пластины, параметры невозмущенного потока вне пограничного слоя, отмеченные индексом “∞”, – скорость U, плотность ρ и температура T, коэффициенты сдвиговой η и объемной ηb вязкостей, коэффициент теплопроводности, обусловленный переносом энергии в поступательных и вращательных степенях свободы ${{k}_{\infty }} = {{k}_{{t\infty }}} + {{k}_{{r\infty }}}$, коэффициент теплопроводности, описывающий диффузионный перенос энергии колебательных квантов, kv. Для обезразмеривания давления и времени используются комбинированные величины ${{{\rho }}_{\infty }}U_{\infty }^{2}$ и L/U соответственно.

В обезразмеренных таким образом переменных исходная система уравнений имеет вид

(1.1)
$\frac{{\partial {\rho }}}{{\partial t}} + {{u}_{i}}\frac{{\partial {\rho }}}{{\partial {{x}_{i}}}} + {\rho }\frac{{\partial {{u}_{i}}}}{{\partial {{x}_{i}}}} = 0$
(1.2)
${\rho }\frac{{\partial {{u}_{i}}}}{{\partial t}} + {\rho }{{u}_{j}}\frac{{\partial {{u}_{i}}}}{{\partial {{x}_{j}}}} = - \frac{{\partial p}}{{\partial {{x}_{i}}}} + \frac{1}{{\operatorname{Re} }}\frac{{\partial {{{\sigma }}_{{ij}}}}}{{\partial {{x}_{j}}}}$
(1.3)
$\begin{gathered} {\rho }\frac{{\partial T}}{{\partial t}} + {\rho }{{u}_{i}}\frac{{\partial T}}{{\partial {{x}_{i}}}} + {\gamma }({\gamma } - 1){{{\text{M}}}^{2}}p\frac{{\partial {{u}_{i}}}}{{\partial {{x}_{i}}}} = \\ = \;\frac{{\gamma }}{{\operatorname{Re} \Pr }}\frac{\partial }{{\partial {{x}_{i}}}}\left( {k\frac{{\partial T}}{{\partial {{x}_{i}}}}} \right) + \frac{{{\gamma }({\gamma } - 1){{{\text{M}}}^{2}}}}{{\operatorname{Re} }}{{{\sigma }}_{{ij}}}\frac{{\partial {{u}_{i}}}}{{\partial {{x}_{j}}}} - \frac{{{{{\gamma }}_{\text{v}}}{\rho }(T - {{T}_{\text{v}}})}}{{\tau }} \\ \end{gathered} $
(1.4)
${{{\gamma }}_{\text{v}}}\left[ {{\rho }\frac{{\partial {{T}_{\text{v}}}}}{{\partial t}} + {\rho }{{u}_{i}}\frac{{\partial {{T}_{\text{v}}}}}{{\partial {{x}_{i}}}}} \right] = \frac{{{\gamma }{{{\mu }}_{2}}}}{{\operatorname{Re} \Pr }}\frac{\partial }{{\partial {{x}_{i}}}}\left( {k\text{v}\frac{{\partial {{T}_{\text{v}}}}}{{\partial {{x}_{i}}}}} \right) + \frac{{{{{\gamma }}_{\text{v}}}{\rho }(T - {{T}_{\text{v}}})}}{{\tau }}$

Здесь тензор вязких напряжений выражается как

(1.5)
${{{\sigma }}_{{ij}}} = {\eta }(T)\left( {\frac{{\partial {{u}_{i}}}}{{\partial {{x}_{j}}}} + \frac{{\partial {{u}_{j}}}}{{\partial {{x}_{i}}}} - \frac{2}{3}\frac{{\partial {{u}_{k}}}}{{\partial {{x}_{k}}}}{{\delta }_{{ij}}}} \right) + {{{\eta }}_{b}}(T)\frac{{\partial {{u}_{k}}}}{{\partial {{x}_{k}}}}{{{\delta }}_{{ij}}}$
$i,j = 1,\;2,\quad {{x}_{1}} = x,\quad {{x}_{2}} = y$

Продольная и поперечная составляющие скорости задаются равенствами u1 = u, u2 = v. Система (1.1)–(1.5) замыкается уравнением состояния идеального газа

(1.6)
${\gamma }{{{\text{M}}}^{2}}p = {\rho }T$

Коэффициенты и параметры, входящие в уравнения (1.1)(1.6), определяются следующим образом. Для температурной зависимости коэффициента сдвиговой вязкости используется формула Сазерленда [19], температурная зависимость коэффициента объемной вязкости описывается формулой ${{{\eta }}_{b}}(T) = {{{\mu }}_{1}}{\eta }(T)$, где ${{{\mu }}_{1}} = {{{{{\eta }}_{{b\infty }}}} \mathord{\left/ {\vphantom {{{{{\eta }}_{{b\infty }}}} {{{{\eta }}_{\infty }}}}} \right. \kern-0em} {{{{\eta }}_{\infty }}}}$. Для двухатомных газов значения параметра μ1 лежат в диапазоне: μ1 = 0–2 [10]. Коэффициенты теплопроводности в (1.3) и (1.4) выражены через коэффициент сдвиговой вязкости с помощью полуэмпирических соотношений Эйкена [10, 13]

(1.7)
$k(T) = {{k}_{t}}(T) + {{k}_{r}}(T),\quad {{k}_{t}}(T) = \frac{5}{2}{{c}_{{Vt}}}{\eta }(T),\quad {{k}_{r}}(T) = \frac{6}{5}{{c}_{{Vr}}}{\eta }(T),\quad {{k}_{\text{v}}}(T) = \frac{6}{5}{{c}_{{V\text{v}}}}{\eta }(T)$
где cVt, cVr, cVv – удельные теплоемкости при постоянных объеме, обусловленные поступательным, вращательным и колебательным движениями молекул соответственно. Предполагается, что коэффициенты теплоемкости постоянны. Принимается также, что поступательные и вращательные степени свободы молекул находятся в состоянии квазиравновесия. При этом для их энергии имеет место равнораспределение по степеням свободы, и соответствующие теплоемкости определяются соотношениями: cVt, = 3R/2 и cVr, = R, где R – газовая постоянная.

Коэффициент γ = cp,/cV, – показатель адиабаты, cV, = cVt, + cVr,, cp, = cV, + R – соответственно удельные теплоемкости при постоянных объеме и давлении при температуре T. Коэффициент ${{{\gamma }}_{v}} = {{c}_{{Vv,\infty }}}{\text{/}}({{c}_{{Vt,\infty }}} + {{c}_{{Vr,\infty }}}) = {{c}_{{Vv,\infty }}}{\text{/}}{{c}_{{V,\infty }}}$ характеризует степень неравновесности колебательной степени свободы, τ – характерное время релаксации последней. Критерии Re = = ${{{\rho }}_{\infty }}L{{U}_{\infty }}{\text{/}}{{{\eta }}_{\infty }}$, ${\text{M}} = {{U}_{\infty }}{\text{/}}\sqrt {{\gamma }R{{T}_{\infty }}} $и $\Pr = {{c}_{{V,\infty }}}{{{\eta }}_{\infty }}{\text{/}}{{k}_{\infty }}$ – соответственно числа Рейнольдса, Маха и Прандтля невозмущенного потока. Параметр ${{{\mu }}_{2}} = {{k}_{{\text{v},\infty }}}{\text{/}}({{k}_{{t,\infty }}} + {{k}_{{r,\infty }}})$. С использованием соотношений (1.7) параметр μ2 выражается в виде

${{{\mu }}_{2}} = \frac{{{{k}_{{\text{v},\infty }}}}}{{{{k}_{{t,\infty }}} + {{k}_{{r,\infty }}}}} = \frac{{12{{\gamma }_{\text{v}}}({{c}_{{Vt,\infty }}} + {{c}_{{Vr,\infty }}})}}{{25{{c}_{{Vt,\infty }}} + 12{{c}_{{Vr,\infty }}}}} = \frac{{20{{{\gamma }}_{\text{v}}}}}{{33}}$

Система (1.1)–(1.6) линеаризуется на стационарном автомодельном погранслойном решении в отсутствие колебательного возбуждения [20]. Чтобы показать, что оно не зависит от объемной вязкости (параметра μ1), запишем стационарную систему уравнений динамики совершенного газа в погранслойном приближении. Для этого введем дополнительное “погранслойное” масштабирование поперечной координаты и скорости

$y{\text{'}} = y\sqrt {\operatorname{Re} } ,\quad \text{v}{\text{'}} = \text{v}\sqrt {\operatorname{Re} } $

Все остальные используемые ниже штрихованные переменные сохраняют обезразмеривание системы (1.1)–(1.6). В результате система (1.1)–(1.6) примет вид

$\frac{{\partial {\rho '}u{\text{'}}}}{{\partial x{\text{'}}}} + \frac{{\partial {\rho '}\text{v}{\text{'}}}}{{\partial y{\text{'}}}} = 0$
$\begin{gathered} {\rho '}\left( {u{\text{'}}\frac{{\partial u{\text{'}}}}{{\partial x{\text{'}}}} + \text{v}{\text{'}}\frac{{\partial u{\text{'}}}}{{\partial y{\text{'}}}}} \right) = - \frac{{\partial p{\text{'}}}}{{\partial x{\text{'}}}} + \frac{\partial }{{\partial y{\text{'}}}}\left( {{\eta '}\frac{{\partial u{\text{'}}}}{{\partial y{\text{'}}}}} \right) + \\ + \;\frac{1}{{\operatorname{Re} }}\frac{\partial }{{\partial x{\text{'}}}}\left[ {2{\eta '}\frac{{\partial u{\text{'}}}}{{\partial x{\text{'}}}} + \left( {{\eta }_{b}^{'} - \frac{{2{\eta '}}}{3}} \right)\left( {\frac{{\partial u{\text{'}}}}{{\partial x{\text{'}}}} + \frac{{\partial \text{v}{\text{'}}}}{{\partial y{\text{'}}}}} \right)} \right] + \frac{1}{{\operatorname{Re} }}\frac{\partial }{{\partial y{\text{'}}}}\left( {{\eta '}\frac{{\partial \text{v}{\text{'}}}}{{\partial x{\text{'}}}}} \right) \\ \end{gathered} $
$\begin{gathered} \frac{{{\rho '}}}{{\operatorname{Re} }}\left( {u{\text{'}}\frac{{\partial \text{v}{\text{'}}}}{{\partial x{\text{'}}}} + \text{v}{\text{'}}\frac{{\partial \text{v}{\text{'}}}}{{\partial y{\text{'}}}}} \right) = - \frac{{\partial p{\text{'}}}}{{\partial y{\text{'}}}} + \frac{1}{{\operatorname{Re} }}\frac{\partial }{{\partial x{\text{'}}}}\left[ {{\eta '}\left( {\frac{{\partial u{\text{'}}}}{{\partial y{\text{'}}}} + \frac{1}{{\operatorname{Re} }}\frac{{\partial \text{v}{\text{'}}}}{{\partial x{\text{'}}}}} \right)} \right] + \\ + \;\frac{1}{{\operatorname{Re} }}\frac{\partial }{{\partial y{\text{'}}}}\left[ {2{\eta '}\frac{{\partial \text{v}{\text{'}}}}{{\partial y{\text{'}}}} + \left( {{\eta }_{b}^{'} - \frac{{2{\eta '}}}{3}} \right)\left( {\frac{{\partial u{\text{'}}}}{{\partial x{\text{'}}}} + \frac{{\partial \text{v}{\text{'}}}}{{\partial y{\text{'}}}}} \right)} \right] \\ \end{gathered} $
$\begin{gathered} {\rho '}\left( {u{\text{'}}\frac{{\partial T{\text{'}}}}{{\partial x{\text{'}}}} + \text{v}{\text{'}}\frac{{\partial T{\text{'}}}}{{\partial y{\text{'}}}}} \right) + {\gamma }({\gamma } - 1){{{\text{M}}}^{2}}p{\text{'}}\left( {\frac{{\partial u{\text{'}}}}{{\partial x{\text{'}}}} + \frac{{\partial \text{v}{\text{'}}}}{{\partial y{\text{'}}}}} \right) = \\ = \;\frac{{\gamma }}{{\Pr }}\left[ {\frac{1}{{\operatorname{Re} }}\frac{\partial }{{\partial x{\text{'}}}}\left( {k{\text{'}}\frac{{\partial T{\text{'}}}}{{\partial x{\text{'}}}}} \right) + \frac{\partial }{{\partial y{\text{'}}}}\left( {k{\text{'}}\frac{{\partial T{\text{'}}}}{{\partial y{\text{'}}}}} \right)} \right] + {\gamma }({\gamma } - 1){{{\text{M}}}^{2}}{\eta '}{{\left( {\frac{{\partial u{\text{'}}}}{{\partial y{\text{'}}}} + \frac{1}{{\operatorname{Re} }}\frac{{\partial \text{v}{\text{'}}}}{{\partial x{\text{'}}}}} \right)}^{2}} + \\ + \;\frac{{2{\gamma }({\gamma } - 1){{{\text{M}}}^{2}}{\eta '}}}{{\operatorname{Re} }}\left[ {{{{\left( {\frac{{\partial u{\text{'}}}}{{\partial x{\text{'}}}}} \right)}}^{2}} + {{{\left( {\frac{{\partial \text{v}{\text{'}}}}{{\partial y{\text{'}}}}} \right)}}^{2}}} \right] + \frac{{{\gamma }({\gamma } - 1){{{\text{M}}}^{2}}}}{{\operatorname{Re} }}\left( {{\eta }_{b}^{'} - \frac{{{2\eta '}}}{3}} \right){{\left( {\frac{{\partial u{\text{'}}}}{{\partial x{\text{'}}}} + \frac{{\partial \text{v}{\text{'}}}}{{\partial y{\text{'}}}}} \right)}^{2}} \\ \end{gathered} $
${\gamma }{{{\text{M}}}^{2}}p{\text{'}} = {\rho '}T{\text{'}}$

Используя разложение по малому параметру

$f = {{f}^{{(0)}}} + \frac{1}{{\sqrt {\operatorname{Re} } }}{{f}^{{(1)}}} + \frac{1}{{\operatorname{Re} }}{{f}^{{(2)}}} + \;...$
где f – вектор текущих параметров газа, в нулевом приближении при $1{\text{/}}\sqrt {\operatorname{Re} } \to 0$ получаем классическую систему уравнений пограничного слоя Прандтля, автомодельное решение которой используется при линеаризации. Очевидно, что объемная вязкость в нее не входит. Выписывая последующие приближения, легко убедиться в том, что влияние объемной вязкости проявляется только во втором приближении по параметру $1{\text{/}}\sqrt {\operatorname{Re} } $, которое, в частности, используется при необходимости учесть непараллельность течения в пограничном слое. Именно в этом случае становится заметной расходимость потока, приводящая к проявлению объемной деформации, с которой связана объемная вязкость.

Автомодельные профили скорости и температуры рассчитывались численным интегрированием системы уравнений в автомодельных переменных

(1.8)
$\frac{d}{{d\xi }}\left( {C\frac{{{{d}^{2}}f}}{{d{{\xi }^{2}}}}} \right) + f\frac{{{{d}^{2}}f}}{{d{{\xi }^{2}}}} = 0$
(1.9)
$\frac{d}{{d\xi }}\left( {k{\text{'}}\frac{{dT{\text{'}}}}{{d\xi }}} \right) + \frac{{(\xi - 1)}}{4}\Pr {{{\text{M}}}^{2}}C{{\left( {\frac{{{{d}^{2}}f}}{{d{{\xi }^{2}}}}} \right)}^{2}} + \Pr f\frac{{dT{\text{'}}}}{{d\xi }} = 0$
с граничными условиями

${{\left[ {\frac{{df}}{{d{\xi }}}} \right]}_{{{\xi } = {\text{0}}}}} = 0,\quad {{\left[ {\frac{{df}}{{d{\xi }}}} \right]}_{{{\xi } = {{{\delta }}_{{\text{c}}}}}}} = 2,\quad {{\left[ {\frac{{dT{\text{'}}}}{{d{\xi }}}} \right]}_{{{\xi } = {\text{0}}}}} = 0,\quad T{\text{'}}({{{\delta }}_{{\text{c}}}}) = 1$

Здесь автомодельная переменная Дородницына–Хоуарта [20]

(1.10)
${\xi } = \frac{{\text{1}}}{{\text{2}}}\int\limits_0^{y{\text{'}}} {\rho } {\text{'}}dy{\text{'}}$
$\frac{{df}}{{d{\xi }}} = 2u({\xi }),\quad C = {\rho '}(T{\text{'}}){\eta '}(T{\text{'}}) = \frac{{{\eta '}(T{\text{'}})}}{{T{\text{'}}}}$
где δс – условная верхняя граница пограничного слоя

$\left| {\frac{{{{U}_{\infty }} - u({{{\delta }}_{{\text{c}}}})}}{{{{U}_{\infty }}}}} \right| \leqslant {{10}^{{ - 6}}}$

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

(1.11)
$u = {{U}_{s}} + u*,\quad \text{v} = \text{v}*,\quad {\rho } = {{{\rho }}_{s}} + {\rho *},\quad T = {{T}_{s}} + T*,\quad {{T}_{\text{v}}} = {{T}_{{\text{v},s}}} + T_{\text{v}}^{*},\quad p = {{p}_{s}} + p{\text{*}}$

Учитывались температурные возмущения коэффициентов переноса

${\eta }(T) = {\eta }({{T}_{s}}) + {{{\eta }}_{{T,s}}}T* = {{{\eta }}_{s}} + {\eta *},\quad k(T) = k({{T}_{s}}) + {{k}_{{T,s}}}T* = {{k}_{s}} + k{\text{*}}$
${{k}_{\text{v}}}(T) = {{k}_{\text{v}}}({{T}_{s}}) + {{k}_{{\text{v},Ts}}}T* = {{k}_{{\text{v},s}}} + k_{\text{v}}^{*}$
${{{\eta }}_{{T,s}}} = {{\left[ {\frac{{d{\eta }}}{{dT}}} \right]}_{{T = {{T}_{s}}}}},\quad {{k}_{{T,s}}} = {{\left[ {\frac{{dk}}{{dT}}} \right]}_{{T = {{T}_{s}}}}},\quad {{k}_{{\text{v}T,s}}} = {{\left[ {\frac{{d{{k}_{\text{v}}}}}{{dT}}} \right]}_{{T = {{T}_{s}}}}}$

Здесь индексом s обозначены значения гидродинамических переменных, относящиеся к сдвиговому течению, а “*” – возмущения этих переменных.

Подстановка (1.11) в исходную систему (1.1)–(1.6) и линеаризация ее с точностью до членов первого порядка малости по возмущениям гидродинамических переменных приводит к следующей системе уравнений для возмущений:

$\frac{{\partial {\rho *}}}{{\partial t}} + {{U}_{s}}\frac{{\partial {\rho *}}}{{\partial x}} + {{{\rho }}_{s}}\left( {\frac{{\partial u{\text{*}}}}{{\partial x}} + \frac{{\partial \text{v}{\text{*}}}}{{\partial y}}} \right) + \text{v}{\text{*}}\frac{{\partial {{{\rho }}_{s}}}}{{\partial y}} = 0,$
$\begin{gathered} {{\rho }_{s}}\frac{{\partial u{\text{*}}}}{{\partial t}} + {{\rho }_{s}}{{U}_{s}}\frac{{\partial u{\text{*}}}}{{\partial x}} + {{\rho }_{s}}\text{v}{\text{*}}\frac{{\partial {{U}_{s}}}}{{\partial y}} = - \frac{{\partial p{\text{*}}}}{{\partial x}} + \frac{{{{\eta }_{s}}}}{{\operatorname{Re} }}\left( {\frac{{{{\partial }^{2}}u{\text{*}}}}{{\partial {{x}^{2}}}} + \frac{{{{\partial }^{2}}u{\text{*}}}}{{\partial {{y}^{2}}}}} \right) + \\ + \;\frac{{{{\eta }_{s}}}}{{\operatorname{Re} }}\left( {{{\mu }_{1}} + \frac{1}{3}} \right)\left( {\frac{{{{\partial }^{2}}u{\text{*}}}}{{\partial {{x}^{2}}}} + \frac{{{{\partial }^{2}}\text{v}{\text{*}}}}{{\partial x\partial y}}} \right) + \frac{1}{{\operatorname{Re} }}\left( {\frac{{\partial u{\text{*}}}}{{\partial y}} + \frac{{\partial \text{v}{\text{*}}}}{{\partial x}}} \right)\frac{{d{{\eta }_{s}}}}{{dy}} + \frac{{\eta {\text{*}}}}{{\operatorname{Re} }}\frac{{{{d}^{2}}{{U}_{s}}}}{{d{{y}^{2}}}} + \frac{1}{{\operatorname{Re} }}\frac{{\partial \eta {\text{*}}}}{{\partial y}}\frac{{d{{U}_{s}}}}{{dy}} \\ \end{gathered} $
$\begin{gathered} {{\rho }_{s}}\frac{{\partial \text{v}{\text{*}}}}{{\partial t}} + {{\rho }_{s}}{{U}_{s}}\frac{{\partial \text{v}{\text{*}}}}{{\partial x}} = - \frac{{\partial p{\text{*}}}}{{\partial y}} + \frac{{{{\eta }_{s}}}}{{\operatorname{Re} }}\left( {\frac{{{{\partial }^{2}}\text{v}{\text{*}}}}{{\partial {{x}^{2}}}} + \frac{{{{\partial }^{2}}\text{v}{\text{*}}}}{{\partial {{y}^{2}}}}} \right) + \frac{{{{\eta }_{s}}}}{{\operatorname{Re} }}\left( {{{\mu }_{1}} + \frac{1}{3}} \right)\left( {\frac{{{{\partial }^{2}}u{\text{*}}}}{{\partial y\partial x}} + \frac{{{{\partial }^{2}}\text{v}{\text{*}}}}{{\partial {{y}^{2}}}}} \right) + \\ + \;\frac{1}{{\operatorname{Re} }}\left[ {\left( {{{\mu }_{1}} - \frac{2}{3}} \right)\frac{{\partial u{\text{*}}}}{{\partial x}} + \left( {{{\mu }_{1}} + \frac{4}{3}} \right)\frac{{\partial \text{v}{\text{*}}}}{{\partial y}}} \right]\frac{{d{{\eta }_{s}}}}{{dy}} + \frac{1}{{\operatorname{Re} }}\frac{{d{{U}_{s}}}}{{dy}}\frac{{\partial \eta {\text{*}}}}{{\partial x}} \\ \end{gathered} $
$\begin{gathered} {{\rho }_{s}}\frac{{\partial T{\text{*}}}}{{\partial t}} + {{\rho }_{s}}{{U}_{s}}\frac{{\partial T{\text{*}}}}{{\partial x}} + {{\rho }_{s}}\text{v}{\text{*}}\frac{{d{{T}_{s}}}}{{dy}} + \gamma {\text{(}}\gamma - 1){{{\text{M}}}^{2}}{{p}_{s}}\left( {\frac{{\partial u{\text{*}}}}{{\partial x}} + \frac{{\partial \text{v}{\text{*}}}}{{\partial y}}} \right) = \\ = \;\frac{{\gamma {{k}_{s}}}}{{\operatorname{Re} \Pr }}\left( {\frac{{{{\partial }^{2}}T{\text{*}}}}{{\partial {{x}^{2}}}} + \frac{{{{\partial }^{2}}T{\text{*}}}}{{\partial {{y}^{2}}}}} \right) + \frac{\gamma }{{\operatorname{Re} \Pr }}\left[ {\frac{{\partial T{\text{*}}}}{{\partial y}}\left( {\frac{{d{{k}_{s}}}}{{dy}}} \right) + \frac{{\partial k{\text{*}}}}{{\partial y}}\left( {\frac{{d{{T}_{s}}}}{{dy}}} \right) + k{\text{*}}\frac{{{{d}^{2}}{{T}_{s}}}}{{d{{y}^{2}}}}} \right] + \\ + \;\frac{{2{{\eta }_{s}}\gamma {\text{(}}\gamma - 1){{{\text{M}}}^{2}}}}{{\operatorname{Re} }}\left( {\frac{{\partial u{\text{*}}}}{{\partial y}} + \frac{{\partial \text{v}{\text{*}}}}{{\partial x}}} \right)\frac{{d{{U}_{s}}}}{{dy}} + \frac{{\gamma {\text{(}}\gamma - 1){{{\text{M}}}^{2}}}}{{\operatorname{Re} }}{{\left( {\frac{{d{{U}_{s}}}}{{dy}}} \right)}^{2}}\eta {\text{*}} - \frac{{{{\gamma }_{\text{v}}}{{\rho }_{s}}(T{\text{*}} - T_{\text{v}}^{*})}}{\tau } \\ \end{gathered} $
$\begin{gathered} {{\rho }_{s}}\frac{{\partial T_{\text{v}}^{*}}}{{\partial t}} + {{\rho }_{s}}{{U}_{s}}\frac{{\partial T_{\text{v}}^{*}}}{{\partial x}} + {{\rho }_{s}}\text{v}{\text{*}}\frac{{d{{T}_{{\text{v},s}}}}}{{dy}} = \frac{{\gamma {{\mu }_{2}}{{k}_{{\text{v},s}}}}}{{{{\gamma }_{\text{v}}}\operatorname{Re} \Pr }}\left( {\frac{{{{\partial }^{2}}T_{\text{v}}^{*}}}{{\partial {{x}^{2}}}} + \frac{{{{\partial }^{2}}T_{\text{v}}^{*}}}{{\partial {{y}^{2}}}}} \right) + \\ + \;\frac{{\gamma {{\mu }_{2}}}}{{{{\gamma }_{\text{v}}}\operatorname{Re} \Pr }}\left[ {\frac{{\partial T_{\text{v}}^{*}}}{{\partial y}}\frac{{d{{k}_{{\text{v},s}}}}}{{dy}} + \frac{{\partial k_{\text{v}}^{*}}}{{\partial y}}\frac{{d{{T}_{{\text{v},s}}}}}{{dy}} + k_{\text{v}}^{*}\frac{{{{d}^{2}}{{T}_{{\text{v},s}}}}}{{d{{y}^{2}}}}} \right] + \frac{{{{\rho }_{s}}(T{\text{*}} - T_{\text{v}}^{*})}}{\tau } \\ \end{gathered} $
$\gamma {{{\text{M}}}^{2}}p{\text{*}} = {{\rho }_{s}}T{\text{*}} + {{T}_{s}}\rho {\text{*}}$

Принималось, что при y = 0 и y = δc все возмущения обращаются в нуль.

Рассматривалась устойчивость периодических по продольной координате x возмущений в форме бегущих плоских волн

${\mathbf{q}}{\text{'}}(x,y,t) = {\mathbf{q}}(y)\exp [i{\alpha }(x - ct)]$
(1.12)
${\mathbf{q}}{\text{*}}(x,y,t) = (u*,\text{v}*,{\rho *},T*,T_{\text{v}}^{*},p*),\quad {\mathbf{q}}(y) = (u,{\alpha }\text{v},{\rho },{\theta },{{{\theta }}_{\text{v}}},p)$
где α – волновое число вдоль периодической переменной x, c = cr + ici – комплексная фазовая скорость, i – мнимая единица.

Подстановка (1.12) в систему уравнений для возмущений дает систему уравнений для их амплитуд

(1.13)
$D{\rho } + \alpha {{\text{v}}_{s}}{\sigma } + {\alpha }\text{v}\frac{{d{{{\rho }}_{s}}}}{{dy}} = 0$
(1.14)
$\begin{gathered} \frac{{{{{\eta }}_{s}}}}{{\operatorname{Re} }}\Delta u - {{{\rho }}_{s}}Du - {\alpha }{{{\rho }}_{s}}\text{v}\frac{{d{{U}_{s}}}}{{dy}} - i{\alpha \varepsilon } + \\ + \;\frac{{{{{\eta }}_{{T,s}}}}}{{\operatorname{Re} }}\frac{{d{{T}_{s}}}}{{dy}}\left( {\frac{{du}}{{dy}} + i{{{\alpha }}^{2}}\text{v}} \right) + \frac{{\theta }}{{\operatorname{Re} }}\frac{d}{{dy}}\left( {{{{\eta }}_{{T,s}}}\frac{{d{{U}_{s}}}}{{dy}}} \right) + \frac{{{{{\eta }}_{{T,s}}}}}{{\operatorname{Re} }}\frac{{d{{U}_{s}}}}{{dy}}\frac{{d{\theta }}}{{dy}} = 0 \\ \end{gathered} $
(1.15)
$\frac{{{\alpha }{{{\eta }}_{s}}}}{{\operatorname{Re} }}\Delta \text{v} - \alpha {{\rho }_{s}}Dv - \frac{{d{\varepsilon }}}{{dy}} + \frac{{{\alpha }{{{\eta }}_{{T,s}}}}}{{\operatorname{Re} }}\frac{{d{{T}_{s}}}}{{dy}}\left( {\frac{{dv}}{{dy}} - iu} \right) + \frac{{i{\alpha }{{{\eta }}_{{T,s}}}}}{{\operatorname{Re} }}\frac{{d{{U}_{s}}}}{{dy}}{\theta } = 0$
(1.16)
$\begin{gathered} \frac{{{\gamma }{{k}_{s}}}}{{\operatorname{Re} \Pr }}\Delta {\theta } - {{{\rho }}_{s}}D\theta - {\alpha }{{{\rho }}_{s}}\frac{{d{{T}_{s}}}}{{dy}}\text{v} - {\alpha (\gamma } - 1){\sigma } + \\ + \;\frac{{2{\gamma (\gamma } - 1){{{\eta }}_{s}}{{{\text{M}}}^{2}}}}{{\operatorname{Re} }}\frac{{d{{U}_{s}}}}{{dy}}\left( {\frac{{du}}{{dy}} + i{{{\alpha }}^{2}}\text{v}} \right) + \left[ {\frac{{\gamma }}{{\operatorname{Re} \Pr }}\frac{d}{{dy}}\left( {{{k}_{{T,s}}}\frac{{d{{T}_{s}}}}{{dy}}} \right) + \frac{{{\gamma (\gamma } - 1){{{\text{M}}}^{2}}{{{\eta }}_{{T,s}}}}}{{\operatorname{Re} }}{{{\left( {\frac{{d{{U}_{s}}}}{{dy}}} \right)}}^{2}}} \right]{\theta } - \\ - \;\frac{{{{{\gamma }}_{\text{v}}}{{{\rho }}_{s}}({\theta } - {{{\theta }}_{\text{v}}})}}{{\tau }} + \frac{{2{\gamma }{{k}_{{T,s}}}}}{{\operatorname{Re} \Pr }}\frac{{d{{T}_{s}}}}{{dy}}\frac{{d{\theta }}}{{dy}} = 0 \\ \end{gathered} $
(1.17)
$\begin{gathered} \frac{{\gamma {{\mu }_{2}}{{k}_{{\text{v},s}}}}}{{\operatorname{Re} \Pr }}\Delta {{\theta }_{\text{v}}} - {{\gamma }_{\text{v}}}{{\rho }_{s}}D{{\theta }_{\text{v}}} - \alpha {{\gamma }_{\text{v}}}{{\rho }_{s}}\frac{{d{{T}_{s}}}}{{dy}}\text{v} + \frac{{{{\gamma }_{\text{v}}}{{\rho }_{s}}(\theta - {{\theta }_{\text{v}}})}}{\tau } + \\ + \;\frac{{\gamma {{\mu }_{2}}}}{{\operatorname{Re} \Pr }}\frac{d}{{dy}}\left( {{{k}_{{\text{v}T,s}}}\frac{{d{{T}_{s}}}}{{dy}}} \right)\theta + \frac{{\gamma {{\mu }_{2}}{{k}_{{\text{v}T,s}}}}}{{\operatorname{Re} \Pr }}\frac{{d{{T}_{s}}}}{{dy}}\left( {\frac{{d\theta }}{{dy}} + \frac{{d{{\theta }_{\text{v}}}}}{{dy}}} \right) = 0 \\ \end{gathered} $
(1.18)
$\gamma {{{\text{M}}}^{2}}p = {{\rho }_{s}}\theta + {{T}_{s}}\rho $
$D = i\alpha ({{U}_{s}} - c),\quad \sigma = iu + \frac{{d\text{v}}}{{dy}},\quad \varepsilon = p - \frac{{\alpha {{\eta }_{s}}}}{{\operatorname{Re} }}\left( {{{\mu }_{1}} + \frac{1}{3}} \right)\sigma ,\quad \Delta = \frac{{{{d}^{2}}}}{{d{{y}^{2}}}} - {{\alpha }^{2}}$

Поскольку в дальнейших расчетах в системе (1.13)–(1.18) использовались профили гидродинамических величин, зависящие от автомодельной переменной Дородницына–Хоуарта (1.10), то производные по y или, что то же самое по y′, заменялись производными по ξ в соответствии с соотношением

$\frac{d}{{dy}} \equiv \frac{d}{{dy{\text{'}}}} = \frac{{{{{\rho }}_{s}}}}{2}\frac{d}{{d{\xi }}}$

Необходимо отметить, что соотношение ε, содержащее параметр µ1 (объемную вязкость), исключается из системы. Для этого достаточно продифференцировать уравнение (1.14) по y, уравнение (1.15) умножить на iα и вычесть одно из другого полученные уравнения. Вместе со сделанным выше замечанием о независимости уравнений пограничного слоя в приближении Прандтля от объемной вязкости это позволяет констатировать, что объемная вязкость, как и прямо связанная с ней вторая вязкость в понимании Стокса [13], не влияют на линейную устойчивость течения в пограничном слое (сравни [11, 12]).

2. СПЕКТРЫ НЕВЯЗКИХ ВОЗМУЩЕНИЙ

Система (1.13)–(1.18) вместе с однородными граничными условиями определяет спектральную задачу, в которой собственными значениями являются комплексные фазовые скорости возмущений c = cr + ici.

Для расчета гидродинамических параметров невозмущенного течения система автомодельных уравнений (1.8), (1.9) приводилась к системе в нормальной форме, которая интегрировалась численно методом “стрельбы” с помощью процедуры Рунге–Кутты четвертого порядка на интервале [0, δc]. Точкой “прицеливания” служила середина интервала, где требовалось совпадение значений вычисляемых величин с точностью до 10–8.

В расчетах использовалось число Прандтля Pr = 0.75, условная верхняя граница пограничного слоя δс = 8, характерное время колебательной релаксации τ = 1, а зависимость сдвиговой вязкости от температуры описывалась формулой Сазерленда в виде

${\eta }(T) = \frac{{1.5{{T}^{{3/2}}}}}{{T + 0.5}}$

Примеры рассчитанных профилей скорости Us, плотности ρs и температуры Ts приведены на рис. 1.

Рис. 1.

Автомодельные профили скорости Us(ξ) (а), температуры Ts(ξ) (б) и плотности ρs(ξ) (б) невозмущенного пограничного слоя: М = 2 (1), 5 (2); 3 – зависимости Ts(ξ), 4 – зависимости ρs(ξ).

Для расчета спектра невязких возмущений в системе (1.13)–(1.18) опускались все диссипативные слагаемые, что эквивалентно переходу в ней к пределу при Re → ∞. При этом полученная “невязкая’’ система сводится к обыкновенному дифференциальному уравнению второго порядка для возмущения поперечной скорости v [10]. В автомодельных переменных оно записывается в виде

(2.1)
$\frac{{{{{\rho }}_{s}}}}{2}\frac{d}{{d{\xi }}}\left[ {\frac{{{{{\rho }}_{s}}}}{{2{\chi }}}\left( {W\frac{{d\text{v}}}{{d{\xi }}} - \text{v}\frac{{d{{U}_{s}}}}{{d{\xi }}}} \right)} \right] = \frac{{{{{\alpha }}^{2}}W\text{v}}}{{{{T}_{s}}}}$
(2.2)
$\text{v}({\text{0}}) = \text{v}({{{\delta }}_{{\text{c}}}}) = 0$

Здесь введены обозначения

${\chi } = {{T}_{s}} - {{m}^{2}}{{{\text{M}}}^{2}}{{W}^{2}},\quad W = {{U}_{s}} - c,\quad {{m}^{2}} = m_{r}^{2} + im_{i}^{2}$
$m_{r}^{2} = \frac{{{{Q}_{0}}(1 + {{{\gamma }}_{\text{v}}} + {\alpha \tau }{{c}_{i}}) + {{Q}_{1}}}}{{Q_{0}^{2} + Q_{1}^{2}}},\quad m_{i}^{2} = - \frac{{{{{\gamma }}_{\text{v}}}({\gamma } - 1){{Q}_{1}}}}{{{\gamma }(Q_{0}^{2} + Q_{1}^{2})}},\quad {{Q}_{0}} = 1 + \frac{{{{{\gamma }}_{\text{v}}}}}{{\gamma }} + {\alpha \tau }{{c}_{i}},\quad {{Q}_{1}} = {\alpha \tau }({{U}_{s}} - {{c}_{r}})$

Волновые числа невязких мод рассчитывались в соответствии со статьей [17]. Из условия обобщенной точки перегиба на автомодельных профилях Us(ξ), Ts(ξ)

$\frac{d}{{d{\xi }}}\left( {\frac{{\text{1}}}{{T_{s}^{2}}}\frac{{d{{U}_{s}}}}{{d{\xi }}}} \right) = 0$
находились координаты точек перегиба ξs в зависимости от числа Маха M и фазовые скорости в них cs = Uss). Затем для cs и числа M в уравнении (2.1) решалась спектральная задача (2.1), (2.2) относительно спектра вещественных волновых чисел αn (где n – номер моды возмущений) для дозвуковых возмущений при

$1 - \frac{1}{{\text{M}}} < {{c}_{s}} < 1,\quad \frac{{{{{\text{M}}}^{2}}}}{{{{T}_{s}}}}{{({{U}_{s}} - {{c}_{s}})}^{2}} > 1$

Задача решалась численно в среде пакета Matlab. Описание алгоритма расчета дано в [21, 22] (см. также [10, 18]).

Рассчитанные зависимости волновых чисел αn (n = 1, 2, 3, ...) для различных невязких мод от M приведены на рис. 2. Кривые 7 соответствуют идеальному газу. Предельное колебательное возбуждение практически не влияет на волновые числа невязких дозвуковых мод. Соответствующие значения даются точками на кривых рис. 2. Расчеты показали, что при M < 2.3 существует лишь первая невязкая мода, а старшие моды (n > 1) отсутствуют.

Рис. 2.

Зависимости αn(M) для первых шести невязких мод: мода I (1), II (2), III (3), IV (4), V (5), VI (6); 7 – идеальный газ, 8 – колебательно возбужденный газ при γv = 0.667.

Для каждого волнового числа αn из полученного набора при фиксированном числе Маха в интервале M = 2–10 решалась спектральная задача (2.1), (2.2) для комплекснозначных фазовых скоростей возмущений c = cr + ici. Графики фазовых скоростей cr и максимальных инкрементов нарастания ${\omega }_{{\,i}}^{{\max }}$ представлены на рис. 3. Возбуждение в данном случае в отличие от течения Куэтта [10, 17] понижает не только инкременты роста, но и фазовые скорости возмущений. Видно, что наиболее неустойчивой как в идеальном, так и в колебательно возбужденном газе является вторая акустическая мода. Можно также заметить, что в обоих случаях для чисел Маха M < 2.3 пограничный слой на пластине устойчив по отношению к двумерным невязким возмущениям.

Рис. 3.

Зависимости cr(M) (а) и ${\omega }_{i}^{{\max }}({\text{M}})$ (б) для первых четырех невязких мод: мода I (1), II (2), III (3), IV (4); 5 – идеальный газ, 6 – колебательно возбужденный газ при γv = 0.667.

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

(2.3)
${{\varepsilon }_{\varphi }} = \left| {\frac{{\varphi (0,{\text{M}}) - \varphi ({{\gamma }_{\text{v}}},{\text{M}})}}{{\varphi (0,{\text{M}})}}} \right| \times 100\% ,\quad {{\gamma }_{\text{v}}} = 0\, - \,0.667{\text{,}}\quad {\text{M}} = {\text{2}}\, - \,{\text{10}}$
Таблица 1.
М Идеальный газ Колебательно возбужденный газ
cr ${\omega }_{i}^{{\max }} \times {{10}^{3}}$ cr εc, % ${\omega }_{i}^{{\max }} \times {{10}^{3}}$ εω, %
Мода I
3 0.4348 0.2201 0.4222 2.8979 0.1932 12.2217
5 0.5362 1.0640 0.5207 2.8907 0.9342 12.1992
7 0.7925 0.7850 0.7696 2.8896 0.6892 12.2038
9 0.8475 0.5220 0.8230 2.8909 0.4583 12.2031
Мода II
3 0.4500 2.1301 0.4365 3.0000 1.8745 11.9994
5 0.6942 4.8801 0.6734 2.9963 4.2945 11.9998
7 0.4115 3.2500 0.3992 2.9891 2.8600 12.0000
9 0.4194 1.9601 0.4068 3.0043 1.7249 11.9994

Как следует из табл. 1, уменьшение фазовых скоростей возмущений, обусловленное колебательной неравновесностью, незначительно, а снижение инкрементов нарастания одинаково для обеих мод и достигает 12%.

3. ВЯЗКИЕ ВОЗМУЩЕНИЯ И КРИВЫЕ НЕЙТРАЛЬНОЙ УСТОЙЧИВОСТИ

Расчеты спектров вязких возмущений велись на основе полной системы уравнений (1.13)–(1.18), в которой число Рейнольдса ${{\operatorname{Re} }_{{\delta }}} = \sqrt {\operatorname{Re} } $ определено по местной толщине пограничного слоя. Собственными значениями служили фазовые скорости вязких возмущений c = cr + ici, а числа Рейнольдса Reδ, Маха M и волновое число α являлись параметрами. Спектральная задача (1.13)–(1.18) с однородными граничными условиями решалась численно в среде пакета Matlab. Использовался метод коллокаций [21, 22]. Подробное описание алгоритма расчета представлено в [2123] (см. также [10, 18]).

Вычисления велись при следующих значениях параметров: μ1 = 0–2, γv = 0–0.667, τ = 1, α = 0–0.3, Reδ = 102–4 × 103, M = 2–5, Pr = 0.75, γ = 1.4, δс = 8, Dα = 10–4, DReδ = 10.

Сравнение инкрементов (декрементов) ωi(α) = αci первой моды для совершенного и колебательно возбужденного газов представлено на рис. 4. Из графиков на рис. 4а видно, что для первой моды при M = 2.4 имеет место аномальное дестабилизирующее влияние вязкости как для совершенного газа, так и для колебательно возбужденного газа. Это влияние немонотонно по числу Рейнольдса. Можно заметить, что при Reδ = 400–1000 инкременты нарастания превышают значение в невязком пределе, затем для Reδ = 2 × 103 оказываются ниже невязкого предела, а при Reδ = 4 × 103 приблизительно совпадают с ним.

Рис. 4.

Зависимости ωi(α) для первой вязкой моды: (а) – M = 2.4, Reδ = 400 (1), 600 (2), 103 (3), 2 × 103 (4), 4 × 103 (5); (б) – M = 3.8, Reδ = 850 (6), 1.5 × 103 (7), 4 × 103 (8); A – кривые ωi(α) для идеального газа, B – совершенный газ, C – колебательно возбужденный газ при γv = 0.667, μ1 = 2.

При этом влияние колебательной неравновесности сводится к смещению зоны неустойчивости в сторону больших волновых чисел и небольшому снижению максимального инкремента. Следует отметить, что такая вязкая дестабилизация первой моды проявляется лишь в узком диапазоне чисел Маха в окрестности M ≈ 2.4. При всех других числах Маха вязкость оказывает на первую моду стабилизирующее воздействие. В частности, это подтверждается представленными на рис. 4б кривыми инкрементов для M = 3.8.

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

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

Расчет нейтральных кривых ωi(α, Reδ) = 0 для n-ой вязкой моды возмущений производился следующим образом. Предварительно устанавливались левая αmin и правая αmax границы расчетного интервала для волнового числа α. Для заданных значений параметров μ1, γv, M подсчитывалось число N попавших в интервал [αmin, αmax] волновых чисел αn невязких мод возмущений. Здесь n = 1, 2, ..., N – номера невязких мод. Далее расчетный интервал [αmin, αmax] разбивался на N + 1 интервал. Границами (m + 1)-го интервала α(m+ 1) = [αm, αm+ 1] служили значения волновых чисел, соответствующих двум соседним m-й и (m + 1)-й невязким модам возмущений. Здесь m = 0, 1, ..., N, а α0 = αmin и αN+ 1 = αmax.

Фиксируя значения параметров μ1, γv, M и номер m + 1 вязкой моды, вычислялись двумерные массивы инкрементов (декрементов) этой моды

(3.1)
${\omega }_{{i,kj}}^{{(m + 1)}} = {\omega }_{i}^{{(m + 1)}}({{{\alpha }}_{{m + 1,k}}},{{\operatorname{Re} }_{{{\delta },j}}}) = {{{\alpha }}_{{m + 1,k}}}{{c}_{{m + 1,k}}}({{{\alpha }}_{{m + 1,k}}},{{\operatorname{Re} }_{{{\delta },j}}})$
где одномерные массивы αm + 1, k и Rej рассчитывались по формулам

${{{\alpha }}_{{m + 1,k}}} = {{{\alpha }}_{{m,0}}} + k\Delta {\alpha },\quad {{{\alpha }}_{{m,0}}} = {{{\alpha }}_{m}},\quad m = 0,1,\; \ldots ,\;N,\quad k = 1,2,\; \ldots ,\;K$
${{\operatorname{Re} }_{{{\delta },j}}} = {{\operatorname{Re} }_{{{\delta },0}}} + j\Delta {{\operatorname{Re} }_{{\delta }}},\quad j = 1,2,\; \ldots ,\;J$

Здесь Δα, ΔReδ – шаги приращений по волновому числу и числу Рейнольдса соответственно.

Массив (3.1) определяет поверхность ${\omega }_{i}^{{(m + 1)}}({\alpha },{{\operatorname{Re} }_{{\delta }}})$ для (m + 1)-й вязкой моды. Координаты в множестве, определяющем геометрическое место точек данной линии уровня ${\omega }_{i}^{{(m + 1)}}({\alpha },{{\operatorname{Re} }_{{\delta }}}) = C$ на плоскости (α, Reδ), находились в соответствии с неравенством

${|\omega }_{i}^{{(m + 1)}}({{{\alpha }}_{{m + 1,k}}},{{\operatorname{Re} }_{{{\delta },j}}}) - C{\text{|}} \leqslant 1{{{\text{0}}}^{{ - 8}}}$
где C – некоторое заданное числовое значение. При C = 0 получается нейтральная кривая ${\omega }_{i}^{{(m + 1)}}({\alpha },{{\operatorname{Re} }_{{\delta }}}) = C$ для (m + 1)-й вязкой моды возмущений.

Кривые нейтральной устойчивости ωi(α, Reδ) = 0 первой и второй мод для различных чисел M, приведенные на рис. 5, дают общее представление об устойчивости течения в пограничном слое. Для Маха М = 2.2 существует только первая неустойчивая мода, известная как обобщение волны Толлмина–Шлихтинга на сжимаемый пограничный слой [17]. Из графиков на рис. 5а видно, что колебательное возбуждение аналогично случаю течения Куэтта [10, 18] сдвигает кривую в область больших волновых чисел и увеличивает критическое число Рейнольдса по сравнению с совершенным газом. На графике в качестве асимптоты нанесено значение волнового числа первой невязкой моды, к которому верхняя ветвь кривой устойчивости стремится снизу. Это означает, что вязкость при данном числе Маха оказывает стабилизирующее воздействие на первую моду.

Рис. 5.

Кривые нейтральной устойчивости ωi(α, Reδ) = 0: M = 2.2 (а), 4.5 (б), 4.8 (в); I – область неустойчивости моды I, II – область неустойчивости моды II; 1 – совершенный газ, 2 – колебательно возбужденный газ при γv = 0.667, μ1 = 2; ${{K}_{1}}$, $K_{1}^{'}$ – критические точки моды I, ${{K}_{2}}$, $K_{2}^{'}$ – критические точки моды II; α1 – асимптота при Reδ → ∞ для моды I, α2 – асимптота при Reδ → ∞ для моды II

При более высоких числах М = 4.5 и 4.8 общий характер влияния колебательного возбуждения сохраняется (см. рис. 5б, 5в). Здесь наиболее неустойчивой, как и в невязком случае, является вторая акустическая мода, для которой критические числа Рейнольдса более чем в два раза меньше, чем для первой моды. Из графиков на рис. 5в следует, что для М = 4.8 области неустойчивости первой и второй мод сливаются, как это отмечалось в [17]. Для обеих мод верхние ветви кривых подходят снизу к соответствующим асимптотическим невязким волновым числам, то есть вязкость в данном случае стабилизирует обе моды.

В табл. 2 приведены числовые значения критических чисел Рейнольдса Reδ, cr и волновых чисел αcr, соответствующие графикам рис. 5. Относительные отклонения εRe вычислялись по формуле (2.3). Как видно из табл. 2, колебательное возбуждение повышает значения Reδ, cr более чем на 13% для наиболее неустойчивой второй моды.

Таблица 2.
M Совершенный газ Колебательно возбужденный газ
αcr Reδ, cr αcr Reδ, cr εRe, %
Мода I
2.2 0.0580 281 0.0635 317 12.8114
4.5 0.0650 462 0.0708 521 12.7706
4.8 0.0701 440 0.0763 496 12.7273
Мода II
4.5 0.2001 221 0.2060 250 13.1222
4.8 0.1710 168 0.1772 191 13.6905

ЗАКЛЮЧЕНИЕ

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

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

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

При числе Маха M ⪢ 2 имеет место локальное дестабилизирующее влияние сдвиговой вязкости, проявляющееся как немонотонное по числу Рейнольдса Reδ поведение инкрементов нарастания вязкой первой моды.

Рассчитанные кривые нейтральной устойчивости для первой и второй наиболее неустойчивых вязких мод при числах Маха М = 2.2, 4.5 и 4.8 показывают, что для обеих мод критические числа Рейнольдса Reδ, cr при максимальном возбуждении превышают на 12–13% соответствующие значения для совершенного газа. При этом наибольшее возрастание Reδ, cr имеет место для наиболее неустойчивой второй моды.

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

  1. Kneser H.O. The interpretation of the anomalous sound-absorption in air and oxygen in terms of molecular collisions // J. Acoust. Soc. America. 1933. V. 5. P. 122–126.

  2. Леонтович М.А. Замечания к теории поглощения звука // Журн. экспер. и теор. физики. 1936. Т. 6. Вып. 6. С. 561–576.

  3. Landau L., Teller E. On the theory of sound dispersion // Collected Papers of L.D. Landau. Oxford: Pergamon Press, 1965. P. 147–153.

  4. Fujii K., Hornung H.G. Experimental investigation of high-enthalpy effects on attachment-line boundary layer // AIAA J. 2003. V. 41. № 7. P. 1282–1291. https://doi.org/10.2514/2.2096

  5. Beierholm A.K.-W., Leyva I.A., Laurence S.J., Jewell J.S., Hornung H.G. Transition Delay in a Hypervelocity Boundary Layer using Nonequilibrium CO2 Injection. GALCIT Technical Report FM 2008.001. Pasadena: California Institute of Technology. October 2008. 132 p.

  6. Bertolotti F.B. The influence of rotational and vibrational energy relaxation on boundary-layer stability // J. Fluid Mech. 1998. V. 372. P. 93–118. https://doi.org/10.1017/S0022112098002353

  7. Johnson H.B., Seipp T., Candler G.V. Numerical study of hypersonic reacting boundary layer transition on cones // Phys. Fluids. 1998. V. 10. P. 2676–2685. https://doi.org/10.1063/1.869781

  8. Wagnild R.M., Candler G.V., Leyva I.A., Jewell J.S., Hornung H.G. Carbon dioxide injection for hypervelocity boundary layer stability // AIAA Paper 2010–1244. January 2010. P. 1–16.

  9. Mack L.M. On the Inviscid Acoustic-Mode Instability of Supersonic Shear Layer. Part 1: Two-dimensional Waves // Theoretical and Computational Fluid Dynamics.1990. V. 2. № 2. P. 97–123. DOI 19910033977

  10. Grigoryev Yu.N., Ershov I.V. Stability and Suppression of Turbulence in Relaxing Molecular Gas Flows. Cham: Springer Intern. Publishing, 2017. 233 p. https://doi.org/10.1007/978-3-319-55360-3

  11. Молевич Н.Е. Асимптотический анализ устойчивости плоскопараллельного пограничного слоя сжимаемого релаксирующего газа // Изв. РАН. МЖГ. 1999. № 5. С. 82–88.

  12. Завершинский И.П., Кнестяпин В.Н. Устойчивость трехмерных возмущений малой амплитуды в неравновесном сжимаемом пограничном слое // Теплофизика высоких температур. 2007. Т. 45. № 2. С. 235–242.

  13. Ferziger J.H., Kaper H.G. Mathematical Theory of Transport Processes in Gases. Amsterdam–London: North Holland publishing company, 1972. 550 p. = Ферцигер Дж., Капер Г.К. Математическая теория процессов переноса в газах. М.: Мир, 1976. 555 с.

  14. Dunn D.W., Lin C.C. On the Stability of the Laminar Boundary Layer in a Compressible Fluid // J. Aeron. Sciences. V. 22. № 7. 1955. P. 455–477. https://doi.org/10.2514/8.3374

  15. Lin C.C. The Theory of Hydrodynamics Stability. Cambridge: Cambridge Univ. Press, 1955. 155 p. = Линь Цзя-цзяо. Теория гидродинамической устойчивости. М: Изд-во иностр. лит., 1958. 194 с.

  16. Гапонов С.А., Маслов А.А. Развитие возмущений в сжимаемых потоках. Новосибирск: Наука. Сиб. отд-ние, 1980. 144 с.

  17. Mack L.M. Boundary Layer Stability Theory. Preprint of JPL Technical Report, Document 900 – 277, Rev. A. Pasadena: California Institute of Technology. May 1969. 272 p.

  18. Григорьев Ю.Н., Ершов И.В. Линейная устойчивость сверхзвукового течения Куэтта молекулярного газа в условиях вязкой стратификации и возбуждения колебательной моды // Изв. РАН. МЖГ. 2017. № 1. С. 11–27. https://doi.org/10.7868/S0568528117010078

  19. Kaye G.W., Laby T.H. Tables of Physical and Chemical Constants. London. New York. Toronto: Longmans, Green & Co., 1958. 248 p. = Кэй Дж., Лэби Т. Таблицы физических и химических постоянных. М.: Физматгиз, 1962. 248 с.

  20. Лойцянский Л.Г. Механика жидкости и газа. М.-Л.: Гос. изд-во технико-теорет. лит., 1950. 676 с.

  21. Canuto C., Hussaini M.Y., Quarteroni A., Zang T.A. Spectral Methods in Fluid Dynamics. Berlin, Heidelberg: Springer-Verlag, 1988. 564 p. https://doi.org/10.1007/978-3-642-84108-8

  22. Trefethen L.N. Spectral methods in Matlab. Philadelphia: Society for Industrial and Applied Mathematics, 2000. 160 p.

  23. Malik M.R. Numerical Methods for Hypersonic Boundary Layer Stability // J. Comp. Phys. 1990. V. 86. P. 376–413. https://doi.org/10.1016/0021-9991(90)90106-B

Дополнительные материалы отсутствуют.