Журнал вычислительной математики и математической физики, 2022, T. 62, № 4, стр. 677-693

1.1. Численный метод

В настоящей работе использована методика численного интегрирования полных уравнений Навье–Стокса, описанная в [3]. Уравнения решаются в криволинейной системе координат $(\xi ,\eta ,\zeta )$ в безразмерном дивергентном виде:

Декартовы координаты $x{\kern 1pt} * = xL$, $y{\kern 1pt} * = yL$, $z{\kern 1pt} * = zL$ отнесены к характерному масштабу длины $L$, время $t{\kern 1pt} * = tL{\text{/}}{{V}_{\infty }}$ – к характерному масштабу времени $L{\text{/}}{{V}_{\infty }}$, компоненты вектора скорости $u{\kern 1pt} * = u{{V}_{\infty }}$, $v{\kern 1pt} * = v{{V}_{\infty }}$, $w{\kern 1pt} * = w{{V}_{\infty }}$ – к модулю вектора скорости набегающего потока ${{V}_{\infty }}$, давление $p{\kern 1pt} * = p({{\rho }_{\infty }}V_{\infty }^{2})$ – к удвоенному скоростному напору в набегающем потоке, остальные газодинамические переменные – к их значениям в набегающем потоке. Звездочка в верхнем индексе обозначает размерную величину; если звездочка отсутствует, переменная полагается безразмерной. Символ “$\infty $” обозначает привязку к набегающему потоку. Для удобства сравнения с результатами других работ будут использоваться размерные переменные.

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

При определении конвективных потоковых величин на гранях ячейки сетки, например, ${\mathbf{E}}_{{i + 1/2,j,k}}^{{n + 1}}$, потоки расщепляются по направлениям. Для каждого направления определяется Матрица Якоби ($A = \partial {\mathbf{E}}{\text{/}}\partial {\mathbf{Q}}$ для направления $\xi $), которая диагонализируется в виде $A = B\Lambda {{B}^{{ - 1}}}$, где $B$ – матрица, составленная из правых собственных векторов, а $\Lambda $ – диагональная матрица собственных значений матрицы $A$. Конвективная составляющая потоковых величин ${\mathbf{E}}{\text{,}}\,\,{\mathbf{G}}{\text{,}}\,\,{\mathbf{F}}$ на грани ячейки аппроксимируется с использованием монотонной схемы типа Годунова (индексы i, j, k, n опущены для читаемости):

где нижними индексами L и R отмечены величины, рассчитанные слева и справа на рассматриваемой грани с использованием восстановленных с помощью процедуры реконструкции значений газодинамических переменных. Например, для грани i + 1/2 индекс L соответствует ячейке i, а индекс R – ячейке i + 1. Нижним индексом LR отмечены величины, вычисляемые с помощью приближенного метода Роу решения задачи Римана о распаде произвольного разрыва. Модификация собственных значений $\varphi (\lambda )$ обеспечивает физически корректное изменение энтропии на разрывах. В работе использована процедура реконструкции по методу WENO-3.

Функция $\Phi $ определяет степень немонотонности схемы и принимает значения от 0 до 1. Для исходной монотонной схемы $\Phi = 1$. Для гибридной схемы в расчетной области $\Phi = \max ({{\Phi }_{0}},\Psi )$, где ${{\Phi }_{0}}$ – некоторая постоянная величина, $\Psi $ – функция-индикатор Джеймсона [4]:

где ${\mathbf{V}}$ – локальный вектор скорости, ${{(\operatorname{rot} {\mathbf{V}})}^{2}}$– квадрат модуля вектора $\operatorname{rot} {\mathbf{V}}$.

Величина $\varepsilon $ выбирается малым положительным числом и позволяет устранить численную особенность деления на ноль в областях, где ${{(\operatorname{div} {\mathbf{V}})}^{2}}$ и ${{(\operatorname{rot} {\mathbf{V}})}^{2}}$ оба равны нулю. Величина $\Psi $ меняется от 0 при ${{(\operatorname{div} {\mathbf{V}})}^{2}} \ll {{(\operatorname{rot} {\mathbf{V}})}^{2}}$ до 1 в случае ${{(\operatorname{div} {\mathbf{V}})}^{2}} \gg {{(\operatorname{rot} {\mathbf{V}})}^{2}}$. В частности, в пограничном слое $\Psi \to 0$, а в окрестности ударных волн $\Psi \to 1$.

Перед выходной границей для подавления возмущений организуется буферная область (см. фиг. 1) с сеточным разрежением. В этой области вниз по потоку плавно восстанавливается монотонная схема ($\Phi $ стремится к единице) по формуле:

где $i$ – продольный сеточный индекс, $is$ – начальный индекс буферной зоны, $k$ – константа. В настоящей работе значение $k = 3$.

В криволинейной системе координат $(\xi ,\eta ,\zeta )$ расчетная сетка является равномерной по каждому из направлений с шагами ${{h}_{\xi }}$, ${{h}_{\eta }}$, ${{h}_{\varsigma }}$. Для аппроксимации диффузионной составляющей векторов ${\mathbf{E}}{\text{,}}\,\,{\mathbf{G}}{\text{,}}\,\,{\mathbf{F}}$ на грани элементарной ячейки применена разностная схема типа центральных разностей второго порядка точности:

где $q$ – любая из неконсервативных (“примитивных”) зависимых переменных задачи $(u,v,w,p,T)$, декартовы компоненты скорости, давление, температура.

После аппроксимаций уравнений Навье–Стокса и граничных условий интегрирование исходных уравнений в частных производных сводится к решению системы нелинейных алгебраических уравнений ${\mathbf{R}}({\mathbf{U}}) = 0$, где ${\mathbf{R}}$ – оператор дискретизации, вычисляющий вектор невязки согласно аппроксимации уравнений, ${\mathbf{U}}$ – вектор искомых неконсервативных переменных $(u,{v},w,p,T)$ (компоненты скорости, давление, температура) во всех узлах расчетной сетки. Длина вектора ${\mathbf{U}}$ составляет ${{n}_{q}}N$, где N – полное число узлов расчетной сетки, включая граничные, ${{n}_{q}}$ – число неизвестных в каждом узле (${{n}_{q}} = 5\left( {u,{v},w,p,T} \right)$ в трехмерной постановке и ${{n}_{q}} = 4\left( {u,{v},p,T} \right)$ в двухмерной). Система сеточных уравнений ${\mathbf{R}}({\mathbf{U}}) = 0$ решается с помощью модифицированного метода Ньютона–Рафсона

где $k$, ${{k}_{0}}$ – номера итераций по нелинейности, ${{k}_{0}} \leqslant k$, ${{{\mathbf{J}}}^{{[{{k}_{0}}]}}} = {{\left( {\partial {\mathbf{R}}{\text{/}}\partial {\mathbf{U}}} \right)}^{{[{{k}_{0}}]}}}$ – матрица Якоби системы нелинейных уравнений (матрица Якоби отображения), ${\mathbf{R}}({{{\mathbf{U}}}^{{[k]}}})$ – вектор невязки, $\tau $ – параметр регуляризации. Здесь выражение ${{\left( {{{{\mathbf{J}}}^{{[{{k}_{0}}]}}}} \right)}^{{ - 1}}}{\mathbf{R}}\left( {{{{\mathbf{U}}}^{{[k]}}}} \right) \equiv {{{\mathbf{Y}}}^{{[k]}}}$ является решением линейной системы уравнений Параметр регуляризации метода Ньютона относительно начального приближения ${{\tau }^{{[k]}}}$ определяется по формуле По мере сходимости итерационного [5] процесса ${{\tau }^{{[k]}}} \to 1$, а скорость сходимости теоретически стремится к квадратичной.

Получение аналитического вида матрицы Якоби ${\mathbf{J}}$ для рассматриваемой численной схемы, включающей решение задачи Римана о распаде разрыва, представляется весьма трудоемким. В данной работе применяется универсальный метод формирования матрицы ${{{\mathbf{J}}}^{{[{{k}_{0}}]}}} = {{\left( {\partial {\mathbf{R}}{\text{/}}\partial {\mathbf{U}}} \right)}^{{[{{k}_{0}}]}}}$ на итерации по нелинейности ${{k}_{0}}$ с помощью процедуры конечных приращений вектора невязки $R$ по вектору искомых переменных ${\mathbf{U}}$. При этом m-й столбец матрицы ${{{\mathbf{J}}}^{{[{{k}_{0}}]}}}$ вычисляется в виде

где ${{{\mathbf{e}}}_{m}}$ – единичный вектор длины ${{n}_{q}}N$, полностью состоящий из 0, кроме единственной 1 на позиции m. Такая методика вычисления Якобиана применима к произвольной системе сеточных уравнений.

1.2. Параметры потока и условия расчетов

Рассматривается номинально двухмерное течение над заостренной плоской пластиной при числе Маха набегающего потока ${{{\text{M}}}_{\infty }} = 3$ и температуре набегающего потока ${{T}_{\infty }} = 103.6$ К. Расчет эволюции возмущений проводится в подобласти; процедура расчета аналогична процедуре, описанной в [6]. Число Рейнольдса составляет ${{R}_{{1,\infty }}} = 2.181 \times {{10}^{6}}$ м^–1. Число Прандтля принимается постоянным: ${\text{Pr}} = \mu {{c}_{p}}{\text{/}}\lambda = 0.71$. Уравнения Навье–Стокса замыкаются уравнением состояния $\gamma {\text{M}}_{\infty }^{2}p = \rho T$, где $\gamma = 1.4$ – показатель адиабаты. Динамический коэффициент молекулярной вязкости рассчитывается по формуле Сазерленда: $\mu = (1 + {{T}_{\mu }}){\text{/}}(T + {{T}_{\mu }}){{T}^{{3/2}}}$, где ${{T}_{\mu }} = T_{\mu }^{*}{\text{/}}T_{\infty }^{*} = $ $ = 110.4\;{\text{K/}}103.6\;{\text{K}} \approx 1.07$.

Численное интегрирование проводится в прямоугольной области, показанной ниже. На входной и верхней границах фиксируются безразмерные параметры набегающего потока: $(u,v,w,p,T) = (1,0,0,1{\text{/}}\gamma M_{\infty }^{2},1)$. Для стационарных расчетов стенка принимается теплоизолированной, на стенке ставятся условия прилипания. Выходной границе предшествует буферная зона с укрупненными ячейками по продольной координате $x$ и нормальной к стенке координате $y$ для демпфирования выходящих через границу возмущений. На выходной границе накладываются мягкие условия как линейная экстраполяция примитивных переменных из расчетной области. На боковых границах ${{{\text{z}}}_{{{\text{min}}}}}$ и ${{{\text{z}}}_{{{\text{max}}}}}$ ставятся условия симметрии.

Расчет проводится следующим образом. Во-первых, двухмерное стационарное невозмущенное течение над плоской пластиной рассчитывается до достижения невязкой величины 10^–8. Во-вторых, из полученного решения вырезается подобласть, в которой далее будет моделироваться развитие возмущений; на новых входных границах подобласти фиксируются газодинамические величины из расчета на первом шаге; стационарное поле устанавливается дополнительно до полной сходимости (величина невязки не превосходит 10^–8). В-третьих, полученное в подобласти стационарное поле дублируется в третьем поперечном направлении $z$; распределение температуры по поверхности фиксируется; в пограничный слой вводятся возмущения типа “вдув–отсос” по процедуре, описанной ниже. Нестационарный расчет проводится до установления квазистационарного режима течения.

1.3. Функция-индикатор Джеймсона

Поле функции-индикатора Джеймсона [4] рассчитано по двухмерному невозмущенному течению и представлено на фиг. 2. “Пятнистость” значения этой функции за слабой ударной волной связана с численными особенностями ее расчета по формуле (1). В этой области значения слагаемых ${{(\operatorname{div} {\mathbf{V}})}^{2}}$ и ${{(\operatorname{rot} {\mathbf{V}})}^{2}}$ являются малыми величинами и могут быть одного порядка).

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

На фиг. 3 приведены профили $u(y)$ и $\Psi (y)$. Следует заметить, что область $\Psi \ll 1$ занимает порядка двух-трех толщин невозмущенного пограничного слоя. Поэтому на начальной стадии турбулизации течения пограничный слой продолжает находиться в области, где для аппроксимации конвективных составляющих потоковых величин используется низкодиссипативная схема типа центральных разностей. При проведении нестационарных расчетов форма профилей $\Psi (y)$ зафиксирована во всех сечениях $x = {\text{const}}$ и не зависит от текущего возмущенного поля течения.

1.4. Определение значения ${{\Phi }_{{\text{0}}}}$

При описанном подходе не удается полностью перейти на схему типа центральных разностей (${{\Phi }_{0}} = 0$), так как при турбулизации течения проявляется численная неустойчивость и процесс расчета перестает сходиться. Эмпирически удается определить наименьшее возможное значение ${{\Phi }_{0}}$, при котором численный процесс интегрирования дифференциальных уравнений остается устойчивым. Это значение зависит от качества сетки и от возмущений. На сетке 20 миллионов узлов ${{\Phi }_{{0,\min }}} = 0.35$, а на сетке 80 миллионов узлов – ${{\Phi }_{{0,\min }}} = 0.21$.

На фиг. 4 приведены изоповерхности Q-критерия, полученные на грубой сетке в 5 млн узлов в моменты времени, когда решение перестает сходиться, где Q – второй скалярный инвариант тензора скоростей деформаций

Как видно, во всех случаях возмущения не успевают дойти до конца расчетной области, причем, чем менее диссипативной является схема (чем меньше величина ${{\Phi }_{0}}$), тем раньше проявляются вычислительные проблемы. На рассматриваемой грубой сетке ${{\Phi }_{{0,\min }}} = 0.3$, однако, в этом случае амплитуда генерируемых возмущений $\varepsilon = 0.003$. Ниже в настоящей работе, как и в [1] рассматривается случай $\varepsilon = 0.00573$.

1.5. Генератор возмущений

Возмущения вводятся в пограничный слой при помощи генератора типа “вдув–отсос” с поверхности пластины. Генератор моделируется при $x{\kern 1pt} * \in [x_{1}^{*},x_{2}^{*}] = [0.394,0.452]$ м в соответствии с работами [1], [2]. В этом диапазоне нормальная составляющая вектора скорости имеет вид:

где где ${{T}_{0}} = 2\pi {\text{/}}{{\omega }_{0}}$, $A(t)$ – амплитуда, $\varepsilon = 0.00573$. Остальные параметры потока в области генератора вычисляются как для случая стенки без генератора. Возмущение с частотой $\omega _{0}^{*}{\text{/}}2\pi = f_{0}^{*} = 6.36$ кГц и волновым числом $\beta _{0}^{*} = 211.52$ м^–1 будем называть фундаментальным.

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

1.6. Расчетная сетка

Характерный масштаб длины выбран как $L = {\text{0}}{\text{.7239}}$ м. Продольный размер буферной зоны, которая ограничена пунктирным прямоугольником на фиг. 5, составляет полторы длины волны фундаментального возмущения, или $1.5{{\lambda }_{х}}$ где ${{\lambda }_{х}} = х_{2}^{*} - х_{1}^{*}$. В настоящей работе используются расчетные сетки из работы [1].

На фиг. 5 показаны расчетная область и расчетная сетка (вид сбоку). Подобласть для проведения основных нестационарных расчетов ограничена сплошным прямоугольником и начинается на расстоянии $x_{0}^{*} = 0.258$ м от передней кромки пластины. Длина подобласти в 14.3 разa больше продольной длины волны фундаментального возмущения. Высота подобласти выбрана $y_{H}^{*} = 0.03$ м, что составляет не менее пяти местных толщин пограничного слоя на выходной границе. Размер подобласти в боковом направлении составляет одну длину волны ${{\lambda }_{z}}$ в поперечном направлении, где ${{\lambda }_{z}} = 2\pi {\text{/}}\beta _{0}^{*} \approx 0.0297$ м.

Расчетная сетка представлена на фиг. 5а, соответствующие сеточные сгущения – на фиг. 5б и 5в. Основные расчеты настоящей работы проведены на сетке 80 миллионов узлов (подробная сетка). Эта сетка соответствует сетке из работы [2] в плоскости $xOy$. Сеточные линии распределены равномерно в поперечном направлении. Количество точек вдоль оси $z$ составляет 201 на подробной сетке. Грубая сетка имеет вдвое меньше узлов по $x$ и по $z$, чем подробная сетка. В вертикальном направлении количество точек для обеих сеток одинаково; поперек пограничного слоя приходится не менее 100 точек. На подробной сетке фундаментальное возмущение разрешено в боковом направлении (по $z$) 201 точками на длину волны по $z$, а в продольном направлении (по $x$) – 320 точками. Стоит отметить, что для распространения монохроматической акустической волны в равномерном потоке требуется около 40 точек на длину волны, чтобы добиться близкого к естественному уровню вязкого затухания волны для используемого диссипативного численного метода. Поэтому на построенных сетках численная диссипация фундаментального возмущения незначительна. Можно принять, что возмущения с длиной волны вчетверо короче длины фундаментальной волны моделируются достаточно достоверно, что позволяет “ухватить” процесс ламинарно-турбулентного перехода даже на грубых расчетных сетках.

1.7. Анализируемые величины

Данные для обработки и сравнения собираются, начиная с момента безразмерного времени $t = 2.261$ после ввода возмущений в пограничный слой, когда установился квазипериодический режим течения. В следующем разделе анализируются свойства переходного течения.

Анализ проводится для спектрального состава возмущений, где сопоставляются амплитуды отдельных гармоник Фурье или максимумы этих амплитуд по нормали к поверхности в рассматриваемом сечении $x = {\text{const}}$. Фурье-анализ и обработка нестационарных результатов выполнены с помощью возможностей языка программирования Python (библиотека numpy). Результат работы процедур быстрого преобразования Фурье по времени и координате $z$ нормирован на ${{N}_{t}} \times {{N}_{z}}{\text{/}}4$, где ${{N}_{t}}$ и ${{N}_{z}}$ – количество точек анализируемого сигнала по времени и по координате $z$ соответственно. В настоящей работе исследованы амплитуды гармоник пульсаций продольной компоненты скорости, давления, температуры, а также максимумы этих величины по нормали к поверхности.

Вихревая структура полей течения визуализируется с помощью -критерия: $Q = 0.5({{\Omega }_{{ij}}}{{\Omega }_{{ij}}} - {{S}_{{ij}}}{{S}_{{ij}}})$, ${{S}_{{ij}}} = 0.5({{\partial }_{j}}{{u}_{i}} + {{\partial }_{i}}{{u}_{j}})$, ${{\Omega }_{{ij}}} = 0.5({{\partial }_{j}}{{u}_{i}} - {{\partial }_{i}}{{u}_{j}})$, ${{u}_{i}}$ – компоненты вектора скорости (для записи использованы тензорные обозначения; предполагается соглашение о суммировании по повторяющемуся в произведении индексу).

Также сопоставляются поля мгновенных значений продольной и поперечной компонент вектора завихренности и средние параметры течения (коэффициент трения, профили продольной компоненты скорости, осредненные по Фавру и по Рейнольдсу).

2.1. Структуры течения

Мгновенная структура возмущенного течения представлена на фиг. 6 в виде изоповерхностей -критерия. Сразу за источником возмущений расположена область линейного развития возмущений, где они формируют X-образные структуры, усиливающиеся вниз по потоку. Вблизи $х{\kern 1pt} * \approx 0.6$ мм появляются первые признаки нелинейного взаимодействия – возмущения начинают искажаться. При $х{\kern 1pt} * \in (0.75,0.85)$ м наблюдается интенсивный нелинейный распад возмущений. За этой областью формируется зона молодой турбулентности с ростом мелкомасштабных вихрей. Эта зона развивается вниз по потоку. Можно видеть, что для гибридной схемы (настоящая работа) размер разрешаемых мелкомасштабных вихрей меньше, чем в случае монотонной схемы (см. [1]). Это различие наблюдается в области нелинейного режима при $х{\kern 1pt} * > 0.85$ м.

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

В каждом поперечном сечении $x{\kern 1pt} * = {\text{const}}$ возмущение можно представить через сумму гармонических колебаний путем преобразования Фурье. Для рассматриваемого течения вдоль плоской пластины разумно делать двухмерное преобразование Фурье по времени и поперечной координате для каждой линии $x{\kern 1pt} * = {\text{const}}$, $y{\kern 1pt} * = {\text{const}}$. Результат такого двухмерного преобразования можно представить в виде амплитуды гармоники $(f{\kern 1pt} *,\beta {\kern 1pt} *) = (hf_{0}^{*},k\beta _{0}^{*})$. Таким образом, результат двухмерного преобразования Фурье можно представить в виде амплитуд двухмерных гармоник пульсации продольной скорости $\hat {u}_{{hk}}^{'}$. Пример представлен на фиг. 7 для начала области молодой турбулентности.

В описанной постановке Фурье-спектр является симметричным, поэтому интерес представляет только часть спектра при $h \geqslant 0$, $k \geqslant 0$. Следует отметить, что пики спектра расположены в шахматном порядке, что объясняется квадратичным (нелинейным) взаимодействием возмущений друг с другом. Например, для стационарного возмущения и других четных частот $h = 0,2,4,...$ наблюдаются только максимумы на четных волновых числах $k = 0,2,4,6,...$, а для нечетных частот $h = 1,3,5,...$ – на нечетных волновых числах $h = 1,3,5,...$ . Такая картина свойственна механизму наклонного распада, когда нелинейно (квадратично) взаимодействуют две гармоники с одинаковыми частотами, но противоположными по знаку волновыми числами. При этом частота удваивается, а волновое число обнуляется $[1,1] + [1, - 1] \to [2,0]$. Чем ближе гармоника к фундаментальной гармонике, тем выше ее амплитуда. Это связано как с тем, что нелинейный распад продвигается в область высоких частот постепенно, так и с тем, что имеется численная пространственно-временная диссипация используемого численного метода. Расчеты настоящей работы выполнены на двух разных сетках, одна из которых вдвое мельче в продольном и боковом направлении. Из фиг. 7а и 7б следует, что спектр возмущений, полученных с применением гибридной схемы, оказывается несколько шире по сравнению со случаем монотонной схемы.

Ниже полученные результаты сопоставляются с результатами работ [1] и [2].

2.2. Линейный режим

Рассмотрим линейный режим развития возмущений, который наблюдается примерно от $x{\kern 1pt} * = 0.4$ м до $x{\kern 1pt} * = 0.6$ м. На фиг. 8 амплитуды фундаментальной моды $[1,\,1]$ пульсаций продольной компоненты вектора скорости $u{\kern 1pt} '$ в сечении $x{\kern 1pt} * = 0.5$ м, полученные в настоящей работе, сопоставляются с результатами работы [1] и [2]. Для данного случая и для других линий $x{\kern 1pt} * = {\text{const}},$ $y{\kern 1pt} * = {\text{const}}$, наблюдается хорошее согласование (фиг. 8).

Среди всех возможных линий $y{\kern 1pt} * = {\text{const}}$ для данного сечения $y{\kern 1pt} * = {\text{const}}$ можно выделить линию $y_{0}^{*}$, на которой амплитуда рассматриваемой гармоники максимальна. Для пульсаций $u{\kern 1pt} '$ или $T{\kern 1pt} '$ эта линия будет находиться в критическом слое пограничного слоя, ${{y}_{0}}{\text{/}}\delta \approx 0.65$.

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

На фиг. 9 представлены структуры возмущений внутри пограничного слоя в сечении $y{\kern 1pt} * = {\text{const}}$. Все поля хорошо согласуются друг с другом. Гибридная схема дает чуть более наполненные контуры по сравнению с монотонной схемой и лучше согласуется с работой [2], но эта разница незначительна.

2.3. Нелинейный режим

Рассмотрим нелинейную стадию развития возмущений. Момент проявления нелинейного взаимодействия можно отметить по картинам Q-критерия, на которых заметно усиление “канатообразных” структур (см. фиг. 10 и 11). Для малых значений $Q$ (15 и 100) (фиг. 10) на сетке 20 и 80 миллионов узлов гибридная и монотонная схемы дают почти одинаковые структуры значений $Q$. Они хорошо согласуются результатами работы [2].

Однако в области молодой турбулентности, где появляются мелкомасштабные структуры и максимальная величина $Q$ растет (10 000 и 40 000) (фиг. 11), диссипативная схема не позволяет идеально воспроизвести результаты низкодиссипативной схемы [2]. Наиболее вероятными причинами этого расхождения являются применение в [2] спектрального метода в боковом направлении и использование высокочастотных гармоник, которые располагаются вблизи частоты (волнового числа) Найквиста для подробной расчетной сетки настоящей работы и поэтому плохо разрешаются на ней.

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

Сделанный выше вывод подтверждается количественно с помощью эволюции амплитуд отдельных гармоник. Рассмотрим, в частности, эволюцию максимальных по $y$ амплитуд. Механизм наклонного резонанса проявляется последовательно. Изначально растет наиболее неустойчивая фундаментальная наклонная волна $[1, \pm 1]$, что обусловлено чисто неустойчивостью пограничного слоя. При достижении некоторой критической амплитуды она начинает нелинейно взаимодействовать сама с собой, порождая кратные гармоники: $h = 0$ и $h = 2$, $k = 0$, $k = 2$, которые начинают расти благодаря нелинейному взаимодействию фундаментальных гармоник. Когда кратные гармоники достигают достаточных амплитуд, они начинают нелинейно взаимодействовать друг с другом и с фундаментальными гармониками, порождая все больше кратных гармоник. Такой процесс и его временная последовательность прослеживаются на фиг. 12 и фиг. 13, на которых изображена эволюция амплитуд гармоник вниз по потоку. Описанный механизм объясняет шахматную структуру спектра, приведенного на фиг. 7.

На фиг. 12а представлена эволюция фундаментальной моды продольной компоненты скорости, полученная на сетке 80 млн узлов с использованием гибридной схемы (линия 5). Решение оказывается ближе к решению из работы [2], чем в случае монотонной схемы [1] (линия 3). На сетке 20 млн узлов решение с использованием гибридной схемы (линия 4) также аналогично ближе к [2] и почти достигает до решения с использованием монотонной схемы и сетки 80 млн узлов (линия 3). На фиг. 12б–г различие в линиях не такое четкое.

На фиг. 13а продемонстрировано, что на сетке 20 млн узлов решение с использованием монотонной схемы (линия 2) сильно отличается от решения [2]. Однако решение, полученное с использованием гибридной схемы (линия 4), значительно лучше совпадает с решением на базе монотонной схемы на сетке 80 млн узлов (линия 3).

Таким образом, гибридная схема дает решения, которые лучше согласуются с решением [2] по сравнению со случаем монотонной схемы [1]. На грубой сетке решение гибридной схемы может достигать точности решения монотонной схемы на мелкой сетке. Разницы в решениях с использованием монотонной и гибридной схем проявляются сильнее в низках модах.

На фиг. 14 приведены мгновенные структуры поперечного вектора завихренности в различных сечениях $z{\kern 1pt} * = {\text{const}}$ и сравнение полученных результатов с результатами работы [1] и [2]. Видно, что вблизи $х{\kern 1pt} * = 0.865$ м появляются мелкомасштабные структуры, свойственные развитой турбулентности. Основные крупномасштабные структуры настоящей работы и работы [1] хорошо совпадают с [2]. Детальное сравнение показывает, что вихри, полученные в работе [1], менее интенсивны по сравнению с вихрями из работы [2], и содержат меньше мелкомасштабных структур, как обсуждалось выше. При этом гибридная схема настоящей работы разрешает вихревые структуры лучше (фиг. 14г, д), чем монотонная схема [1] (фиг. 14б, в). Визуально гибридная схема позволяет достигать точности подробной сетки с использованием грубой сетки [2].

На фиг. 15 показаны мгновенные сечения продольной компоненты вектора завихренности, которые сопоставлены с результатами работы [1] и [2]. Структура течения симметрична относительно плоскости z* = 0; поэтому приводится только половина области по z*. В сечении $х{\kern 1pt} * = 0.862$ м (фиг. 15а) расположен большой вихрь. С ростом x он растет и, начиная с $х{\kern 1pt} * = 0.872$ м (фиг. 15б), распадается с образованием более мелких вихрей. Результаты, полученные на гибридной схеме (настоящая работа), лучше согласуются с результатами работы [2] в случае монотонной схемы (работа [1]).

Рассмотрим зависимость осредненной величины коэффициента трения ${{с}_{f}}$ от продольной координаты x. Локальный коэффициент трения вычисляется по формуле

Как и в работе [1], в настоящей работе рассматривается усреднение по времени в течение пяти периодов фундаментальной гармоники $\Delta t = 10\pi {\text{/}}{{\omega }_{0}}$и по размаху z всей расчетной области:

Начиная с $х{\kern 1pt} * = 0.72$ м, величина ${{с}_{f}}$ резко увеличивается в окрестности точки $х{\kern 1pt} * = 0.86$ м, что соответствует началу ламинарно-турбулентного. В этом диапазоне результаты для диссипативной и гибридной схем совпадают с результатами работы [2] даже на грубой сетке. Ниже по потоку результаты на грубой сетке оказываются ниже по сравнению со случаем подробной сетки. В особенности это проявляется для монотонной схемы, а в случае гибридной схемы соответствующая кривая мало отклоняется от случая подробной расчетной сетки $(x{\kern 1pt} * \geqslant 0.9\;{\text{м}})$.

На фиг. 17 показаны средние профили газодинамических переменных в сечении $х{\kern 1pt} * = 0.996$ м. Чем сильнее эффект нелинейных взаимодействий в рассматриваемом сечении, тем более наполнены средние профили.

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