Журнал вычислительной математики и математической физики, 2019, T. 59, № 6, стр. 1024-1036

1. ПОСТАНОВКА ЗАДАЧИ

Рассматривается прямоугольный канал $D = [0,a] \times [0,b]$, заполненный жидкостью. Двумерные движения идеальной несжимаемой жидкости описываются уравнениями Эйлера в терминах завихренности течения $\omega (t,x,y)$ и функции тока $\psi (t,x,y)$:

Здесь $x$, $y$ – пространственные переменные, а $t$ – время. Компоненты скорости $(u,\text{v})$ выражаются через функцию тока $u = {{\psi }_{y}}$, $\text{v} = - {{\psi }_{x}}$. Вертикальные стороны $x = 0$ и $x = a$ прямоугольника $D$ соответственно вход ${{S}_{ + }}$ и выход ${{S}_{ - }}$ течения, а горизонтальные стороны $y = 0$ и $y = b$ – непроницаемые стенки.

Существует несколько вариантов математически корректной постановки краевых условий для задачи (1.1), (1.2) в канале, см. [13], [14]. В этой статье рассматриваются граничные условия В.И. Юдовича [2], [3].

На ${{S}_{ + }}$ задана завихренность и нормальная к границе компонента скорости втекающей жидкости

где ${{\gamma }_{ + }}$ и ${{\omega }_{ + }}$ – заданные функции, причем ${{\gamma }_{ + }}(y) > 0$ при $0 < y < b$.

На остальной части границы задается нормальная компонента скорости (нормаль внешняя):

При этом ${{\gamma }_{ - }}(y) > 0$ для $y \in (0,b)$.

Условие несжимаемости требует, чтобы заданные функции ${{\gamma }_{ - }}$ и ${{\gamma }_{ + }}$ удовлетворяли соотношению

Условия (1.3), (1.4) в терминах функции тока имеют вид

В начальный момент времени $t = 0$ $\omega $ задана во всем канале:

В статьях [2], [3] установлено, что двумерная задача протекания (1.1)–(1.7) глобально разрешима.

Левая часть уравнения (1.1) представляет собой материальную производную $\tfrac{{D\omega }}{{Dt}}$, откуда следует, что завихренность пассивно переносится жидкими частицами. Значение завихренности в частице жидкости сохраняется во времени и является интегралом системы.

Стационарные решения задачи определяются системой уравнений

с граничными условиями (1.6).

Первое уравнение в (1.8) удовлетворяется при существовании функциональной зависимости $\omega = F(\psi )$, т.е. поиск стационарных решений можно свести к решению уравнения

где $F(\psi )$ – некоторая функция, которую можно разыскивать численно или аналитически. Альтернативным подходом являются задание зависимости $F(\psi )$ и решение уравнения (1.9) относительно $\psi $.

Если известно некоторое стационарное решение $\hat {\psi }(x,y)$, $\hat {\omega }(x,y)$, то его устойчивость можно изучать, рассматривая динамику малых возмущений $\xi (x,y,t)$ и $\zeta (x,y,t)$. После подстановки $\psi = \hat {\psi } + \xi (x,y,t)$ и $\omega = \hat {\omega } + \zeta (x,y,t)$ в (1.1) и (1.2) получим систему уравнений относительно возмущений:

Граничные условия для (1.10) определяются видом предполагаемых возмущений. В этой статье рассматриваются возмущения, удовлетворяющие условиям:

Неустойчивость стационарных режимов можно показать, анализируя линеаризованную задачу (1.10), (1.11). В силу идеальности жидкости можно предположить существование нейтральной части спектра у линеаризованного оператора. Для исследования устойчивости решений в нейтральном случае необходимо решать полную задачу (1.10), (1.11) с ненулевыми начальными условиями.

Динамика частиц жидкости в $D$ описывается уравнениями

Функция тока $\psi $ является гамильтонианом системы (1.12), ее значение сохраняется для каждой жидкой частицы в области $D$, а траекториями частиц являются линии постоянного значения $\psi $. Это  свойство во многом определяет динамику идеальной жидкости с граничными условиями (1.3), (1.4).

В стационарном режиме протекания область течения разделяется на проточную зону ${{D}_{0}}$ и, возможно, застойные. Жидкая частица, входящая в область ${{D}_{0}}$, в точке $(0,{{y}_{ + }}) \in {{S}_{ + }}$, выходит из нее в точке ${{y}_{ - }} = R({{y}_{ + }})$, причем функция $R$ корректно определена на $y \in (0,b)$. При этом в силу того, что вихрь $\omega $ и функция тока $\psi $ являются интегралами уравнений движения частиц, должны выполняться равенства

Полученные соотношения приводят к равенству (см. [10])

которое выполняется во всей проточной зоне ${{D}_{0}}$ даже при наличии застойных зон. Таким образом, в проточной зоне однозначно определяется функция $F(\psi )$, которая задается условиями на границе (1.6). В каждой из рециркуляционных зон стационарных течений выполняется свое соотношение типа (1.14).

В этой статье расчеты проводились для конкретных функций, определяющих граничные условия (1.6). На входе ${{S}_{ + }}$ и выходе ${{S}_{ - }}$ из канала $D$ предполагается постоянная по времени и вертикальной координате нормальная компонента скорости, чему соответствуют функции

Здесь ${{Q}_{1}}$ – параметр, в расчетах принималось ${{Q}_{1}} = 0.05$.

Функция завихренности на входе в канал в условии (1.6) имела вид

В расчетах использовались значения параметров ${{y}_{c}} = \tfrac{b}{2}$, $\sigma = 50$. Параметр $A$ рассматривался в качестве бифуркационного для исследования возникновения автоколебаний.

При $A = 0$ существует тривиальное стационарное решение $\omega = 0$, $\psi = {{Q}_{1}}{\kern 1pt} y$ задачи (1.8), (1.6). Этому решению соответствует асимптотически устойчивое (см. [9]) проточное течение без застойных зон. При увеличении параметра $A$ при численном решении нестационарной задачи на больших временах устанавливаются стационарные режимы, а при дальнейшем росте $A$ – автоколебания (см. [7], [9]).

Сохранение значения завихренности в частицах жидкости позволяет реализовать эффективный бессеточный метод вихрей в ячейках. Метод детально описан в [8], [15], а его реализация для многопроцессорных вычислительных комплексов приведена в [16]. В основе метода расчета лежат следующие положения: поле завихренности $\omega (x,y)$ задано значениями в $N$ частицах; частицы пассивно переносят завихренность; динамика частиц описывается системой обыкновенных дифференциальных уравнений (1.12), которая решается псевдосимплектическим методом Рунге-Кутты; $\psi $ в каждый момент $t$ приближается отрезком ряда Фурье

где ${{g}_{{ij}}}(x,y) = sin\tfrac{{i\pi x}}{a}sin\tfrac{{j\pi y}}{b}$; коэффициенты ${{\psi }_{{ij}}}(t)$ находятся в результате применения метода Бубнова-Галеркина к уравнению (1.2); распределение $\omega (x,y)$ приближается двумерной кусочной аппроксимацией кубическими полиномами. Этот метод использовался для решения нестационарной задачи в данной работе.

Для изучения сценария возникновения автоколебаний в канале при росте параметра $A$ необходимо численно решать все перечисленные в этом разделе задачи. Далее описаны используемые методы и полученные с их помощью результаты. Приведены результаты расчетов для канала $D$ со сторонами $a = 3$, $b = 1$.

2. СТАЦИОНАРНЫЕ РЕЖИМЫ

Для поиска стационарных режимов численно решалась задача (1.9) с граничными условиями (1.6). Функциональная зависимость $\omega (\psi )$ в проточной зоне стационарного режима определяется граничными условиями функциями ${{\omega }_{ + }}(y),\;{{\psi }_{ + }}(y)$ и ${{\psi }_{ - }}(y)$. Для рассматриваемых в статье функций (1.15) и (1.16) эта зависимость имеет вид

Стационарные решения с функциональной зависимостью (2.17) являются решениями Уравнение (2.18) решалось численно методом конечных разностей. В области $D$ вводилась прямоугольная равномерная сетка размерности ${{n}_{x}} \times {{n}_{y}}$, где ${{n}_{x}}$, ${{n}_{y}}$ – количество узлов по координатам $x$ и $y$. Для аппроксимации производных применялись центральные разности второго порядка точности. В результате поиск стационарных режимов сводился к решению нелинейной системы алгебраических уравнений: Здесь ${{\psi }_{{i,j}}}$ – значения $\psi $ в узле с номером $i$, $j$, а ${{h}_{x}}$, ${{h}_{y}}$ – шаги сетки. Отметим, что неизвестными являются ${{\psi }_{{i,j}}}$ для $i = 2,\; \ldots ,\;{{n}_{x}} - 1$, $j = 2,\; \ldots ,\;{{n}_{y}} - 1$, а остальные задаются граничными условиями. Для решения системы (2.19) применялся итерационный метод Ньютона с аппроксимацией матрицы Якоби конечными разностями. В качестве начального приближения использовалась функция (1.15) при $A \ll 1$ и алгоритм продолжения по параметру $A$ для поиска решений при увеличении $A$. В результате решения находилось сеточное решение для функции тока $\{ {{\psi }_{{i,j}}}i = $ $ = \;1,\; \ldots ,\;{{n}_{x}},j = 1,\; \ldots ,\;{{n}_{y}}\} $, а значение завихренности находилось как $\{ {{\omega }_{{i,j}}} = {{F}_{1}}({{\psi }_{{i,j}}})i = 1,\; \ldots ,\;{{n}_{x}},$ $j = 1,\; \ldots ,\;{{n}_{y}}\} $.

Найденные сеточные решения приближались отрезками ряда Фурье по пространственным переменным $x$ и $y$:

Предполагается, что количество членов отрезка ряда (2.20) меньше, чем число узлов в (2.19). В расчетах принималось ${{k}_{x}} = {{n}_{x}} - 1$, ${{k}_{y}} = {{n}_{y}} - 1$.

Для вычисления коэффициентов ${{\hat {\psi }}_{{i,j}}}$, ${{\hat {\omega }}_{{i,j}}}$ отрезков ряда (2.20) использовалась процедура численного интегрирования произведения сеточного решения уравнения (2.19) на соответствующую тригонометрическую функцию. Такой способ аппроксимации стационарного режима позволяет вычислять значения $\hat {\psi }(x,y)$, $\hat {\omega }(x,y)$ и их частных производных в любой точке области $D$. Это необходимо для вычисления невязки уравнения (1.1) при контроле вычислений, анализа возмущенной задачи (1.10) и формирования начального условия (1.7), соответствующего стационарному режиму, для решения нестационарной задачи методом вихрей в ячейках.

Для контроля корректности найденного режима $\hat {\psi }(x,y)$, $\hat {\omega }(x,y)$ и проверки стационарности численного решения вычислялись следующие характеристики.

1. Невязка численного решения для (1.1)

2. Строилась функциональная зависимость ${{\omega }_{c}}({{\psi }_{c}})$, где ${{\omega }_{c}}$, ${{\psi }_{c}}$ – численные решения нестационарной задачи (1.1)–(1.7), полученные методом вихрей в ячейках с начальным условием ${{\omega }_{0}} = \hat {\omega }$. Вычислялось значение

Отметим, что малость величины ${{\delta }_{f}}$ при начальном условии ${{\omega }_{0}} = \hat {\omega }$ на больших временах ${{T}_{s}}$ говорит о возможной устойчивости стационарного режима. Для неустойчивого решения ${{\delta }_{f}}$ может быть малой только на небольших временах ${{T}_{s}}$ и должна расти с ростом ${{T}_{s}}$.

При расчетах для поиска стационарных режимов при $a = 3$, $b = 1$ и зависимости (2.17) использовалась сеточная аппроксимация (2.19) с числом узлов ${{n}_{x}} = 180$, ${{n}_{y}} = 60$. На фиг. 1 представлены режимы для четырех значений параметра $A$. Величины ${{\delta }_{\omega }}$ для этих режимов находятся в интервале от $3 \times {{10}^{{ - 6}}}$ (для $A = 0.1$) до $3 \times {{10}^{{ - 3}}}$ ($A = 1.5$). Характеристика ${{\delta }_{f}}$ для представленных режимов при ${{T}_{s}} = 5$ изменялась в интервале ${{10}^{{ - 3}}} < {{a}_{f}} < {{10}^{{ - 2}}}$.

Результаты вычислений показывают, что при малых $A$ стационарные режимы не имеют застойных зон, т.е. являются проточными (см. фиг. 1) для $A = 0.1,\;0.4$. При $A \approx 0.69$ возникает застойная зона, которая увеличивается в размерах при росте $A$ (см. фиг. 1). Для $A = 0.8,\;1.4$ у найденных режимов во всей области выполняется зависимость ${{F}_{1}}$ (см. (2.17)) а распределения завихренности на входе и выходе канала совпадают. Проточная зона включает струю, в центре которой завихренность является максимальной, а течение наиболее интенсивное. При росте $A$ струя приближается к границе $y = 0$ области $D$.

Можно предположить, что увеличение параметра $A$ ведет не только к изменению и усложнению топологии стационарных режимов, но и к изменению характера их устойчивости. Анализ устойчивости стационарных режимов проведен в следующем разделе.

3. АНАЛИЗ УСТОЙЧИВОСТИ СТАЦИОНАРНЫХ РЕЖИМОВ

Для анализа устойчивости найденных стационарных решений $\hat {\psi }(x,y)$, $\hat {\omega }(x,y)$ численно исследовалась задача (1.10). Для дискретизации системы уравнений (1.10) с граничными условиями (1.11) применялся метод прямых. Метод состоял в приближении решения $\xi (x,y,t)$, $\zeta (x,y,t)$ значениями ${{\xi }_{{i,j}}}(t)$ и ${{\zeta }_{{i,j}}}(t)$ в узлах пространственной прямоугольной сетки размерности ${{m}_{x}} \times {{m}_{y}}$ и аппроксимации производных $\xi (x,y,t)$ и $\zeta (x,y,t)$ по $x$ и $y$ конечными разностями второго порядка точности. В результате задача анализа устойчивости стационарных режимов сводится к исследованию устойчивости нулевого равновесия системы обыкновенных дифференциальных уравнений большого порядка.

В каждый момент времени $t$ значения ${{\xi }_{{i,j}}}$ в узлах сетки определяются из дискретного аналога второго уравнения системы (1.10) с граничными условиями (1.11):

Здесь ${{\xi }_{{i,j}}}$, ${{\zeta }_{{i,j}}}$ – значения $\xi $ и $\zeta $ в узле с номером $i$, $j$, а ${{\tau }_{x}}$, ${{\tau }_{y}}$ – шаги сетки. Неизвестными являются ${{\zeta }_{{i,j}}}$ для $i = 2,\; \ldots ,\;{{m}_{x}} - 1$, $j = 2,\; \ldots ,\;{{m}_{y}} - 1$, а значения ${{\zeta }_{{i,j}}}$ в граничных узлах равны нулю в силу (1.11).

Для внутренних узлов области $D$, т.е. при $i = 2,\; \ldots ,\;{{m}_{x}} - 1$, $j = 2,\; \ldots ,\;{{m}_{y}} - 1$, дискретизированные уравнения динамики возмущения завихренности ${{\zeta }_{{i,j}}}(t)$ имеют вид

Для аппроксимации пространственных производных для возмущения завихренности ${{\zeta }_{{{{m}_{x}},j}}}$ в граничных узлах, соответствующих выходу из канала, применялись разности против потока по переменной $x$, что обусловлено граничными условиями (1.11). Соответствующие уравнения имеют вид

В остальных граничных узлах ${{\xi }_{{i,j}}}(t) = {{\zeta }_{{i,j}}}(t) = 0$ в силу (1.11). Правые части уравнений (3.22) и (3.23) в качестве коэффициентов при переменных содержат значения частных производных $\hat {\psi }(x,y)$, $\hat {\omega }(x,y)$ по $x$ и $y$ в узлах сетки, которые вычисляются аналитически с использованием выражений (2.20).

В результате описанной дискретизации устойчивость стационарных режимов $\hat {\psi }(x,y)$, $\hat {\omega }(x,y)$, представленных на фиг. 1, исследовалась с помощью изучения свойств системы обыкновенных дифференциальных уравнений (3.22), (3.23) при изменении параметра $A$. Стационарному режиму задачи (1.1)–(1.7) соответствует нулевое равновесие системы (3.22), (3.23).

В первую очередь изучалась устойчивость нулевого равновесия системы (3.22), (3.23) по линейному приближению. Для этого вычислялась матрица Якоби системы и находились ее собственные значения ${{\lambda }_{i}}$, $i = 1,\; \ldots ,\;({{m}_{x}} - 1) \times ({{m}_{y}} - 2)$. Результаты вычислений представлены на фиг. 2. При описанных в этом разделе расчетах использовалась сетка с ${{m}_{x}} = 60$, ${{m}_{y}} = 20$ узлами. Для вычисления спектра матрицы Якоби применялся HQR-алгоритм.

При малых значениях $A$ спектр сконцентрирован в окрестности некоторой кривой на плоскости $({\text{Re}}{{\lambda }_{i}},{\text{Im}}{{\lambda }_{i}})$ и состоит из нейтральной части (чисто-мнимые собственные значения) и устойчивой (собственные значения с отрицательной действительной частью), см. фиг. 2, $A = 0.05$. При росте параметра $A$ структура спектра усложняется, но вплоть до $A \approx 0.8$ он состоит из устойчивой и нейтральной частей, см. фиг. 2 для $A = 0.15,\;0.4,\;0.7$. Наличие нейтрального спектра не позволяет сделать выводы об устойчивости или неустойчивости стационарного режима. При $A \approx 0.8$ стационарный режим становится неустойчивым, так как пара комплексно-сопряженных собственных значений переходят из левой в правую полуплоскость. Вычисленный спектр устойчивости стационарного режима для $A = 0.8$ изображен на фиг. 2. Это означает, что для возмущений стационарных режимов при $A = 0.8$ возникает неустойчивая колебательная мода, которая может приводить к возникновению автоколебаний.

При дальнейшем увеличении $A$ происходят аналогичные переходы, каждый из которых приводит к возникновению новых неустойчивых колебательных мод. В результате может сформироваться сложный нестационарный режим. Спектр устойчивости стационарного режима при $A \geqslant 0.8$ состоит уже из устойчивой, нейтральной и неустойчивой частей. Например, при $A = 1$, см. фиг. 2, неустойчивая часть найденного численно спектра состоит из восьми пар комплексно-сопряженных собственных значений с положительной действительной частью.

Для уточнения вопроса об устойчивости режимов при наличии нейтрального спектра численно решалась задача Коши для системы уравнений (3.22)–(3.23) с различными ненулевыми начальными значениями. Стремление этих решений к нулю является аргументом в пользу асимптотической устойчивости режима, а их рост говорит о неустойчивости. Вычислительные экспериментами проводились с различными начальными данными. Для всех вариантов расчетов при одном значении $A$ результаты качественно совпадали. Для оценки величины возмущения стационарного режима удобно использовать величину ${{\delta }_{u}} = ma{{x}_{{i,j}}}\left| {{{\zeta }_{{i,j}}}} \right|$.

Результаты вычислений для двух значений $A$ даны на фиг. 3, на которой изображена динамика величины ${{\delta }_{u}}$. При $A = 0.6$ часть собственных значений линеаризованной системы (3.22), (3.23) лежит на мнимой оси, а остальные имеют отрицательную действительную часть. При этом величина ${{\delta }_{u}}$ стремится к нулю, что говорит о затухании возмущений стационарного режима. Такой результат получается как при решении задачи Коши для сеточной аппроксимации (3.22), (3.23), так и при использовании метода вихрей в ячейках. Аналогичные результаты получены для всех рассматриваемых $A < 0.8$. То есть при $A < 0.8$ можно предположить асимптотическую устойчивость стационарных режимов.

Для неустойчивых стационарных режимов решение задачи (3.22)–(3.23) демонстрирует рост возмущений при росте $t$. При решении нестационарной задачи (1.1)–(1.7) методом вихрей в ячейках с начальными данными ${{\omega }_{0}} = \hat {\omega }(x,y)$ устанавливался нестационарный режим (см. следующий раздел). Эти расчеты подтверждают результаты численного анализа устойчивости по линейному приближению.

Отметим, что в рассматриваемой в этом разделе задаче устойчивости стационарных режимов возможен непрерывный спектр. Показать это или опровергнуть не представляется возможным, т.к. решения, устойчивость которых анализируется, получены численно. В силу дискретизации и использования численных методов поиска собственных значений матриц непрерывность спектра дифференциального оператора нарушается. Однако это не мешает получать качественно верные численные результаты об устойчивости или неустойчивости решений с непрерывным спектром (см. [17], [18]). Анализ спектра стационарных решений вида (20) с помощью дискретизации (3.22), (3.23) проводился для различных значений $10 < {{m}_{x}},{{m}_{y}} < 60$ при всех $A$. При этом бифуркационное значение появления неустойчивого спектра слегка “сдвигалось” по параметру $A$ при изменении количества узлов сетки, но сохранялось возникновение неустойчивых колебательных мод.

4. ВОЗНИКНОВЕНИЕ АВТОКОЛЕБАНИЙ

Переход при изменении параметра через мнимую ось пары комплексно-сопряжeнных собственных значений в спектре устойчивости равновесия конечномерной динамической системы сопровождается бифуркацией Андронова–Хопфа. При этом либо от равновесия ответвляется устойчивый предельный цикл (мягкая потеря устойчивости), либо неустойчивый предельный цикл в момент бифуркации схлопывается в равновесие (жeсткая потеря устойчивости). При наличии дополнительного нейтрального спектра при бифуркационном значении параметра разнообразие сценариев возрастает, но во многих случаях в окрестности равновесия возможно возникновение нестационарных режимов [19], [20]. Сказанное справедливо и при проведении численного анализа уравнений в частных производных, так как в этом случае в результате дискретизации возникают конечномерные динамические системы большой размерности.

Таким образом, описанный в разд. 3 сценарий возникновения неустойчивости стационарных режимов, при котором пара комплексно-сопряженных собственных значений переходит в положительную полуплоскость, может сопровождаться буфуркацией ответвления нестационарного режима. В рассматриваемой задаче ситуация сильно осложняется существованием большого числа чисто мнимых собственных значений в спектре устойчивости при бифуркационном значении, см. фиг. 2. Это и большая размерность систем (1.12) или (3.22), (3.23) затрудняет теоретическую интерпретацию происходящих бифуркаций и требует численного описания.

Для поиска и анализа нестационарных режимов численно решалась задача (1.1)–(1.7) с начальными условиями, близкими к стационарным режимам ${{\omega }_{0}} = \hat {\omega } + 0.001({{g}_{{1,1}}} + {{g}_{{2,3}}})$ на достаточно большом интервале $t \in [0,\;1000]$. Такой временной интервал соответствует 15-кратному прохождению жидкой частицей всего канала при $A = 0$. Для решения использовался метод вихрей в ячейках, который кратко описан в разд. 1, а подробно в [8], [15], [16]. В расчетах использовалось $N = 100\,000$ частиц, $60 \times 60$ членов отрезка ряда Фурье, шаг по времени $h = 1{\text{/}}200$ в методе Рунге–Кутты. Вычисления проводились для $0.05 \leqslant A \leqslant 5$.

В качестве критерия возникновения автоколебаний использовалось нарушение функциональной зависимости $\omega (\psi )$, нестационарность нормы $\left\| \text{v} \right\|$ скорости жидкости в точках вблизи выхода из канала $D$ с координатами $x = 295{\text{/}}300a$, $y \in [0.25b,0.3b]$ и интеграла завихренности $\Omega = \int_D \,\omega (x,y)dxdy$. Результаты представлены на фиг. 4 и фиг. 5. На фиг. 4 для каждого значения $A$ приведены три фрагмента. Ряд (а) изображает $\omega (\psi )$ во всей области $D$, ряд (б) в увеличенном масштабе для окрестности значения $\omega = A$ вблизи входа в канал ($x < 1{\text{/}}3$), а (в) – для $x > 2a{\text{/}}3$, т.е. вблизи выхода из $D$.

При значениях $A < 0.8$ начальные возмущения затухают со временем и при $t = 1000$ устанавливается стационарный режим. На фиг. 4 изображена численно построенная функция $\psi (\omega )$ для $A = 0.7$ и $t = 1000$. Четко прослеживается строгая функциональная зависимость $\psi (\omega )$, что говорит о стационарности установившегося режима. Это подтверждают и построенные графики зависимости $\left\| \text{v} \right\|$ от $t \in [900,\;1000]$ для шести точек, изображенные на фиг. 5, которые иллюстрируют постоянность во времени $\left\| \text{v} \right\|$. Данные вычисления подтверждают результаты анализа устойчивости стационарных режимов разд. 2 для $A < 0.8$, где предполагалась их асимптотическая устойчивость.

Увеличение $A$ до $0.8$ приводит к неустойчивости стационарного режима по линейному приближению. При $A = 0.8$ значение инкремента положительной действительной части собственных значений мало (см. фиг. 3), что обусловливает медленную и плохо обнаруживаемую в численном эксперименте нестационарную динамику. При росте $A$ нестационарные процессы нарастают вместе с увеличением количества положительных собственных значений в спектре устойчивости стационарного режима и их инкремента. Эти явления иллюстрируют результаты расчетов, представленные на фиг. 4 и фиг. 5 для $A = 1,\;1.5,\;2.5$. При $A = 1$ функциональная зависимость $\omega (\psi )$ незначительно нарушается вблизи выхода из канала $D$, см. нижний рисунок для $A = 1$, но сохраняется вблизи входа в канал (средний рисунок). При увеличении $A$ нарушение возрастет у выхода из канала и распространяется к зоне входа в канал, см. фиг. 4 для $A = 1.5$. Для $A = 2.5$ эти процессы проявляются еще сильнее. Во всех описанных случаях разрушение $\omega (\psi )$ наиболее сильно в окрестности максимального значения $\omega = A$ завихренности, что соответствует центру “струи” в проточной области течения.

Построенные численно графики зависимости нормы вектора скорости $\left\| \text{v} \right\|$ в шести точках вблизи выхода из $D$, см. фиг. 5, также иллюстрируют возникновение и развитие нестационарных режимов. Видно, что с ростом $A$ увеличивается амплитуда колебаний и усложняется характер динамики. Сложность динамики $\left\| \text{v} \right\|$, которая является квазипериодической или хаотической, объясняется тем, что для значений параметров, при которых проводились расчеты в спектре устойчивости стационарного режима, существует несколько пар комплексно-сопряженных собственных значений, см. фиг. 2. Комбинацией соответствующих им неустойчивых колебательных мод может объясняться непериодичность динамики.

При достаточно больших значениях $A$ автоколебания хорошо демонстрируют не только величина $\left\| \text{v} \right\|(t)$, но и поле завихренности, линии тока и значение интеграла завихренности $\Omega (t)$. На фиг. 6 для $A = 5$ в верхнем ряду изображено изменение поля $\omega $ во времени, а на нижнем графике сложная динамика $\Omega (t)$. Видно, что вблизи входа в канал распределение $\omega $ близко к стационарному и нарастает к выходу из $D$. Для $\Omega (t)$ проводился фурье-анализ сигнала с применением алгоритма БПФ, который демонстрирует непрерывный спектр, что говорит о хаотичности нестационарного режима при $A = 5$.

5. ЗАКЛЮЧЕНИЕ

В статье описаны расчеты для протекания жидкости сквозь плоский канал $D$ со сторонами $a = 3$, $b = 1$, но качественно аналогичные результаты получаются и для более длинных и коротких каналов. Во всех случаях наблюдалось не типичное для идеальной жидкости явление – возникновение автоколебаний при росте бифуркационного параметра с одинаковым сценарием. В качестве бифуркационного параметра в вычислительных экспериментах использовалась интенсивность завихренности жидкости на входе в канал, но аналогичные результаты могут быть получены и для других параметров, например, длины канала $D$.

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

Автоколебания возникают в результате мягкой потери устойчивости стационарным режимом при росте интенсивности завихренности на входе в канал (параметр $A$).

В докритической области $A$ (до потери устойчивости) можно предположить асимптотическую устойчивость стационарных режимов.

Неустойчивость возникает в результате перехода пары комплексно-сопряженных собственных значений из отрицательной в положительную полуплоскость. При этом в докритической и сверхкритической областях значений параметра $A$ в спектре устойчивости существуют дополнительные собственные значения на мнимой оси.

Потеря устойчивости стационарного режима, возможно, связана с возникновением застойной зоны в течении.

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

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

При больших надкритичностях реализуется развитый нестационарный режим с хаотической динамикой.

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