Применение методов бифуркационного анализа для решения задач гидромеханики Текст научной статьи по специальности «Математика»

УДК 532.517,571.938

А. Н. Нуриев, А. Г. Егоров

ПРИМЕНЕНИЕ МЕТОДОВ БИФУРКАЦИОННОГО АНАЛИЗА

ДЛЯ РЕШЕНИЯ ЗАДАЧ ГИДРОМЕХАНИКИ

Ключевые слова: система уравнений Навье-Стокса, бифуркационный анализ, метод продолжения по параметру, метод

гомотопии.

В работе представлены численные методы бифуркационного анализа для исследования стационарной системы уравнений Навье-Стокса. Исследование включает локализацию решений системы, построение ветвей и идентификацию локальных бифуркаций системы. Рассмотрено применение методов для решения задачи о вторичном циркуляционном течении вокруг цилиндра, осциллирующего с большой частотой и малой амплитудой в вязкой жидкости.

Keywords: Navier-Stokes system of equations, bifurcation analysis, continuation method, homotopy method.

In the paper the bifurcation analysis methods for the study of stationary Navier-Stokes system of equations are presented. The study includes methods for location of different solutions, continuation of solution branches and identification of local bifurcations of the system. Application of these methods for solving the problem of the steady streaming around a circular cylinder oscillating with high frequency and small amplitude in a viscous fluid is considered.

Введение

Методы бифуркационного анализа (БА) представляют эффективный аппарат для исследования нелинейных систем, возникающих в различных предметных областях [1-4]. В настоящей работе рассматривается их применение в задачах гидромеханики: для локализации, продолжения и анализа решений дискретизированной стационарной системы уравнений Навье-Стокса.

Система уравнений Навье-Стокса является нелинейной однопараметрической системой с единственным безразмерным параметром - числом Рейнольдса (Яе). Влияние нелинейности определяется величиной этого параметра. При малых значениях Яе линейные члены доминируют, и система имеет единственную ветвь решения. С увеличением Яе влияние нелинейности возрастает, что приводит к появлению новых ветвей решения. Для многих задач уже при умеренных числах Рейнольдса существуют десятки различных по структуре решений. Основным методом исследования системы уравнений Навье-Стокса является сеточный. Дискретизация уравнений Навье-Стокса на подробных сетках приводит к дискретным системам большой размерности. Применение для поиска решений таких систем классических методов в большинстве случаев приводит к ограниченным результатам. Они не позволяют изучать все разнообразие решений. Помимо большой размерности, лимитирующими факторами здесь являются отсутствие начальных приближений для решений и наличие точек бифуркации у системы. Использование специализированных методов БА позволяет существенно расширить возможности численного исследования стационарных уравнений Навье-Стокса.

БА применяется для исследования задач гидромеханики достаточно редко, поэтому остановимся более подробно на его компонентах,

методах и целях. В данной работе БА реализуется в рамках классического подхода (описанного например в [4]) для анализа однопараметрических нелинейных систем. Его целью является построение возможных ветвей решения и их стратификации относительно бифуркаций фазового портрета системы. Анализ проводится в расширенном фазово-параметрическом пространстве, которое является прямым произведением фазового и параметрического пространств системы. Основными компонентами анализа являются задачи локализации решений, продолжения по параметру, исследования точек бифуркаций.

Задача локализации служит для поиска вещественных решений системы при фиксированном значении Яе. Прямое выделение решений невозможно, т.к. неизвестны ни их число, ни хорошие начальные приближения к ним. Как правило, однако, одно из возможных решений, непрерывно продолжающееся от малых чисел Рейнольдса, строится без особых затруднений. Это обстоятельство позволяет использовать для решения задачи локализации метод БРМ гомотопии [2]. Его реализация применительно к задачам вычислительной гидродинамики представлена в разделе 3 работы.

Задача продолжения по параметру призвана обеспечить движение по выделенным локализацией ветвям решения в фазово-параметрическом пространстве системы. Для ее решения рассматривается расширенная система уравнений, в которой помимо дискретных переменных, неизвестным считается также и параметр Яе системы. Метод решения основан на подходе «предиктор-корректор». В качестве предиктора в работе используется метод касательной, в качестве корректора метод Мура-Пенроуза [4]. Такая связка позволяет эффективно аппроксимировать ветвь даже в окрестностях точек бифуркаций. Подробно метод продолжения описан во втором разделе данной работы.

Исследование точек бифуркаций позволяет провести заключительную стратификацию ветвей и построить бифуркационную диаграмму. Основными локальными бифуркациями системы являются бифуркация складки и бифуркация Андронова-Хопфа [4]. Точки бифуркации первого типа диагностируются на этапе построения ветвей, так как соответствуют частным экстремумам по параметру системы. Более точная идентификация этих точек, а также локализация точек бифуркации Андронова-Хопфа проводится при решении спектральной задачи. Постановка этой задачи и схема её решения изложены в четвертом разделе работы.

Описанный в работе аппарат БА был использован для исследования вторичных стационарных течений «steady streaming» вокруг цилиндра, осциллирующего с большой частотой и малой амплитудой в вязкой жидкости [5-8]. Применение методов БА позволило обнаружить в области умеренных чисел Рейнольдса несколько решений и построить соответствующие бифуркационные диаграммы. Постановка задачи приведена в разделе 1, а основные результаты представлены в разделе 5.

1. Постановка и дискретизация задачи

Задача о вторичном стационарном течении вокруг осциллирующего цилиндра описывает вторичные стационарные течения, возникающие при периодических колебаниях цилиндра в вязкой несжимаемой жидкости в условиях малой амплитуды и высокой частоты колебаний [8].

Уравнения, описывающие стационарное вторичное течение в терминах вихрь у - функция

тока о получены в результате асимптотического анализа [6], предполагающего малую амплитуду и большую частоту колебаний цилиндра. Число Рейнольдса Re, построенное по линейной скорости колебаний и толщине стоксовского нестационарного пограничного слоя при этом полагается величиной порядка единицы. Вне тонкого погранслоя на цилиндре стационарное вторичное течение в главном члене удовлетворяет стационарной системе уравнений Навье-Стокса. В полярной системе координат (г, в) она имеет вид

Re\дв\удго - дгудво\/г = А® А у = -о

(1)

(2)

Граничные условия на цилиндре определяются [6] посредством сращивания с аналитическим решением в пограничном слое. На бесконечности задаются условия затухания: ш\ = 0, д гш\ = 0

1г=1 ' 1г=ю

дгщ| , = 3^т2Д, дгю\ = 0

' ' 1г=1 ' 1г=ю

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

к = (Н2,кД). Обозначая через (щк\) сеточную функцию тока и завихренность, получим

Re

е 2 z (дв) кУк(д z) к(с) - (д z) кУк (дв) к®к е 2 z

е 2Z (д z )Н (®) + (дв)Н®н

е

(3а)

Здесь,

z J к

аппроксимации

\=-(\,в)щк ■ (3б) (дг)к, (дв)к - центрально разностные соответствующих частных пPоизвоДных, (Аг Д =(дг )2 + (дД - дискретный

аналог оператора Лапласа.

На границе цилиндра г = 0 ив условной бесконечности г = записываются дискретные

аналоги граничных условий (2). Для сеточной функции тока они получаются из первой пары условий (2)

=0 = 0 (д. )Гщ|г=ги= 0. (4)

Здесь верхний индекс «и» указывает на использование направленной разности второго порядка точности для аппроксимации нормальной производной. Вторая пара условий в (2) преобразуется в граничные условия для завихренности:

Ю4=0 =[{К «1п2Д-Щк1=к

Iz=H^ J i 2 Z ) Hz

(д z )к

|(м),

- 2®h|z

= 0,

(5)

первое из которых по существу представляет условие Тома [9].

2. Задача продолжения

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

При малых числах Рейнольдса основная ветвь легко находятся с помощью метода Ньютона с нулевого начального приближения. Для продвижения по ветви в область больших чисел Рейнольдса решается задача продолжения [4]. Опишем общий алгоритм решения этой задачи. Далее он используется не только для продвижения по основной, но и любой другой локализованной ветви решения.

Представим исходную систему уравнений рассматриваемой задачи в следующей операторной форме:

Р(у) = 0 (6)

Здесь Р - нелинейный оператор задачи, у = (х, а) - расширенный вектор переменных, состоящий из вектора неизвестных дискретных значений функции тока и завихренности

zz

х = (ук,ак) и параметра задачи а = Яе. Задача продолжения состоит в построении последовательности точек у1, у2,—, Ут,

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

Построение к-ой точки по (к -1) -ой проводится на основе схемы «предиктор-корректор». На этапе «предиктор» начальное приближение для к-ой точки вычисляется с помощью экстраполяции по методу касательных:

У0 = Ук-1 +Аак ■ аУк-1.

Здесь - текущее значение шага,

выбираемое равным ак-1 - ак-2, а

dyk-1 = (аХк-1, Зак-1) - вектор касательной к кривой равновесия в точке Ук-1. Он находится из системы уравнений

Л (Ук-1 )аУк-1 = С1, а-1 =!, (7)

где Л (Ук-1) - якобиан расширенной системы в к-1 точке. Для коррекции начального приближения используется модификация итерационного метода Ньютона - метод Мура-Пенроуза [4]. На каждой итерации этого метода решается линейная система уравнений вида

л(уГ1)УГ = Р(уГ-1) - Л(УГ1)УГ-1.

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

1 п-1 п-1 ^ ^

вектору аук в точке ук возмущенной кривой

(уГ-УГ-1, аУП-1 ) = о. (9)

Здесь скобки означают скалярное произведение; нормализованный касательный вектор определяется решением системы уравнений

J(уГ'Ж-1 = о, ааГ-1 = 1. (10)

При расходимости итераций метода Мура-Пенроуза

шаг Аад- уменьшается вдвое, и процедура

коррекции проводится заново. Необходимость в этом возникает лишь в окрестности точек бифуркации складки.

При реализации описанного метода необходимо решать однотипные системы линейных уравнений (7), (8)-(9), (10) . Все они представимы в виде

J(У*)У = 0, (у+, у) = я, (11)

где величины у , у+, /, я являются заданными. Для решения системы линейных уравнений (11) используется метод обобщенных минимальных невязок (ОЫЯБ8) [10]. Этот подход часто применяется в комбинации с методом Ньютона для решения задач большой размерности (см. например [11]). В нашем случае выбор метода обусловлен нессиметричностью и отсутствием диагонального преобладания матрицы системы (11). Эффективность вМЯБ8 во многом зависит от

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

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

3. Задача локализации

Осуществить непосредственный переход от основной ветви к другим ветвям решения в фазово-параметрическом пространстве возможно, если они имеют точки пересечения (см. например [4]). В общем случае для выявления других ветвей необходимо применять специальные методы. Для этих целей, на втором этапе БА, решается задача локализации, которая выполняет поиск единичных решений других ветвей при фиксированном значении числа Рейнольдса Re = Re0.

Для решения задачи локализации используется модифицированный FPN (Fixed point -Newton) метод гомотопии. Оригинальный метод был предложен в работе [2] для поиска всех вещественных решений одного уравнения в задачах химической инженерии. Реализацию аналогичного метода для системы из нескольких уравнений можно найти и в работе [4], где она рассматривается в общей теории БА как метод разделения ветвей в окрестности точек ветвления (точек пересечения ветвей). Модификация метода для исходной задачи основана на выборочном применении функции гомотопии к уравнениям исходных систем.

Метод FPN гомотопии позволяет находить решения задачи, начиная с одного начального приближения и выполняя последовательное отображение одного решения в другое. Функция гомотопии, реализующая такое отображение для системы (3)-(5), строится следующим образом

F (x,Re0 ) = t -1 G (x) = 0

Здесь оператор F содержит уравнения движения (3а) и граничные условия (4), (5), а оператор G - уравнение связи между завихренностью и функцией тока (3б). Таким образом, построенное отображение (12) не нарушает закона сохранения массы, и удовлетворяет уравнению связи между функцией тока и завихренностью.

Система (12) является однопараметрической нелинейной системой, зависящей от параметра гомотопии t. При t=1 её решения удовлетворяют исходной задаче, среди которых есть решение основной ветви, известное для заданного Re0 . Для перехода к другим решениям рассмотрим задачу продолжения по параметру гомотопии. Представим задачу (12) в форме (6), где оператор P содержит все уравнения системы и сс = t. Используя в качестве начального приближения у0 решение основной ветви, построим ветвь задачи (12) по алгоритму из

раздела 2 и выделим все решения при t= 1.

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

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

4. Исследование точек бифуркаций

Последний этап БА заключается в стратификации фазового портрета исследуемой системы относительно точек бифуркаций. Как было отмечено ранее, основными локальными бифуркациями однопараметрической системы являются бифуркация складки и бифуркация Андронова-Хопфа. Точки бифуркаций первого типа приближено находятся на втором этапе БА, так как соответствуют точкам экстремума ветви относительно параметра системы. В терминах раздела 2 это означает, что da компонента касательного вектора dy в этих точках меняет свой знак. Более точная локализация точек бифуркации складки, а также локализация точек бифуркации Андронова-Хопфа проводится при решении спектральной задачи.

Спектральная задача для решения () системы (3)-(5) записывается в виде:

Яан + Яе

+ Яе

е 2 г (5,) к¥0(д г) й (-£) - (5 г) У0(дв) е

е2 г (д,) (д г) н (4) - (д г) н^н (дв)

е 2 2

+

(д г )Н (4) + (д,) 2 4

= -(А г,в\у/н ■

а

Здесь 1 - собственные числа, (¥н,4) -собственный вектор. Точка бифуркации складки соответствует переходу через ноль вещественного собственного числа спектральной задачи, точка бифуркация Андронова-Хопфа - переходу через ноль вещественной части у пары комплексно-сопряженных собственных чисел. Так как задача исследуется в области умеренных чисел Рейнольдса, то для идентификации обеих бифуркаций достаточно найти несколько собственных чисел с наибольшей вещественной частью. Для этого в работе используется метод ортогональных проекций Арнольди, реализованный в пакете АЯРАСК [13].

5. Результаты

Задача о вторичном стационарном течении вокруг осциллирующего цилиндра рассматривается в диапазоне [0,150] чисел Рейнольдса. Выбранный диапазон охватывает интервал чисел Рейнольдса,

который обычно используется при исследовании основной ветви решения 1<Яе<100 [14, 15], и включает большие числа Рейнольдса, что объясняется повышенным интересом к дополнительным ветвям. Сеточная задача записывается в области г е [0,3.4] (в физической области соответствует г е И30]) на сетке с числом узлов п = 128 х128. Используемое преобразование координат г = 1п г приводит к сгущению узлов вокруг цилиндра в физической области, что позволяет разрешать характерные особенности течения во всем рассматриваемом диапазоне на выбранной сетке.

Локализация решений задачи проводится при Яе = 144. Основные значения параметров метода продолжения, используемые при решении задачи локализации, приведены в следующей таблице 1.

Таблица 1 - Значение параметров метода продолжения для задачи локализации

параметр Да Птр п %тгез п т^ £

шт 0.00002 3 10 4 10-10

тах 0.002 6 20 7 10-10

Здесь Да - шаг по параметру гомотопии между двумя точками аппроксимирующей последовательности, птр - число итераций метода Мура-Пенроуза, необходимых для вычисления каждой точки последовательности, п^тгех -количество итераций метода вМКЕЗ, используемых для решения линейных задач на каждой итерации метода Мура-Пенроуза, ^ - количество итераций многосеточного предобуславливателя на каждой итерации ОМЯБ8, £ - допустимая норма невязки. Для каждого из этих параметров приведено максимальное и минимальное используемое значение. Наиболее критический с точки зрения времени вычисления набор параметров используется в окрестности точек бифуркации складки. Для построения диаграммы задачи локализации требуется найти более 3000 точек аппроксимирующей последовательности.

Задача локализации позволяет построить 9 различных решений задачи. Бифуркационная диаграмма задачи представлена на рис. 1; по вертикальной оси отложена норма вектора

неизвестных

( 2п \ 2

ъ

/(2п -1)

а

по

горизонтальной оси параметр гомотопии - t .

Бифуркационная диаграмма симметрична относительно оси t = 1. Найденные решения исходной задачи выделены на диаграмме квадратными маркерами и имеют следующие обозначения: точка соответствующая основному решению задачи обозначена как 50, в остальных точках пересечения кривой гомотопии с прямой t = 1 содержится по два дополнительных решения; сверху вниз: 5„ и 521 , 512 и $22 , 531 и ^4! , 5*32 и

е

X =

1=1

S42 • Картины течений локализуемых решений представлены на рис. 2.

Рис. 1 - Бифуркационная диаграмма задачи локализации при Яе=144

Все представленные решения являются я -периодичными. Каждое дополнительное решение имеет антисимметричную (по функции тока) пару. Парные решения отображаются в одной точке на бифуркационной диаграмме, так как имеют одинаковое значение нормы ||*||.

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

5ц ¡¡и

Рис. 2 - Картины течения при Яе=144

Ко втором типу относятся решения с индексами 11, 12, 21, 22. У этого типа выделяются

три зоны течения: центральная, содержащая два вихря в плоскости течения, и две боковых, в которых наблюдается циркуляция около бесконечно удаленных вихрей. Режимы 11, 12 и 21, 22 отличаются размером и формой вихрей в центральной области. Решения с индексами 31, 32, 41, 42 относятся к третьему типу. Как и в предыдущем случае, здесь в области течения расположены два вихря, но в отличие от течений второго типа они находятся в непосредственной близости от цилиндра. Для режимов 31,32 и 41,42 эти вихри имеют различную мощность. Вдали цилиндра имеет место циркуляционное течение с затухающей интенсивностью.

Для исследования эволюции найденных решений во всем рассматриваемом диапазоне решаются задачи продолжения. Основные значения параметров для задач аналогичны приведенным в таблице 1, за исключением шага Аа, который для данного случая варьируется в следующих границах [0.5, 10]. Учитывая относительно небольшой диапазон изменения числа Рейнольдса, для аппроксимации ветвей в фазово-параметрическом пространстве требуется построить на порядок меньше точек, чем для задачи локализации.

Бифуркационная диаграмма фазово-параметрического пространства изображена на рис 3; по вертикальной оси отложена норма вектора неизвестных - ||*||, а по горизонтальной оси число Рейнольдса - Яе.

_1_I_1_I_I_I__

20 40 60 80 100 120 140 Яе

Рис 3 - Бифуркационная диаграмма фазово-параметрического пространства

На диаграмме представлены 5 ветвей, антисимметричные дополнительные ветви, содержащие решения с индексами 11, 12 и 21, 22, а также 31, 32 и 41, 42 в указанной параметризации совпадают. Основная ветвь не имеет точек пересечения с дополнительными ветвями в фазово-параметрическом пространстве, в свою очередь найденные дополнительные ветви в рассматриваемом диапазоне также не пересекаются. Первые из найденных дополнительных решений задачи появляются при Яе«42 в результате бифуркаций складки (отмечены круглыми

маркерами). При этом однотипные по структуре решения принадлежат одним ветвям. Третья и четвертая бифуркация складки происходят при Яе«92, в результате образуя ещё две дополнительные ветви, содержащие четыре различных решения.

Исследование спектральной задачи для выявленных решений показало наличие собственных чисел с положительной вещественной частью для всех решений дополнительных ветвей и для решений основной ветви при Яе > 7.2. Критическое значения числа Рейнольдса Яе « 7.2 на основной ветви отвечает бифуркации Андронова-Хопфа. Возникновение бифуркаций складки при Яе « 42, Яе « 92 связано с переходом через ноль одного из вещественных собственных чисел спектральной задачи. Следует подчеркнуть, что наличие собственных чисел с положительной вещественной частью не связано непосредственно с физической неустойчивостью нестационарных течений вокруг цилиндра, осциллирующего в вязкой жидкости. Это связано с тем, что изучаемые в этой работе вторичные стационарные течения имеют следующий порядок малости [6] по сравнению с основным осцилляционным потоком при высоких частотах колебания цилиндра. В частности, основная ветвь вторичного течения в экспериментах наблюдается до значений Рейнольдса порядка 100 [13,14]. При столь больших числах Рейнольдса решение спектральной задачи для основной ветви дает несколько собственных чисел с положительной вещественной частью. Поэтому «неустойчивость» выявленных в данной работе дополнительных ветвей решений не говорит еще об их физической нереализуемости. Этот вопрос нуждается в

дальнейшем изучении, хотя уже сейчас можно сказать, что структуры вторичных течений с индексами 41, 42 в большой степени подобны экспериментально наблюдаемым в режиме «F» работы [15].

Литеаратура

1. K. Johannsen, Computational Geosciences, 7, 169-182 (2003).

2. S. K. Rahimian, et al., Computers Chemical Engineering, 35, 3, 403-411 (2011).

3. А.Н. Нуриев, Вестн. Каз. технологического ун-та., 16, 334-336 (2011).

4. Y.A. Kuznetsov, Elements of applied bifurcation theory. Springer, Berlin, 1995, 591 p.

5. N. Riley, Annual Review of Fluid Mechanics., 33, 43-65 (2001).

6. N. Riley, IMA Journal of Applied Mathematics, 3, 4, 419434 (1967).

7. C.-Y. Wang, J. FluidMech, 32. 55-68 (1968).

8. Г. Шлихтинг, Теория пограничного слоя. Наука, Москва, 1974, 712 с.

9. К. Флетчер, Вычислительные методы в динамике жидкостей: Т.2. Мир, Москва, 1991, 552 с.

10. R. Barrett et al. Templates for the solution of linear systems: building blocks for iterative methods. Society for Industrial Mathematics, 1987, 143 p.

11. J. Ju, G. Lapenta, Lecture Notes in Computer Science, 3514/2005, 601-666 (2005).

12. А.Г. Егоров, А.Н. Нуриев, Учeн. зап. Каз. гос. ун-та. Сер. Физ.-матем. Науки, 151, 3, 130-143 (2009).

13. R.B. Lehoucq et al. ARPACK users' guide - solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods. SIAM. Software, environments, tools, 1998, 142 p.

14. A.F. Bertelsen, J. Fluid Mech., 64, 3, 589-697 (1974).

15. M. Tatsuno, P.W. Bearman, J. Fluid Mech, 211, 157-182 (1990).

© А. Н. Нуриев - асс. каф. информатики и прикладной математики КНИТУ, асп. каф. аэрогидромеханики КФУ, 2апипеу@к8и.га; А. Г. Егоров - проф., зав. каф. аэрогидромеханики КФУ, [email protected].