Известия РАН. Механика жидкости и газа, 2019, № 3, стр. 3-15
ЛИНЕЙНАЯ УСТОЙЧИВОСТЬ СВЕРХЗВУКОВОГО ПОГРАНИЧНОГО СЛОЯ РЕЛАКСИРУЮЩЕГО ГАЗА НА ПЛАСТИНЕ
Ю. Н. Григорьев a, b, *, И. В. Ершов a, c, **
a Институт вычислительных технологий СО РАН
Новосибирск, Россия
b Новосибирский государственный университет
Новосибирск, Россия
c Новосибирский государственный аграрный университет
Новосибирск, Россия
* E-mail: grigor@ict.nsc.ru
** E-mail: ivershov1969@gmail.com
Поступила в редакцию 16.07.2018
После доработки 18.10.2018
Принята к публикации 18.10.2018
Аннотация
На основе линейной теории устойчивости исследовано развитие невязких и вязких двумерных дозвуковых возмущений в сверхзвуковом пограничном слое колебательно возбужденного газа на плоской пластине. В качестве исходной модели использовалась система двухтемпературной газовой динамики, включающая релаксационное уравнение Ландау–Теллера. Невозмущенный поток описывался автомодельным погранслойным решением для совершенного газа. Показано, что в невязком приближении возбуждение снижает максимальные инкременты нарастания наиболее неустойчивой второй моды на 10–12% по сравнению с идеальным газом. Рассчитаны кривые нейтральной устойчивости для первой и второй наиболее неустойчивых мод для чисел Маха М = 2.2, 4.5, 4.8. Для обеих мод критические числа Рейнольдса при максимальном возбуждении превышают на 12–13% соответствующие значения для совершенного газа.
Эффект аномального поглощения ультразвука в молекулярных газах, связанный с конечным временем столкновительной релаксации внутренних степеней свободы молекул, был изучен и получил физическую интерпретацию в 30-х годах прошлого столетия [1–3]. В последние два десятилетия в ряде экспериментальных [4, 5] и теоретических исследований [6–8] было показано, что этот диссипативный эффект может существенно повлиять на устойчивость течений, если характерные частоты релаксации близки к частотам наиболее неустойчивых акустических мод возмущений.
В этом направлении наиболее систематические экспериментальные результаты, подтвержденные адекватным численным моделированием, были получены группой профессора Г. Хорнунга [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}}}$Продольная и поперечная составляющие скорости задаются равенствами u1 = u, u2 = v. Система (1.1)–(1.5) замыкается уравнением состояния идеального газа
Коэффициенты и параметры, входящие в уравнения (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)$Коэффициент γ = 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 выражается в виде
Система (1.1)–(1.6) линеаризуется на стационарном автомодельном погранслойном решении в отсутствие колебательного возбуждения [20]. Чтобы показать, что оно не зависит от объемной вязкости (параметра μ1), запишем стационарную систему уравнений динамики совершенного газа в погранслойном приближении. Для этого введем дополнительное “погранслойное” масштабирование поперечной координаты и скорости
Все остальные используемые ниже штрихованные переменные сохраняют обезразмеривание системы (1.1)–(1.6). В результате система (1.1)–(1.6) примет вид
Используя разложение по малому параметру
Автомодельные профили скорости и температуры рассчитывались численным интегрированием системы уравнений в автомодельных переменных
(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$Здесь автомодельная переменная Дородницына–Хоуарта [20]
(1.10)
${\xi } = \frac{{\text{1}}}{{\text{2}}}\int\limits_0^{y{\text{'}}} {\rho } {\text{'}}dy{\text{'}}$При выводе линеаризованных уравнений для возмущений мгновенные значения гидродинамических переменных представлялись в виде
(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{*}}$Учитывались температурные возмущения коэффициентов переноса
Здесь индексом s обозначены значения гидродинамических переменных, относящиеся к сдвиговому течению, а “*” – возмущения этих переменных.
Подстановка (1.11) в исходную систему (1.1)–(1.6) и линеаризация ее с точностью до членов первого порядка малости по возмущениям гидродинамических переменных приводит к следующей системе уравнений для возмущений:
Принималось, что при y = 0 и y = δc все возмущения обращаются в нуль.
Рассматривалась устойчивость периодических по продольной координате x возмущений в форме бегущих плоских волн
(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)$Подстановка (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.13)–(1.18) использовались профили гидродинамических величин, зависящие от автомодельной переменной Дородницына–Хоуарта (1.10), то производные по y или, что то же самое по y′, заменялись производными по ξ в соответствии с соотношением
Необходимо отметить, что соотношение ε, содержащее параметр µ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, а зависимость сдвиговой вязкости от температуры описывалась формулой Сазерленда в виде
Примеры рассчитанных профилей скорости Us, плотности ρs и температуры Ts приведены на рис. 1.
Для расчета спектра невязких возмущений в системе (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}}}}$Здесь введены обозначения
Волновые числа невязких мод рассчитывались в соответствии со статьей [17]. Из условия обобщенной точки перегиба на автомодельных профилях Us(ξ), Ts(ξ)
Задача решалась численно в среде пакета Matlab. Описание алгоритма расчета дано в [21, 22] (см. также [10, 18]).
Рассчитанные зависимости волновых чисел αn (n = 1, 2, 3, ...) для различных невязких мод от M приведены на рис. 2. Кривые 7 соответствуют идеальному газу. Предельное колебательное возбуждение практически не влияет на волновые числа невязких дозвуковых мод. Соответствующие значения даются точками на кривых рис. 2. Расчеты показали, что при M < 2.3 существует лишь первая невязкая мода, а старшие моды (n > 1) отсутствуют.
Для каждого волнового числа αn из полученного набора при фиксированном числе Маха в интервале M = 2–10 решалась спектральная задача (2.1), (2.2) для комплекснозначных фазовых скоростей возмущений c = cr + ici. Графики фазовых скоростей cr и максимальных инкрементов нарастания ${\omega }_{{\,i}}^{{\max }}$ представлены на рис. 3. Возбуждение в данном случае в отличие от течения Куэтта [10, 17] понижает не только инкременты роста, но и фазовые скорости возмущений. Видно, что наиболее неустойчивой как в идеальном, так и в колебательно возбужденном газе является вторая акустическая мода. Можно также заметить, что в обоих случаях для чисел Маха M < 2.3 пограничный слой на пластине устойчив по отношению к двумерным невязким возмущениям.
Числовые значения фазовых скоростей и максимальных инкрементов нарастания первых четырех мод невязких возмущений для идеального и колебательно-возбужденного газов, а также их относительные отклонения приведены в табл. 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]. Подробное описание алгоритма расчета представлено в [21–23] (см. также [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 приблизительно совпадают с ним.
При этом влияние колебательной неравновесности сводится к смещению зоны неустойчивости в сторону больших волновых чисел и небольшому снижению максимального инкремента. Следует отметить, что такая вязкая дестабилизация первой моды проявляется лишь в узком диапазоне чисел Маха в окрестности 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}}})$Здесь Δα, ΔReδ – шаги приращений по волновому числу и числу Рейнольдса соответственно.
Массив (3.1) определяет поверхность ${\omega }_{i}^{{(m + 1)}}({\alpha },{{\operatorname{Re} }_{{\delta }}})$ для (m + 1)-й вязкой моды. Координаты в множестве, определяющем геометрическое место точек данной линии уровня ${\omega }_{i}^{{(m + 1)}}({\alpha },{{\operatorname{Re} }_{{\delta }}}) = C$ на плоскости (α, Reδ), находились в соответствии с неравенством
Кривые нейтральной устойчивости ωi(α, Reδ) = 0 первой и второй мод для различных чисел M, приведенные на рис. 5, дают общее представление об устойчивости течения в пограничном слое. Для Маха М = 2.2 существует только первая неустойчивая мода, известная как обобщение волны Толлмина–Шлихтинга на сжимаемый пограничный слой [17]. Из графиков на рис. 5а видно, что колебательное возбуждение аналогично случаю течения Куэтта [10, 18] сдвигает кривую в область больших волновых чисел и увеличивает критическое число Рейнольдса по сравнению с совершенным газом. На графике в качестве асимптоты нанесено значение волнового числа первой невязкой моды, к которому верхняя ветвь кривой устойчивости стремится снизу. Это означает, что вязкость при данном числе Маха оказывает стабилизирующее воздействие на первую моду.
При более высоких числах М = 4.5 и 4.8 общий характер влияния колебательного возбуждения сохраняется (см. рис. 5б, 5в). Здесь наиболее неустойчивой, как и в невязком случае, является вторая акустическая мода, для которой критические числа Рейнольдса более чем в два раза меньше, чем для первой моды. Из графиков на рис. 5в следует, что для М = 4.8 области неустойчивости первой и второй мод сливаются, как это отмечалось в [17]. Для обеих мод верхние ветви кривых подходят снизу к соответствующим асимптотическим невязким волновым числам, то есть вязкость в данном случае стабилизирует обе моды.
В табл. 2 приведены числовые значения критических чисел Рейнольдса Reδ, cr и волновых чисел αcr, соответствующие графикам рис. 5. Относительные отклонения εRe вычислялись по формуле (2.3). Как видно из табл. 2, колебательное возбуждение повышает значения Reδ, cr более чем на 13% для наиболее неустойчивой второй моды.
ЗАКЛЮЧЕНИЕ
Исследование развития невязких и вязких двумерных дозвуковых возмущений в сверхзвуковом пограничном слое колебательно возбужденного газа на пластине позволяет сделать следующие выводы.
В рамках линейной теории устойчивости слабая термическая неравновесность внутренних степеней свободы молекул, описываемая объемной вязкостью, не влияет на устойчивость течения.
Возбуждение колебательных степеней свободы практически не влияет на волновые числа и фазовые скорости невязких мод, в то время как максимальные инкременты нарастания наиболее неустойчивой невязкой второй моды снижаются на 10–12% по сравнению с идеальным газом.
При числе Маха M ⪢ 2 имеет место локальное дестабилизирующее влияние сдвиговой вязкости, проявляющееся как немонотонное по числу Рейнольдса Reδ поведение инкрементов нарастания вязкой первой моды.
Рассчитанные кривые нейтральной устойчивости для первой и второй наиболее неустойчивых вязких мод при числах Маха М = 2.2, 4.5 и 4.8 показывают, что для обеих мод критические числа Рейнольдса Reδ, cr при максимальном возбуждении превышают на 12–13% соответствующие значения для совершенного газа. При этом наибольшее возрастание Reδ, cr имеет место для наиболее неустойчивой второй моды.
Список литературы
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.
Леонтович М.А. Замечания к теории поглощения звука // Журн. экспер. и теор. физики. 1936. Т. 6. Вып. 6. С. 561–576.
Landau L., Teller E. On the theory of sound dispersion // Collected Papers of L.D. Landau. Oxford: Pergamon Press, 1965. P. 147–153.
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
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.
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
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
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.
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
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
Молевич Н.Е. Асимптотический анализ устойчивости плоскопараллельного пограничного слоя сжимаемого релаксирующего газа // Изв. РАН. МЖГ. 1999. № 5. С. 82–88.
Завершинский И.П., Кнестяпин В.Н. Устойчивость трехмерных возмущений малой амплитуды в неравновесном сжимаемом пограничном слое // Теплофизика высоких температур. 2007. Т. 45. № 2. С. 235–242.
Ferziger J.H., Kaper H.G. Mathematical Theory of Transport Processes in Gases. Amsterdam–London: North Holland publishing company, 1972. 550 p. = Ферцигер Дж., Капер Г.К. Математическая теория процессов переноса в газах. М.: Мир, 1976. 555 с.
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
Lin C.C. The Theory of Hydrodynamics Stability. Cambridge: Cambridge Univ. Press, 1955. 155 p. = Линь Цзя-цзяо. Теория гидродинамической устойчивости. М: Изд-во иностр. лит., 1958. 194 с.
Гапонов С.А., Маслов А.А. Развитие возмущений в сжимаемых потоках. Новосибирск: Наука. Сиб. отд-ние, 1980. 144 с.
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.
Григорьев Ю.Н., Ершов И.В. Линейная устойчивость сверхзвукового течения Куэтта молекулярного газа в условиях вязкой стратификации и возбуждения колебательной моды // Изв. РАН. МЖГ. 2017. № 1. С. 11–27. https://doi.org/10.7868/S0568528117010078
Kaye G.W., Laby T.H. Tables of Physical and Chemical Constants. London. New York. Toronto: Longmans, Green & Co., 1958. 248 p. = Кэй Дж., Лэби Т. Таблицы физических и химических постоянных. М.: Физматгиз, 1962. 248 с.
Лойцянский Л.Г. Механика жидкости и газа. М.-Л.: Гос. изд-во технико-теорет. лит., 1950. 676 с.
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
Trefethen L.N. Spectral methods in Matlab. Philadelphia: Society for Industrial and Applied Mathematics, 2000. 160 p.
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
Дополнительные материалы отсутствуют.
Инструменты
Известия РАН. Механика жидкости и газа