Прикладная математика и механика, 2019, T. 83, № 5-6, стр. 749-769

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

Влияние сильной термической неравновесности колебательных степеней свободы молекул и диссоциации молекул на ламинарно-турбулентный переход при гиперзвуковом обтекании тонких конусов систематически изучалось в серии экспериментальных и численных исследований [3–5].

Первые численные расчеты устойчивости пограничного слоя термически неравновесного газа на пластине в полной постановке были выполнены Bertolotti [6]. Неравновесность соответствовала стационарным условиям полета в атмосфере на высоте $H$ = 12 км при числе Маха M = 4.5, когда происходит возбуждение колебательных мод, и условиям эксперимента в аэродинамической трубе, где, наоборот, энергия колебательных мод, “замороженных” при температуре торможения в расширяющемся потоке, в процессе релаксации повышает статическую температуру газа.

Устойчивость сверхзвукового пограничного слоя на пластине рассчитывалась [7] в более широком диапазоне чисел Маха M = 1.6–4.8 и для предельного уровня колебательного возбуждения при отсутствии диссоциации.

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

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

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

Асимптотические оценки критического числа Рейнольдса в пограничном слое колебательно возбужденного газа на пластине для случая двумерных дозвуковых возмущений впервые были получены в работе [9]. В близкой постановке [10] рассматривался случай трехмерных возмущений. Неравновесность создавалась введением в уравнение Ландау–Теллера для колебательной энергии (температуры) стационарного источника накачки, что представляет интерес в свете возможных приложений. Однако из-за несовершенства и противоречивости использованной в этих работах модели полученные результаты вызывают серьезную критику.

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

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

В отличие от более полной постановки [6, 7], в [9] рассматривается пространственно однородное уравнение Ландау–Теллера, не включающее конвективный перенос возбуждения стационарным течением. При этом процесс диффузии колебательного возбуждения, существенный из-за градиента колебательной энергии поперек пограничного слоя не учитывается [6, 7, 9]. Введенный в это уравнение источник никак не компенсируется отводом энергии из поступательных степеней свободы, что можно было бы сделать либо теплоотводом с поверхности пластины, либо включением в уравнение энергии однородной по пространству эндотермической реакции [13]. Очевидно, что это противоречит стационарности невозмущенного течения. Можно привести и ряд других более мелких, но неочевидных недостатков.

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

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

1. Исходные уравнения спектральной задачи. На основе линейной теории гидродинамической устойчивости рассматривается развитие двумерных дозвуковых возмущений при сверхзвуковом обтекании колебательно возбужденным газом плоской полубесконечной пластины. Начало декартовой системы координат $\left( {x,~y} \right)$ совпадает с носиком пластины, координата $x$ ориентирована по направлению невозмущенного потока, координата $y$ направлена в поток по нормали к пластине. Течение описывается системой уравнений двухтемпературной газовой динамики [5, 12].

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

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

где ${\mathbf{q}}\left( y \right)$ – вектор амплитуд возмущений; $u,v,\rho ,\theta ,{{\theta }_{v}}$, p – амплитуды возмущений продольной и поперечной скоростей, плотности, статической и колебательной температур и давления соответственно; α – вещественное волновое число вдоль периодической переменной x; $c = {{c}_{r}} + i{{c}_{i}}$ – комплексная фазовая скорость возмущения; $i$ – мнимая единица.

Подстановка ${\mathbf{q}}{\text{'}}$ в систему уравнений для возмущений приводит к уравнениям для их амплитуд, которые имеют вид [7]

где

Индексом “s” здесь и далее обозначены значения соответствующих величин в стационарном невозмущенном течении. В качестве граничных условий принималось, что при $y = 0$ и на условной верхней границе пограничного слоя $y = {\delta }$ все возмущения обращаются в нуль.

В монографиях [2, 12] показано, что для двухатомных газов при умеренном уровне термической неравновесности поведение температурных зависимостей коэффициентов сдвиговой ${\eta }$ и объемной ${{{\eta }}_{b}}$ вязкостей практически совпадает, при этом их отношение меняется в диапазоне: = 0–2. Ввиду этого, а также для последующего сопоставления полученных результатов расчетов с данными [7], в уравнениях (1.1)–(1.6) в качестве температурной зависимости для коэффициента объемной вязкости аналогично [7] используется формула ${{{\eta }}_{b}}\left( T \right) = {{{\mu }}_{1}}{\eta }\left( T \right)$, где ${{{\mu }}_{1}} = {{{\eta }}_{{b\infty }}}{\text{/}}{{{\eta }}_{\infty }}$ = 0–2, а температурная зависимость коэффициента сдвиговой вязкости ${\eta }\left( T \right)$ описывается законом Сазерленда [17].

Коэффициенты теплопроводности $k = {{k}_{t}} + {{k}_{r}}$ и ${{k}_{v}}$, обусловленные поступательно-вращательным и колебательным движением молекул, выражены через коэффициент сдвиговой вязкости с помощью полуэмпирических соотношений Эйкена [7, 16]:

где ${{c}_{{Vt}}}$, ${{c}_{{Vr}}}$ и ${{c}_{{Vv}}}$ – удельные теплоемкости при постоянном объеме, связанные соответственно с поступательным, вращательным и колебательным движением молекул. Предполагается, что коэффициенты теплоемкости постоянны. Принимается также, что поступательные и вращательные степени свободы молекул находятся в состоянии квазиравновесия. При этом для их энергии имеет место равнораспределение по степеням свободы, и соответствующие теплоемкости определяются соотношениями: ${{c}_{{Vt\infty }}}$ = = 3R/2 и ${{c}_{{Vr\infty }}} = R$, где $R$ – газовая постоянная.

В уравнениях системы (1.1)–(1.6) коэффициент ${\gamma } = {{c}_{{p\infty }}}{\text{/}}{{c}_{{V\infty }}}$ – показатель адиабаты, ${{c}_{{V\infty }}} = {{c}_{{Vt\infty }}} + {{c}_{{Vr\infty }}}~$ и ${{c}_{{p\infty }}} = {{c}_{{Vt\infty }}} + {{c}_{{Vr\infty }}} + R~$ – соответственно удельные теплоемкости при постоянных объеме и давлении при температуре ${{T}_{\infty }}$. Параметр ${{{\gamma }}_{v}}$ = ${{c}_{{V{v}\infty }}}{\text{/}}\left( {{{c}_{{Vt\infty }}} + {{c}_{{Vr\infty }}}} \right)$ = = ${{c}_{{V{v}\infty }}}{\text{/}}{{c}_{{V\infty }}}$ характеризует степень неравновесности колебательной моды молекул, ${\tau }$ – характерное время релаксации колебательной степени свободы молекул. Критерии подобия ${\text{Re}} = {{{\rho }}_{\infty }}L{{U}_{\infty }}{\text{/}}{{{\eta }}_{\infty }}$, ${\text{M}} = {{U}_{\infty }}{\text{/}}\sqrt {{\gamma }R{{T}_{\infty }}} $ и ${\text{Pr}} = {{c}_{{p\infty }}}{{{\eta }}_{\infty }}{\text{/}}\left( {{{k}_{{t\infty }}} + {{k}_{{r\infty }}}} \right)$ = ${{c}_{{p\infty }}}{{{\eta }}_{\infty }}{\text{/}}{{k}_{\infty }}$ – соответственно числа Рейнольдса, Маха и Прандтля невозмущенного потока. С использованием формул Эйкена параметр ${{{\mu }}_{2}}$ для двухатомных газов выражается как μ₂ = = ${{k}_{{v\infty }}}{\text{/}}\left( {{{k}_{{t\infty }}} + {{k}_{{r\infty }}}} \right)$ = $20{{{\gamma }}_{v}}{\text{/}}33$.

Система линейных уравнений (1.1)–(1.6) имеет восьмой порядок и вместе с однородными граничными условиями определяет несамосопряженную спектральную задачу для развивающихся во времени возмущений, где собственными значениями являются комплексные фазовые скорости возмущений $c = {{c}_{r}} + i{{c}_{i}}$. Необходимо отметить, что соотношение ε, содержащее параметр (объемную вязкость), исключается из системы. Для этого достаточно продифференцировать уравнение (1.2) по y, уравнение (1.3) умножить на $i\alpha $ и вычесть одно из другого полученные уравнения. Уравнения пограничного слоя в приближении Прандтля, которые используются для расчета невозмущенного стационарного течения, также не зависят от объемной вязкости. Это позволяет констатировать, что объемная вязкость, как и прямо связанная с ней вторая вязкость в понимании Стокса [11], не влияют на линейную устойчивость течения в пограничном слое.

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

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

Здесь ${\xi }$ – автомодельная переменная Дородницына–Хоуарта [11]

а ${{{\delta }}_{c}}$ – условная верхняя граница пограничного слоя, причем

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

В расчетах здесь и далее использовалось число Прандтля ${\text{Pr}}$ = 0.75, за условную верхнюю границу пограничного слоя принималось значение ${{{\delta }}_{c}}$ = 8, зависимость сдвиговой вязкости от температуры описывалась формулой Сазерленда в виде [17]

Как известно, константа Сазерленда слабо зависит от температуры и в основном определяется родом газа. При этом относительная константа $\bar {C} = C{\text{/}}{{T}_{\infty }}$ меняется в пределах $\bar {C}$ = 0.3–1.0 [17]. В данном случае выбор $\bar {C}$ = 0.5 связан с необходимостью сравнения результатов с данными работ [7, 8].

Расчеты были выполнены для чисел Маха ${\text{M}}$ = 2.2, 4.5, 4.8. Примеры рассчитанных профилей скорости ${{U}_{s}}$, плотности ${{{\rho }}_{s}}$ и температуры ${{T}_{s}}$ приведены на рис. 1.

3. Асимптотика “невязких” решений. Асимптотические решения системы (1.1)–(1.6) для больших значений чисел Рейнольдса ищутся в виде ряда теории возмущений

В нулевом приближении, эквивалентном предельному переходу ${\text{Re}} \to \infty $ в уравнениях (1.1)–(1.6), приходим к системе уравнений для невязких возмущений, которая рассматривалась ранее в [15]. Было показано, что она может быть сведена к линейному уравнению второго порядка для возмущения давления. Таким образом, в этом приближении могут быть построены только два линейно независимых решения. Остальные шесть решений должны быть получены с учетом вязких слагаемых. При этом для вывода секулярного уравнения здесь надо использовать только затухающие при $y \to \infty $ решения [19].

Так как в коэффициентах уравнения для возмущения давления используются автомодельные профили скорости ${{U}_{s}}$ и температуры ${{T}_{s}}$, удобно перейти в нем к дифференцированию по автомодельной переменной (2.3) в соответствии с соотношением

Преобразованное таким образом уравнение в несамосопряженной форме записывается в виде

Для дальнейших вычислений в уравнении (3.1) проводится линейная замена независимой переменной ${\psi } = {\xi } - {{{\xi }}_{c}}$, которая, не меняя форму уравнения, переносит начало отсчета поперечной координаты в критический слой, где реальная фазовая скорость возмущения равна местной скорости невозмущенного потока ${{c}_{r}} = {{U}_{{sc}}}$. Для рассматриваемых ниже нейтральных возмущений примем ${{c}_{r}} \equiv c$. Индексом “c” будут обозначаться величины, вычисляемые в критическом слое при ${\psi } = 0$.

В уравнении (3.1) введены обозначения

Для нейтральных возмущений, как и для растущих, необходимо выполнение первого условия (теоремы) Рэлея [18]

Поэтому ${\psi } = 0$ представляет собой регулярную особую точку уравнения (3.1). Так как коэффициент $a\left( {\psi } \right)$ имеет в этой точке полюс первого порядка, а $b\left( {\psi } \right)$ является гладкой функцией, решение в окрестности особой точки строится методом Фробениуса [20] в виде обобщенного степенного ряда

Предварительно разложим коэффициенты уравнения (3.1) в степенные ряды по переменной ${\psi }$ в окрестности особой точки ${\psi } = 0$. Очевидно, что разложение функции $W$ имеет вид

При этом разложение коэффициента $a\left( {\psi } \right)$ строится исходя из соотношения

Получаем

Первые члены разложения коэффициента $b\left( {\psi } \right)$ записываются как

Подставляя представление (3.2) в уравнение (3.1) с учетом разложений (3.3), (3.4) и приравнивая нулю выражение при низшей степени ${{{\psi }}^{{r - 2}}}$, получаем следующее определяющее уравнение для показателя r [20]

Так как ${{c}_{0}} \ne 0$, то отсюда следует, что (3.5) имеет два корня ${{r}_{1}} = 3$ и ${{r}_{2}} = 0$. Поскольку ${{r}_{1}} \ne {{r}_{2}}$ и ${{r}_{1}} - {{r}_{2}} = 3$ – целому числу, то в соответствии с общей теорией [20] два линейно независимых решения уравнения (3.1) строятся в виде

В силу однородности уравнения (3.1) без ограничения общности можно положить $c_{0}^{{\left( 1 \right)}} = c_{0}^{{\left( 2 \right)}} = 1$. Подстановка разложения ${{p}_{1}}\left( {\psi } \right)$ в (3.1) позволяет получить рекуррентные выражения для коэффициентов $\{ c_{k}^{{\left( 1 \right)}}\} $. В результате имеем

При подстановке ${{p}_{2}}\left( {\psi } \right)$ в (3.1) сингулярная часть отделяется, а неизвестные коэффициенты $\{ c_{k}^{{\left( 2 \right)}}\} $ также определяются по рекуррентным соотношениям, что позволяет записать второе решение в виде

где

Фундаментальное решение строится численно в виде

Интегрирование ведется в двух интервалах – $\left[ {{{{\psi }}_{{\varepsilon }}},{\delta } - {{{\xi }}_{c}}} \right]$ для ${\psi } > 0$ и $\left[ {{{{\psi }}_{{ - {\varepsilon }}}}, - {{{\xi }}_{c}}} \right]$ для ${\psi } < 0$. В первом из интервалов краевая задача решается методом стрельбы. На верхней границе ${\psi } = - {{{\xi }}_{с}}$ ставится условие $p = 0$. В окрестности критического слоя в обоих случаях в качестве условий берутся асимптотические разложения фундаментального решения (3.8) при малых значениях $~{{{\psi }}_{{ \pm {\varepsilon }}}}$. Константа $C$ заранее неизвестна и определяется итеративно. Задается ее некоторое значение при фиксированных значениях $c$ и $\alpha $. Уравнение интегрируется к внешней границе $~{{{\psi }}_{{\delta }}} = {\delta } - {{{\xi }}_{c}}~$, и константа меняется до тех пор, пока не будет удовлетворено граничное условие $p = 0$. С найденной константой уравнение (3.1) интегрируется с асимптотическим условием при ${\psi } < 0$ до нижней границы ${\psi } = ~ - {{{\xi }}_{c}}~$.

Расчеты выполнялись для чисел Маха ${\text{M}}$ = 2.2, 4.5 и 4.8. Для каждого числа Маха ${\text{M}}$ находилась координата обобщенной точки перегиба ${{{\xi }}_{s}}$ из уравнения

и соответствующее значение фазовой скорости ${{c}_{s}} = {{U}_{s}}\left( {{{{\xi }}_{s}}} \right)$. Диапазон изменения фазовой скорости для первой моды выбирался аналогично [8] в соответствии с неравенством а для остальных мод из интервала

Шаг по $c$ был выбран $\Delta c = {{10}^{{ - 4}}}$. Волновые числа изменялись в диапазоне $\alpha \in \left[ {0,0.3} \right]$ с шагом $\Delta \alpha = {{10}^{{ - 4}}}$.

При интегрировании в интервале ${\psi } > 0$ в выражении ${{m}^{2}}$ пренебрегалось мнимой частью $m_{i}^{2}$ в силу ее малости, как это было сделано в асимптотике. При этом отпадает необходимость расщеплять решение уравнения (3.1) на реальную и мнимую части. При интегрировании в интервале ${\psi } < 0$ решение расщеплялось на реальную и мнимую части.

После расчета массива значений $p~\left( { - {{{\xi }}_{c}}} \right) = P~\left( {c,{\alpha }} \right)$ граничные значения “невязкой” части фундаментального решения, входящие в секулярное уравнение, рассчитывались по формулам для возмущений скоростей и температур, следующих из системы (1.1)–(1.6) при ${\text{Re}} \to \infty $:

4. Асимптотика “вязких” решений. Построение асимптотики “вязких” решений осуществляется на основе упрощения полной линеаризованной системы (1.1)–(1.6), как это было сделано в [15, 19]. Упрощенная система имеет вид

Следует отметить, что упрощенная система, как и исходная, имеет восьмой порядок. Если перейти здесь к совершенному газу: ${{{\gamma }}_{v}} = {{{\theta }}_{v}} = 0$, то система трансформируется в “вязкую” систему Дана–Линя шестого порядка [21].

Надо отметить, что у системы (4.1)–(4.4) есть два тривиальных решения u = = $v = {\theta } = {{{\theta }}_{v}} = 0$ и $u = 1$, $v = - iy$, ${\theta } = {{{\theta }}_{v}} = 0$. Поскольку здесь необходимы только убывающие (ограниченные) при $y \to \infty $ решения, то единственное удовлетворяющее этому условию нулевое тривиальное решение будет заменяться “невязким” решением (3.10):

Процесс построения “вязких” решений покажем на примере уравнения импульсов (4.2), которое преобразуем следующим образом. Введем замены независимой [22] и зависимой [21] переменных по формулам

В новых переменных входящие в (4.2) производные выражаются как

Подстановка этих выражений в (4.2) дает

где

Поскольку строится асимптотическая теория для ${\text{Re}} \gg 1$, можно пренебречь, как это было сделано в [21], слагаемым $P\left( Y \right)$. Заметив еще, что

приходим к уравнению Эйри [22]

Линейно независимые решения (4.6) могут быть представлены в нескольких эквивалентных формах, например, функциями Эйри первого и второго рода ${\text{A}}{{{\text{i}}}_{{1,2}}}\left( z \right)$. В нашем случае удобно воспользоваться аппаратом обобщенных функций Эйри ${{A}_{k}}\left( {z,p} \right)$ [23], для которых имеется ряд необходимых в дальнейшем анализе соотношений. В частности,

где $j = 0,~\; \pm 1,\;~ \pm 2, \ldots $ и $k = 1,~2$.

Асимптотическое представление для обобщенных функций Эйри имеет вид [23]

При этом ${{A}_{1}}\left( {z,j~} \right)\sim {{A}_{ - }}\left( {z,j~} \right)$, ${{A}_{2}}\left( {z,j~} \right)\sim {{A}_{ + }}\left( {z,j~} \right)$.

Первое из соотношений (4.7) позволяет записать два линейно независимых решения (4.6) как

Из выражений (4.5), (4.11) для производной возмущения продольной скорости следует

Чтобы найти решение для возмущения продольной скорости, необходимо проинтегрировать последнее. Из перехода к независимой переменной ${\zeta }$ имеем для дифференциалов

Интегрируя $u_{{1,2}}^{{\text{'}}}$, получаем два “вязких” решения

Здесь при интегрировании из-под знака интеграла были вынесены медленно меняющиеся функции, как это сделано в [21].

Только одно из полученных решений удовлетворяет необходимому асимптотическому условию убывания при $y \to \infty $. Как следует из выражений (4.10), (4.12), оно выражается через обобщенную функцию Эйри первого рода первого порядка. Таким образом, используемое в дальнейшем анализе решение записывается в виде

Чтобы построить решения уравнений (4.3) и (4.4) для возмущений температур, выполним в них следующие преобразования. Приближенно примем, что ${\text{33/}}\left( {20\gamma } \right) \approx 1$. Введем вспомогательные функции

Переход обратно к температурам, получаем

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

где

Как и выше, делаем в уравнении (4.15) замену зависимых и независимых переменных

В результате приходим к уравнениям

которые с точностью до постоянного множителя в коэффициенте совпадают с уравнением (4.6). Как и выше, необходимо выбрать только убывающие при $y \to \infty $ решения, которые выражаются через функции Эйри первого рода

Отсюда в соответствии с соотношениями (4.14) получаем решение уравнения (4.3) для возмущения статической температуры

Поскольку в связанные уравнения (4.1), (4.2) входит только ${\theta }\left( y \right)$, то найденное решение позволяет отщепить уравнение (4.4), понизив таким образом порядок определителя в секулярном уравнении.

Два линейно независимых решения для возмущения поперечной скорости строятся следующим образом. Полагая в (4.1) ${\theta } = 0$, $u \ne 0$ и интегрируя с учетом (4.13), получаем

где, как и в (4.12), для записи решения через обобщенную функцию Эйри ${{A}_{1}}\left( {{\zeta },2} \right)$ из-под интеграла были вынесены функции, медленно меняющиеся по сравнению с ${{A}_{1}}\left( {{\zeta },1} \right).$

Аналогично, при ${\theta } \ne 0$, $u = 0$ выполняя интегрирование в (4.1) с подстановкой (4.16), находим

В результате два линейно независимых решения усеченной “вязкой” системы (4.1)–(4.3) записываются в виде

5. Секулярное уравнение. Для перехода к спектральной задаче необходимо потребовать, чтобы линейная комбинация независимых решений (3.10) и (4.19) удовлетворяла однородным граничным условиям системы (1.1)–(1.6), тогда

Чтобы однородная система (5.1) имела нетривиальное решение $\left( {{{c}_{1}},~{{c}_{2}},{{c}_{3}}} \right)$, ее определитель должен быть равен нулю

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

Используя систему уравнений для невязких возмущений и формулы (3.9), получим выражение

Соотношение (5.4) позволяет разделить в уравнении (5.3) “невязкую” и “вязкую” части. В результате имеем

Для преобразования “вязкой” правой части уравнения (5.5) введем функцию Титьенса

и вспомогательную функцию значения которых табулированы [24, 25].

Функция Эйри первого рода связана с функцией Ханкеля первого рода порядка 1/3 соотношением [23]

которое позволяет выразить функции $F\left( z \right)~$ и $G\left( z \right)$ через используемые в статье обобщенные функции Эйри в виде

Выполнив в решениях (4.13) и (4.17) замену ${\zeta } = - z$, с учетом соотношений (5.6) получаем

Другое отношение “вязких” решений, входящее в уравнения (5.5), преобразуем следующим образом:

На основе численных расчетов кривых нейтральной устойчивости [7] имеет место оценка

В результате получаем

Вычислив производные при $y = 0$, множители перед функциями в (5.7) и (5.8) можно записать в виде

Используя соотношения (5.7)–(5.9), преобразуем уравнение (5.5), разрешив полученное выражение относительно ${\Pi }\left( 0 \right)$. В результате получаем секулярное уравнение в виде

которое используется в расчетах нейтральных кривых $N = \frac{{\left( {{\gamma } - 1} \right){{c}^{2}}{{{\text{M}}}^{2}}}}{{{{T}_{s}}\left( 0 \right)}},\quad 1 + {\lambda } = \frac{K}{c}.$

6. Асимптотики секулярного уравнения. Функции $F\left( z \right)$ и $G\left( z \right)$ табулированы [24, 25] до значений $z \leqslant 10$. В проводимых расчетах требуется знать решение секулярного уравнения для больших значений $z$. Для расчета функций $F\left( z \right)$ и $G\left( z \right)$ для больших значений аргумента использовались асимптотические представления [19, 21]:

где $q = {{z}^{{ - 3/2}}}{\text{/}}\sqrt 2 $.

Представляет интерес приближенная оценка минимального критического числа Рейнольдса ${\text{R}}{{{\text{e}}}_{{{\text{cr}}}}}$, подобная полученной ранее [15, 25] для совершенного газа. Необходимо отметить, что для оценки отдельных членов секулярного уравнения (5.10) при асимптотиках $\alpha \operatorname{Re} \gg 1$, $\alpha \ll 1$ в [15, 25] использовались профили параметров несущего потока, близкие к автомодельным профилям, получаемым из системы (2.1), (2.2). В частности, профиль скорости аппроксимировался сочленением прямой и дуги параболы, а температура получалась из интеграла Крокко аналогичной системы. Это позволяет воспользоваться с минимальными изменениями ранее полученными результатами [26], детально изложенными в [15, 25]: $1 + {\lambda } \cong 1$ и $N\Delta G \ll 1$. При этом уравнение (5.10) преобразуется к виду

где $\tilde {F}\left( z \right)$ – модифицированная функция Титьенса, которая также табулирована [24, 25]. Соответственно, левая “невязкая” часть уравнения (5.10) представляется в форме

При выводе соотношения (6.2) использовано представление “невязких” решений системы (1.1)–(1.6) рядами по степеням ${{\alpha }^{2}}$, в котором удержаны только первые члены разложения [15, 25], а функция ${{K}_{1}}\left( c \right)$ здесь имеет вид

При малых волновых числах ${\alpha }$

и мнимая часть ${{K}_{1}}\left( c \right)$ остается неизменной, а вещественная часть также пренебрежимо мала [15, 25].

Вершине кривой нейтральной устойчивости, дающей значение ${\text{R}}{{{\text{e}}}_{{{\text{cr}}}}}$, соответствует максимум ${{\tilde {F}}_{i}}\left( z \right)$ = 0.58, при значении аргумента $z$ = 3.24. При этом ${{\tilde {F}}_{r}}\left( z \right)$ = 1.46. Из соотношений (6.1) и (6.2) следует система для определения критических значений ${{c}_{{{\text{cr}}}}}$ и ${{{\alpha }}_{{{\text{сr}}}}}$:

Здесь

а индексом “c”, как и выше, обозначаются величины, вычисляемые в критическом слое при $y = {{y}_{c}}$.

После нахождения значений ${{c}_{{{\text{cr}}}}}$ и ${{\alpha }_{{{\text{сr}}}}}$ критическое число Рейнольдса, определенное по текущей длине $L$ пластины, вычисляется как

а критическое число Рейнольдса, определенное по местной толщине пограничного слоя ${\delta }$, по формуле

В табл. 1 представлены полученные из решения уравнений (6.3)–(6.5) числовые оценки критических чисел Рейнольдса ${\text{R}}{{{\text{e}}}_{{{\delta },{\text{cr}}}}}$, волновых чисел ${{{\alpha }}_{{{\text{cr}}}}}$ и фазовых скоростей ${{c}_{{{\text{cr}}}}}$ для совершенного (${{{\gamma }}_{v}} = 0$) и колебательно возбужденного (${{{\gamma }}_{v}} = 0.667$) газов. Видно, что при максимальном уровне колебательного возбуждения оценки значений критических чисел Рейнольдса ${\text{R}}{{{\text{e}}}_{{{\delta },{\text{cr}}}}}$ примерно на 7% выше, чем их значения в совершенном газе.

7. Численные расчеты нейтральных кривых. В расчетах спектров вязких возмущений число Рейнольдса ${\text{Re}}$, определенное по текущей длине $L$, было заменено на число ${\text{R}}{{{\text{e}}}_{{\delta }}} = \sqrt {{\text{Re}}} $, определенное по местной толщине пограничного слоя. Секулярное уравнение (5.10) имеет структуру, аналогичную структуре уравнений, рассматривавшихся ранее [15, 19]. Левая “невязкая” часть этого уравнения зависит от фазовой скорости $c$ и волнового числа ${\alpha }$, а правая “вязкая” часть выражается через табулированные функции переменной $z$ и также зависит от $c$. Поэтому координаты точки на кривых нейтральной устойчивости ${\text{R}}{{{\text{e}}}_{{\delta }}}\left( {{\alpha },~{{{\gamma }}_{{v}}},{\text{M}}} \right)$ на плоскости $\left( {{{{\operatorname{Re} }}_{\delta }},~\alpha } \right)$ рассчитывались в той же последовательности, что и в [15, 19]. Все вычисления проводились для времени колебательной релаксации $\tau $ = 1.

Для фиксированных значений числа Маха M и степени колебательной неравновесности ${{{\gamma }}_{v}}$ задавалось очередное значение фазовой скорости $c$ в интервале $c$ = [0, 1] с шагом ${\Delta }c = {{10}^{{ - 4}}}$. Интеграл в правой части секулярного уравнения (5.10) рассчитывался по формуле Симпсона [27]. Действительные и мнимые значения правой части уравнения (5.10) вычислялись для интервала $z$ = [0, 10] с помощью таблиц [24, 25], а для интервала $z$ = [10, 50] – с помощью асимптотической формулы [19, 21]. В расчетах принимался в обоих случаях табличный шаг ${\Delta }z$ = 0.04.

После расчета массивов правой части секулярного уравнения (5.10) вычислялись массивы его левой части для заданного значения $c$ в зависимости от волнового числа $c$ = [0, 0.5] с шагом $\Delta \alpha = {{10}^{{ - 4}}}$. Рассчитанные массивы правой и левой частей (5.10) сравнивались между собой до получения равенства с точностью до ${{10}^{{ - 8}}}$, если оно достигалось при фиксированном значении $c$. Далее расчет повторялся для очередного значения фазовой скорости. В результате были сформированы массивы значений волновых чисел ${{\alpha }_{k}}$, фазовых скоростей ${{c}_{k}}$ и переменной ${{z}_{k}}$, соответствующие точкам на нейтральной кривой. С использованием формулы

полученной путем замены независимой переменной в уравнении (4.9), вычислялись значения числа Рейнольдса ${\text{R}}{{{\text{e}}}_{{{\delta },k}}}$ на нейтральной кривой.

Нейтральные кривые для различных мод возмущений при фиксированном значении числа Маха M определялись по интервалам значений фазовой скорости ${{c}_{k}}$ на них. При этом в случае ${{c}_{k}} = \left( {1 - {{{\text{M}}}^{{ - 1}}},~{{c}_{{s1}}}} \right]~$ точки принадлежат кривой нейтральной устойчивости моды I, а в случае ${{c}_{k}} = \left( {{{c}_{{s1}}},{{c}_{{s2}}}} \right]~$ кривым нейтральной устойчивости моды II. Для моды $n$ (где $n$ = 1, 2, …) значения фазовых скоростей ${{c}_{{sn}}}$ при заданном числе Маха определялись из равенства ${{c}_{{sn}}} = {{U}_{s}}\left( {{{{\xi }}_{{sn}}}} \right)$, а координаты точек перегиба ${{{\xi }}_{{sn}}}$ – из уравнения (3.9).

На рис. 2 приведены результаты расчетов кривых нейтральной устойчивости ${\text{R}}{{{\text{e}}}_{{\delta }}}\left( {\alpha } \right)$ для совершенного (${{{\gamma }}_{{v}}} = 0$) и колебательно возбужденного (${{{\gamma }}_{{v}}}$ = 0.667) газов при числах Маха ${\text{M}}$ = 2.2 и 4.5. Сплошные кривые соответствуют решению исходной спектральной задачи (1.1)–(1.6), а штриховые – асимптотическому решению (5.10). Точки ${{K}_{1}}$, $K_{1}^{'}$ и ${{K}_{2}}$, $K_{2}^{'}$ на нейтральных кривых представляют собой критические точки для мод I и II, а прямые ${{\alpha }_{1}}$ и ${{\alpha }_{2}}$ – асимптоты при ${\text{R}}{{{\text{e}}}_{{\delta }}} \to \infty $ для первой и второй мод, полученные ранее Мэком в работе [8]. Области неустойчивости мод I и II на фигуре обозначены, соответственно, римскими цифрами I и II.

Расчеты показали, что при ${\text{M}}$ = 2.2 существует только первая неустойчивая мода, представляющая собой обобщение волны Толлмина–Шлихтинга на сжимаемый пограничный слой [8]. Из графиков видно, что асимптотические кривые удовлетворительно согласуются с результатами, полученными для исходной спектральной задачи (1.1)–(1.6). При этом переход к асимптотической теории приводит к расширению областей неустойчивости обеих мод и уменьшению критических значений их чисел Рейнольдса.

Из поведения кривых нейтральной устойчивости следует, что учет колебательной релаксации в уравнениях исходной спектральной задачи (1.1)–(1.6) и секулярном уравнении (5.10) (асимптотической теории) приводит к уменьшению размеров областей неустойчивости обеих мод и росту критических значений их чисел Рейнольдса.

В табл. 2 сравниваются значения критических чисел Рейнольдса ${\text{R}}{{{\text{e}}}_{{{\delta },{\text{cr}}}}}$ и волновые числа ${{\alpha }_{{{\text{cr}}}}}$, полученные из численного решения исходной спектральной задачи и даваемые асимптотической теорией. Видно, что для обеих мод критические значения чисел Рейнольдса, вычисленные по асимптотической теории, приблизительно на 14–15% меньше соответствующих значений, полученных при численном решении полной спектральной задачи. При этом учет колебательной релаксации приводит к возрастанию значений критических чисел Рейнольдса ${\text{R}}{{{\text{e}}}_{{{\delta },{\text{cr}}}}}$ обеих мод примерно на 12–13% относительно их значений в совершенном газе.

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

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

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

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