CC BY
91
УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА Том 151, кн. 3 Физико-математические пауки 2009
УДК 532.517.3
НЕЕДИНСТВЕННОСТЬ СТАЦИОНАРНОГО ТЕЧЕНИЯ ВЯЗКОЙ ЖИДКОСТИ В КВАДРАТНОЙ КАВЕРНЕ
А.Г. Егоров, А.Н. Нуриев
Аннотация
Рассмотрена классическая тестовая задача вычислительной гидромеханики о циркуляционном течении вязкой несжимаемой жидкости в квадратной каверне. Предложен эффективный алгоритм решения уравнений Навье Стокса, позволивший получить стационарные картины течения на очень подробных (до 107 узлов) сетках, для больших (до 10Б) чисел Рейнольдса (йе). Показано, что при больших йе решение задачи не единственно. Особое внимание уделено анализу основной и одной из дополнительных ветвей решения, начинающейся при относительно малых (~ 14000) числах Рейнольдса.
Ключевые слова: уравнение Навье Стокса, стационарное решение, квадратная каверна, мпогосеточпый метод, неединственность, устойчивость.
Введение
Структура двумерного стационарного течения вязкой несжимаемой жидкости в квадратной каверне с подвижной границей при больших числах Рейнольдса является предметом всестороннего изучения в последние сорок лет. Вряд ли этот интерес стоит объяснять лишь практической значимостью: довольно давно из физических и численных экспериментов [1. 2] известно, что двумерное течение в каверне теряет устойчивость по третьему измерению уже при умеренных числах Рейнольдса (Ие < 1000). Научная значимость проблемы определяется двумя обстоятельствами. С одной стороны, решения данной модельной задачи содержат важнейшие структурные элементы вязких течений такие, как сдвиговый поток, нарастающие на стенках пограничные слои, первичный, присоединенный и вторичные вихри. С другой стороны, разномасштабность и сложность этих структур предъявляют растущие с ростом Т1с требования к качеству применяемых численных схем и размерности используемых сеток. Последнее делает задачу о стационарном циркуляционном течении в каверне классическим объектом для тестирования новых методов решения уравнений Навье Стокса. Существование решения данной задачи при любых числах Рейнольдса гарантируется известной теоремой О.А. Ладыженской [3]. Вопросы единственности и устойчивости решения при больших Т1с остаются открытыми.
К сожалению, в огромной массе литературы, посвященной задаче течения в каверне, вопрос об отыскании нескольких ветвей ее решения, насколько нам известно. не ставился. Это объясняется тем. что тестовые расчеты в большинстве работ проводились при некритично больших числах Рейнольдса (Ие < 104), когда различные методы на достаточно подробных сетках приводят к одному и тому же решению. Отметим в этой связи эталонные результаты Гиа [4]. исследовавшего задачу о течении в каверне вплоть до Т1е = 104 на сетках с максимальным числом узлов п = 2572 и обнаружившего единственную стационарную ветвь решения. В дальнейшем будем ссылаться на нее как на основную ветвь. Лишь небольшое
число авторов рассматривало более широкий диапазон чисел Рейнольдса. Так. Бар-раги и Кэри [5]. использовав в своих расчетах высокоточную конечноэлементную схему с числом узлов n = 2Б72, продолжили основную ветвь до Re = 12500. В [6] Эртурк проводил расчеты со вторым порядком точности на равномерных сетках с числом узлов n = 6002. Ему удалось продолжить основную ветвь вплоть до Re = = 20000. В следующей работе [7] Эртурк подтвердил свои результаты используя компактную схему высокого порядка.
Особо выделим работы С.А. Исаева с соавторами [8. 9]. в которых с помощью специализированного (VP2/3) и универсального (FLUENT) гидродинамических пакетов анализируется течение в квадратной каверне с подвижной границей при высоких числах Рейнольдса. Расчеты проводились на структурированных и неструктурированных сетках со сгущением к стенкам каверны: максимальное число узлов составило n = 4002. В [8] авторы исследовали диапазон Re < 40000, а
= 60000
< 20000
с другом, так и с основной ветвыо стационарного решения. С ростом числа Re. однако, наблюдается некоторое расхождение ветвей. Оно выражается, например, в деталях картин течения в правом нижнем углу каверны (рис. 4к [8] и рис. 8е [9]) и
= 40000
Упомянутые различия все же не столь велики, чтобы их определенно нельзя было списать на погрешность расчетов. Видимо, по этой причине авторы не акцентируют на этом внимание. Тем не менее вопрос о том. не наблюдались ли в [8] и [9] различные ветви стационарного решения, остается открытым. Если это так. то. насколько нам известно, это первое наблюдение такого рода.
Основной целыо настоящей работы является установление факта неединственности стационарного решения задачи о циркуляционном течении в каверне с подвижной границей. Указанная цель требует проведения многочисленных расчетов на подробных сетках при больших (десятки тысяч) числах Рейнольдса. Это. в свою очередь, предполагает использование высокоэффективных численных схем решения уравнений Навье Стокса. Очевидно, что использование метода установления, опирающегося на решение нестационарной задачи, недопустимо в силу малой скорости сходимости и, что более важно, неустойчивости искомых стационарных решений. Численная схема должна опираться на стационарную формулировку рассматриваемой задачи. Такая итерационная схема решения стационарной системы Навье Стокса со вторым порядком точности в терминах скорость давление представлена в первых двух разделах статьи. Высокая скорость сходимости, точность и устойчивость этой схемы обеспечиваются двумя факторами. Во-первых, специальная дискретизация конвективных слагаемых позволяет гарантировать диагональное преобладание матрицы линеаризованной системы на каждом итерационном шаге, одновременно обеспечивая второй порядок точности схемы. Во-вторых, использование одного цикла многосеточного метода при решении промежуточных линейных задач позволяет прийти к желаемому компромиссу между необходимой точностью и скоростью вычислений.
Использование в разделе 3 разработанной численной схемы для нахождения основной ветви стационарного решения позволило провести систематические расчеты для Re < 80000 та очень подробных сетках с числом узлов до n = 40962. Сравнение полученных результатов с имеющимися результатами других авторов подтвердили высокую точность и эффективность предложенной схемы. Далее, в разделе 4 с помощью специального алгоритма поиска обнаруживаются несколько других ветвей решения рассматриваемой задачи. Особое внимание уделяется анализу ветви, появляющейся при минимальных (Re ~ 14000) числах Рейнольдса.
О х
О х
х о
х о
х О
Рис. 1. Разнесенное расположение узлов. Жирной линией выделен блок сетки, в точках х определено давление, в о и • - скороети и и V соответственно
1. Постановка и дискретизация задачи
Стационарная система уравнений Навье Стокса, описывающая двухмерное установившееся движение вязкой несжимаемой жидкости, является основой при математическом моделировании гидродинамических процессов. Обозначая через и, V, р горизонтальную и вертикальную компоненты безразмерной скорости и безразмерное давление, запишем ее в виде
иих + vuy + рх — Ке-1Ди = 0,
и^х + vvy + ру — Re-1Дv = 0, (1)
их + Vy = 0.
Областью течения является квадратная каверна П = [0,1] х [0,1] с границей дП, состоящей из твердых стенок, верхняя из которых движется с постоянной скоростью и=1
дПх : у =1: и = 1, V = 0;
дПх : у = 0 : и = 0, V = 0; (2)
дПу : х = 0,1 : и = 0, V = 0.
Для дискретизации задачи (1), (2) в области П строится равномерная сетка с
шагом Н = п-1. Используется разнесенное расположение узлов, когда дискретные значения давления определяются в центрах ячеек сетки, а дискретные составляю-
дП
через точки, в которых определены нормальные к границе составляющие скорости. При этом все узлы, в которых определено давление, лежат внутри расчетной области.
При разнесенном расположении узлов дискретизация уравнения неразрывности проводится в центрах ячеек. Производные в третьем из уравнений (1) аппроксимируются со вторым порядком точности центральными разностями. Используя для них обозначения (дх)|!с\ (ду)^ , запишем дискретный аналог уравнения неразрывности в виде
(дх)1с) ин + (ду^ ^ = 0. (3)
Дискретные аналоги уравнений сохранения количества движения записываются в точках, помеченных на рис. 1 символами о и • (о-узды и •-узлы). Для дискретизации со вторым порядком точности вязких слагаемых, градиента давления и конвективных слагаемых используются стандартная пятиточечная дискретизация оператора Лапласа и центральные разности:
ил, (дх)^с) и^ + ^ (ду)^с) и^ — Ке-1Д^и^ + (дх)^с) рл, = 0,
(4)
ил, (дх)!с) ^ ^ (ду)|!с) ^ — Ие-1Д^л + (ду)|!с) рл = 0.
Обратим внимание на то, что вертикальная скорость ^ в первом из уравнений (4) должна быть определена в о-узлах, в то врет как она известна лишь в соседних • -узлах. Для ее определения с точпостью О (Н2)в о-узле применяется осреднение
•
ляция для и ^ используется и во втором из уравнений (4).
На каждой из границ расчетной области необходимо удовлетворить два граничных условия (2). Одно из них равенство нулю нормальной к границе компоненты скорости удовлетворяется точно. Для удовлетворения второго равенства нулю тангенциальной компоненты скорости используются граничные условия Тома [13], которые, как известно [14], согласованы по точности О(Н2) с представленной выше дискретизацией уравнений.
Система (3), (4), как и исходная система уравнений Навье Стокса, существенно нелинейна. Поэтому для ее решения необходимо применять специально сконструированные итерационные методы.
2. Итерационный процесс
Для решения нелинейной системы (3), (4) строится итерационная процедура с решением на каждом шаге линейной задачи. Так как нелинейность рассматриваемой системы заключена в конвективных членах уравнений (4), то проблема заключается в выборе способа их линеаризации. Предлагается провести ее следующим образом:
(5)
(«і (дж)1“р) «і + V0 (ду)|г“р) «і) - Ие + (дх)^ Рі =
«і (дх)1“Р) «і + (ду)І“Р) <) - (< (дх)1С) «і + (ду)ІС) «і
і (дх)і“Р) «і + (ду)І“Р) «і) - Ие-1Д 1«1 + (ду)!С) Рі =
= (и і (дх)і“Р) V0 + V0 (ду^ V0) - (и і (дх)1С) V0 + V0 (ду)1С) V0) .
Здесь верхними индексами «0» отмечены переменные, известные с предыдущей итерации, переменные без индексов относятся к текущей итерации. Верхний индекс «ир» указывает на использование направленных разностей при вычислении соответствующих дискретных производных. Применение направленных разностей в (5) обеспечивает диагональное преобладание матрицы линеаризованной системы уравнений.
Используемая итерационная процедура состоит из следующих шагов.
Шаг 0. Выбирается начальное приближение (и0,«0,р0).
Шаг 1. Решается линейная задача вида
Ь(«”+1/2,-у”+1/2,р”+1; Ие) = /,
где Ь - матрица коэффициентов линеаризованной системы (3), (5),
(и"+1/2, ®"+1/2,р"+1) - вектор неизвестных, / - вектор, состоящий из правых частей уравнений (5) и подсчитываемый по и",®".
Шаг 2. Вычисляется новое приближение для скоростей по формулам и"+1 = аи"+1/2 + (1 - а)и", ®"+1 = а®"+1/2 + (1 - а)®",
где а - коэффициент релаксации (0 < а < 1);
Шаг 3. Рассчитывается значение невязки в уравнениях (3), (4). Если норма невязки г больше наперед заданного значения е0, то осуществляется переход к шагу 1.
Наиболее трудоемким шагом данной процедуры является решение линейной системы уравнений (3), (5). Учитывая, что его необходимо проводить на каждой итерации, эффективность предложенного подхода можно было бы поставить под сомнение. Однако, как показывают расчеты, для сходимости итерационного процесса достаточно найти лишь приближенное решение линейной задачи. Оно строится с помощью одной итерации геометрического многосеточного метода [15].
Конкретизируем базовые операции многосеточного метода. Для сглаживания использовалась так называемая «блочная релаксация» [15], когда на каждом релаксационном шаге модифицируются 5 узловых значений, входящих в данный блок сетки: давление в центре и четыре значения скорости на гранях блока. Для этого решается линейная система пяти уравнений с неизвестными модифицируемыми значениями. При проходе по всем блокам, очевидно, каждая дискретная компонента скорости обновляется дважды, а дискретные давления в центрах ячеек по одному разу. Релаксационный коэффициент сглаживателя принимался равным
0.75, что является теоретически оптимальным выбором при решении уравнения Стокса [15]. Операция сужения, используемая для переноса невязки с точной на грубую сетку, проводилась по естественным шаблонам:
1 2 1 'V =1 г2 к = 8 1 1 = 1 '2 к = 4 1 1
о 'к, 2 • 2 ' к, х
1 2 1 8 1 1 4 1 1
Здесь индексы и, ®, р указывают на уравнения сохранения импульса и массы, записанные в точках о, •, х соответственно. Аналогичным образом сужались коэффициенты ик, ®0 в уравнении (5). Обратная операция продолжения с грубой сетки на точную проводилась с помощью билинейной интерполяции.
Помимо базовых операций метода необходимо выбрать структуру многосеточного цикла и определить оптимальное количество релаксаций при его прямом и обратном ходе. Тестировались два наиболее популярных цикла: V и \¥ (см. [15]). Численные эксперименты показали, что для рассматриваемой задачи предпочтителен \У-цикл. Уже при низких числах Рейнольдса скорость сходимости численной схемы с У-циклом в четыре раза ниже, чем у аналогичной схемы с ^У-циклом, а при умеренных числах Рейнольдса схема с У-циклом вообще не сходится к стационарному решению. Это является следствием того, что более «дешевый» У-цикл не позволяет за одну итерацию добиться достаточно хорошего приближения к точному решению задачи (3), (5).
Модификация \У(1,1) многосеточного метода \У-цикл с одной релаксацией при прямом и обратном ходе показала высокую эффективность в диапазоне чисел Рейнольдса 0 < Т1е < 25000. Здесь она сходится по времени СРи почти вдвое
Рис. 2. Зависимость оптимального релаксационного параметра а от йе
г
Рис. 3. Зависимость невязки г от числа итераций к при йе = 35000
быстрое, чем альтернативная модификация \¥(2,2) с двумя релаксациями на прямом и обратном ходе. Однако в диапазоне высоких чисел Рейнольдса 30000 < Т1е <
< 80000 устойчивую сходимость к стационарному решению демонстрирует лишь схема "\¥(2,2). Именно она и выбиралась в качестве основной в проводимых расчетах. Отметим, что смена порядка обхода ячеек при первой и второй релаксациях в Л¥(2,2) может существенно ускорить сходимость итерационного процесса. Нами использовался лексикографический порядок при первой релаксации (слева направо, сверху вниз) и обратный (сверху вниз, слева направо) при второй релаксации.
Важным параметром, обеспечивающим сходимость рассматриваемой итерационной процедуры, является внешний релаксационный параметр а. Если его выбрать слишком малым, то скорость сходимости итерационного процесса снизится, если слишком большим, то процесс разойдется. Расчеты показывают, что опти-а
образом числом Рейнольдса. На рис. 2 изображена зависимость оптимального внешнего релаксационного параметра от числа Рейнольдса, построенная по результатам расчета на сетке 1024 х 1024. Вертикальные отрезки па этом рисунке
а а
диапазоне своего изменения не является критичным.
3. Основная ветвь решения
Для малых чисел Рейнольдса решение рассматриваемой задачи единственно, и при любом выборе начального приближения описанная итерационная процедура сходится к этому решению. Для больших И.о сходимость численной схемы можно гарантировать только при специальном выборе начального приближения, которое должно находиться в «окрестности» точного решения. В расчетах основной ветви это достигалось методом «последовательных нагружений», когда полученное для меньших чисел Рейнольдса решение использовалось в качестве начального приближения для больших чисел Рейнольдса. Обычно шаг нагружения составлял ДИе = 5000 . При попадании в «окрестность» скорость сходимости от конкретного выбора начального приближения практически не зависела. Это проиллюстрировано на рис. 3, где представлена история сходимости итерационного процесса на сетке 1024 х 1024 = 35000
с ДИе = 15000, 10000, 5000, а норма невязки определяется как максимальная из трех поточечных невязок в уравнениях (3), (5).
Рис. 4. Профили скоростей в средних сечепиях каверны. Кривые слева отвечают значениям йе = 1000, 5000, 50000. Ось абсцисс: х и (и + 1)/2; ось ординат: у и (V +1)/2. Кривые справа построены для йе = 50000 на различных сетках
Табл. 1
Число итераций и время расчета при использовании численной схемы с оптимальным выбором параметров
Re а Число итераций
n = 512 n = 1G24 n = 2G48 n = 4G96
5000 0.9 89 68 59 59
20000 0.8 399 274 240 245
50000 0.55 1094 700 523 479
80000 0.4 2650 1575 1255 1030
Время 1-й итерации, с 0.61 2.6 12 74
Эффективность предложенного алгоритма решения рассматриваемой задачи можно оценить по данным табл. 3. для одного шага нагружения с ARe = 5000. Расчеты проводились до достижения невязкой величины 10-9 на стандартном персональном компьютере на базе процессора Coro(TM)2 Duo Processor Е6550. Как и следовало ожидать, с ростом Re число итераций и время расчета растет. Тем не менее время счета остается приемлемым даже для предельно больших чисел Re ^ ~ 105. Это позволяет эффективно использовать предложенный алгоритм расчета в широком диапазоне изменения чисел Рейнольдса.
Подчеркнем принципиальную важность использования в расчетах подробных сеток. Они должны позволять разрешать топкие пограничные слон на стенках каверны н малые структурные элементы течения в ее углах. Использование «грубых» сеток (n = 128, 256) приводит при больших Re к расходимости используемой итерационной процедуры. В то же время на более подробных сетках итерационная процедура сходится за меньшее число итераций (см. табл. 1). На рис. 4 показаны вертикальная компонента скорости в среднем горизонтальном сечении и горизонтальная компонента в среднем вертикальном сечении каверны. Видно, что при больших числах Рейнольдса на стенках каверны формируются очень тонкие пограничные слои, которые, тем не менее, вполне разрешаются на сетках использованных в расчетах .
(а)
(Ь)
(с)
= 1000
Ь - йе = 5000, с - йе = 10000, (1 - йе = 20000, е - йе = 40000, £ - йе = 80000
Наблюдаемую в каверне картину течения удобно представлять с помощью линий тока. Для различных чисел Рейнольдса они приведены на рис. 5. Гидродинамический анализ представленной на этом рисунке эволюции картины циркуляционного течения в каверне неоднократно проводился различными авторами (см., например, работу [9] и цитированную в ней литературу). Не будем останавливаться на этом анализе подробно, тем более, что наши наблюдения и гидродинамические выводы вполне согласуются со сделанными в [9]. Отмстим лишь, что указанные в [9] для Т1е « 60000 тенденции эволюции картины течения сохраняются и при дальнейшем росте числа Рейнольдса. Основная из них состоит в растяжении угловых вторичных вихрей, их пережатии и образования на их месте каскада мелкомасштабных вихревых ядер.
В табл. 2 приведены основные экстремальные характеристики полученных решений при Т1е ^ 20000. Помимо максимального значения функции тока фтах -величины, традиционно используемой для описания мощности первичного вихря, в таблице представлены также суммарная кинетическая энергия течения
К =
1
(и
I йх<іу,
максимальная и минимальная вертикальные скорости «тах, ут
и минимальная
горизонтальная скорость мт;п = тт и(0.5, у) в среднем сечении каверны. Аналогичные таблицы для Т1е ^ 20000 приведены в [6, 7, 9]. Полученные нами данные отличаются от результатов [6. 7] не более чем в четвертой значащей цифре.
Представленные в табл. 2 результаты получены на наиболее подробной из использованных сеток с шагом Н = 4096-1. В скобках показано, как изменяются
Табл. 2
Экстремальные характеристики базового решения при больших И.е
И.е ^шах К Ушіп ^шах ишіп
2 • 104 0.12212(19) 0.1218 0.04764(9) -0.7722(34) -0.7733 0.4773(6) 0.4767 -0.4638(42) -0.4644
5 4 3 111 ООО ^ ^ ^ 0.12185(96) 0.12160(75) 0.12138(96) 0.04750(8) 0.04735(47) 0.04722(38) -0.7843(60) -0.7924(48) -0.7983(8028) 0.4827(32) 0.4859(65) 0.4880(8) -0.4637(41) -0.4622(8) -0.4603(10)
6 • 104 0.12120(43) 0.1233 0.04710(29) -0.8028(86) -0.8287 0.4897(906) 0.5312 -0.4581(91) -0.5015
8 7 11 О О 0.12103(30) 0.12087(120) 0.04700(22) 0.04691(716) -0.8056(134) -0.8093(186) 0.4909(21) 0.4919(33) -0.4561(72) -0.4540(53)
последние значащие цифры при проведении экстраполяции по схеме Ричардсона [16] па случай Н = 0. Схема Ричардсона основывается на представлении искомой величины (например, ^тах) в виде
^тах(Н) = <ах + с1н2 + С2Н4 + • • • .
Выполняя расчеты па трех различных сетках, можно найти константы С1, С2 и, что самое важное, определить предельное при Н ^ 0 значение ^тах- Данные табл. 2 получены при использовании сеток с Н = 4096-1, Н = 2048-1, Н = 1024-1. Аналогичная экстраполяция выполнялась и по сеткам с Н = 2048-1, Н = 1024 1, Н = 512-1. Результаты двух экстраполяций отличались не более чем в пятой значащей цифре, что дополнительно свидетельствует о высокой точности проведенных расчетов.
Знаменатели в табл. 2 соответствуют данным, полученным в работе [9]. Как видно, при Ие = 20000 они хорошо согласуются с нашими расчетами. Однако с ростом числа Рейнольдса наблюдается расхождение, достигающее 8% для «тах при И.е = 60000. Различен характер поведения величины ^тах с ростом Т1е: она убывает в наших расчетах и возрастает согласно [9]. Отличны также и детали картины течения в правом иижиом угле каверны (см., например, рис. 5о и рис. 8о [9]). Это, на наш взгляд, может свидетельствовать о том, что в настоящей работе и в работе [9] получены различные ветви стационарного решения рассматриваемой задачи при больших числах Рейнольдса. В этой связи интересно отметить, что в [9] использовалась отличная от нашей стратегия расчетов. Вместо нагружения при фиксированной сетке по числу Рейнольдса проводилось нагружение по размеру сотки при фиксированном И.о. Именно этот прием мы используем в следующем разделе настоящей статьи для того, чтобы идентифицировать дополнительные ветви стационарного решения.
В заключение этого раздела остановимся на вопросе линейной устойчивости базового стационарного решения по отношению к двумерным возмущениям. На основе имеющихся в литературе данных можно утверждать, что основная ветвь теряет свою устойчивость с ростом числа Рейнольдса в результате бифуркации Хопфа с переходом от стационарного к периодическому режиму течения. Имеющиеся оценки критического числа Рейнольдса Т1с(.г лежат в достаточно широком
диапазоне 6500 < Recr < 9000. Результаты последних лет существенно сужают его. Так. в работе [10] на основе метода Дженнингса получено и прямым численным моделированием подтверждено значение Recr = 8000. Расчеты, проведенные в [11] с использованием метода Арнольди, дали значение Recr = 8030. В [12] использовался метод собственных ортогональных проекций (POD): критическое число Рейнольдса оценено величиной Recr = 7819. Наши расчеты [17], как и в [11], опирались на метод ортогональных проекций Арнольди для разреженных матриц, реализованный в пакете ARPACK [18]. Они проводились па сетке с числом узлов n = 5122
= 8051
=
T = 2.25, что вполне согласуется с оценкой T = 2.21, данной в [10].
4. Дополнительные ветви решения
Помимо основной ветви стационарного решения в задаче о циркуляционном течении в каверне при больших числах Рейнольдса (Re > 14000) удается обнаружить и другие, дополнительные, ветви. Они оказываются линейно неустойчивыми, как, впрочем, и основная ветвь в этом диапазоне чисел Re. Неустойчивость дополнительных ветвей делает их расчет трудоемкой задачей. К тому же число ветвей растет с ростом числа Рейнольдса. По этим причинам мы ограничиваемся в расчетах сеткой с числом узлов 512 х 512 и диапазоном Re < 20000, в котором используемая сетка позволяет уверенно разрешать основные структурные элементы течения.
Построение дополнительной ветви включает в себя две подзадачи. Первая из них заключается в нахождении отличных от основного решений задачи при заданном числе Рейнольдса. Вторая состоит в отслеживании эволюции этих решений во всем диапазоне чисел Рейнольдса.
Первая подзадача решалась на основе процедуры нагружения по размеру сетки при фиксированном числе Рейнольдса. Выбиралось максимальное из расчетных = 20000 64 х 64
полной картины течения. Однако их использование позволяло с помощью описанной в разделе 3 итерационной процедуры при малых значениях параметров релаксации а = 10-1 ^ 10-2 обеспечить сходимость в зависимости от выбора начального приближения к нескольким стационарным решениям задачи. Низкая при а
размерности задачи. Обнаруженные решения последовательно интерполировались на более подробные сетки и итерационно уточнялись на них. Построенные в результате данной процедуры решения на подробной сетке определяют единичные точки па дополнительных ветвях бифуркационной диаграммы.
Решение второй подзадачи призвано обеспечить продвижение по дополнительным ветвям с изменением числа Рейнольдса. К сожалению, метод нагружения по числу Рейнольдса, использованный для вычисления главной ветви, оказался в данной ситуации практически непригодным. Это связано с двумя обстоятельствами. Во-первых, сама сходимость к дополнительному решению чрезвычайно чувствительна к качеству начального приближения, что требует использования неприемлемо малых шагов нагружения ARe. При увеличении ARe итерационная процедура сходится к основной ветви решения. Во-вторых, метод нагружения по числу Рейнольдса принципиально неприменим в окрестности точки бифуркации Rcb минимального значения Re, вплоть до которого продолжима дополнительная ветвь. В этой окрестности при каждом Re существует по крайней мере пара близких дополнительных решений.
Обе указанные трудности обходятся с помощью специального метода продолжения (“pseudo-arclength continuation method with secant prediction”), предложенного
для анализа динамических систем [19] и успешно использованного в [20] при построении бифуркационной диаграммы задачи Элдера. Согласно этому методу продвижение по дополнительной ветви реализуется решением расширенной системы дискретных уравнений (3), (4), в которой наряду со скоростями и, V и давлением р неизвестным считается также число Рейнольдса. Дополнительное условие для нахождения Т1ек па к-м шаге метода записывается как
(ий - и&_1, и&_1 - ий_2) + ^ - Vfc_1, Vfc_1 - ^ — 2 ) = 6, (6)
где скобки означают скалярное произведение, 6 - фиксированный малый параметр.
Для решения расширенной задачи (3), (4), (6) использовалась следующая модификация алгоритма, описанного в разделе 3.
Шаг 0. С помощью линейной экстраполяции вычисляется начальное приближение (и£,^,р£,Е е£)
г0 = 2^— - 2^—2, г =(и^,р,К,е).
Шаг 1. Проводится одна итерация многосеточного метода для решения пары линейных задач
Г(и”+1/4 ч,”+1/4 Р”+1/4-Р?<=>”)_ / Г(и”+1/2 ч,”+1/2 Р”+1/2-Р?<=>” + Л)_ /
Мик > Рк ;Кек) = 7 > Мик > Рк +Л) = 7 •
Л
тах величиной порядка единицы. По найденной паре решений приближенно вычисляются производные и, v,p по числу Рейнольдса, после чего значения указанных функций уточняются по формулам
и”+3/4 = и”+1/4 + Т (и”+1/2 - и”+1/4) ,
..Г3/4 = »”+‘/4 + т (<'”+1/2 - <'”+‘/4) ,
”+1 ”+1/4 , ( ”+1/2 ”+1/4\
Р к = Р к + Т(ЧР к - Р к ):
Ее”+1 = 11 е” + Т Л.
Константа т находится из уравнения (6):
е / ”+1/4 \ / ”+1/4 \
= 6 - (и к - ик_1,ик_1 - ик_2 ) - (V к - V к_1, Vк_1 - V к_2 )
/ ”+1/2 ”+1/4 \ . / ”+1/2 ”+1/4 \
(и к - и к ,ик_1 - ик_2 ) + (V к - V,, .Vк_1 - Vk_2 )
Шаг 2. Вычисляется новое приближение для скоростей по формулам
и”+1 = аи”+3/4 + (1 - а)и”, v”+1 = ™”+3/4 + (1 - а>”.
В расчетах коэффициент релаксации а принимался равным 10_1 ^ 10_2.
Шаг 3. Рассчитывается значение невязки в уравнениях (3), (4). Если норма невязки г больше наперед заданного значения е0, то осуществляется переход к шагу 1.
Построение полной бифуркационной диаграммы для рассматриваемой задачи выходит за рамки настоящей статьи, основной целыо которой является демонстрация существования отличных от базового решения. Поэтому ограничимся в дальнейшем представлением единственной дополнительной ветви, начинающейся при относительно небольших числах Рейнольдса. Ее отображение в плоскости Ке^5
Рис. 6. Основная (1) и дополнительная (2) ветви стационарного решения задачи о циркуляционном течении в квадратной каверне
Рис. 7. Картины течения для дополнительной ветви решения. Маркеры а, Ь. с соответствуют отмеченным точкам па бифуркационной диаграмме рис. 6
представлено на рис. 6. Через S обозначено среднее касательное напряжение на нижней стейке каверны:
S = [ — (0, ж) йж.
Дополнительная ветвь образуется при Т1е = Т1еь ~ 13800 в результате бифуркации складки. Общая картина течения для найденной ветви решения представлена па рис. 7. Как видно, в нижней части каверны она существенно отличается от основной картины течения (рис. 5). Бывшие ранее изолированными, вторичные угловые вихри сливаются здесь в единую вихревую систему, что характерно обычно для вытянутой по вертикали прямоугольной каверны [21].
Исследование дополнительной ветви на устойчивость проводится так же. как и исследование основной ветви, и сводится к отысканию собственных чисел соответствующей матрицы. Расчеты показывают наличие нескольких комплексно сопряженных собственных чисел с положительной вещественной частью на протяжении всей дополнительной ветви а—6—с (рис. 6). Это свидетельствует о ее неустойчивости. Наибольшее вещественное собственное число является простым, оно отрица-
тельно на дуге а—6, положительно на Ь—с и равно нулю в точке бифуркации Ь. Это позволяет идентифицировать соответствующую бифуркацию как бифуркацию складки.
Работа выполнена при финансовой поддержке РФФИ (проект X- 08-01-00548-а).
Summary
A.G. Egorov, A.N. Nuriev. Non-uniqueness of a Stationary Viscous Flow in the Square Lid-driven Cavity.
The article considers the classical benchmark problem in computational hydromechanics regarding a viscous incompressible flow in a square lid-driven cavity. An effective algorithm for solving a Navier Stokes system of equations is proposed, that allows to construct stationary solutions on very detailed grids (up to 107 grid points) for large (up to 105) Reynolds numbers. Non-uniqueness of the stationary solution at large Reynolds numbers is shown. Special attention is given to the analysis of the main branch and one of the additional branches of the solution, appearing at relatively small (« 14000) Reynolds numbers.
Key words: Navier Stokes equation, stationary solution, lid-driven cavity, multigrid,
non-uniqueness, stability.
Литература
1. Prasad A.K.,Koseff J.R. Reynolds number and end-wall effects on a lid-driven cavity
flow // Phys. Fluids. 1989. V. 1, No 2. P. 208 218.
2. Shankar P.N., Deshpan.de M.D. Fluid Mechanics in the Driven Cavity // Annu. Rev.
Fluid Mech. 2000. V. 32. P. 93 136.
3. Ладыженская О.А. Математические вопросы динамики вязкой несжимаемой жидкости. М.: Наука, 1970. 205 с.
4. Ghia U., Ghia K.N., Shin С. Т. Higli-Re Solutions for Incompressible Flow Using the Navier-St.okes Equations and a Multigrid Method // J. Comp. Phys. 1982. V. 48.
P. 387 411
5. Barragy E., Carey G.F. Stream Function-Vorticity Driven Cavity Solutions Using p Finite Elements // Computers and Fluids. 1997. V. 26. P. 453 468.
6. Erturk E., Corke T.C., Gokcol C. Numerical Solutions of 2-D Steady Incompressible
Driven Cavity Flow at High Reynolds Numbers // Int. J. Numer. Met.li. Fluids. 2005.
V. 48. P. 747 774.
7. Erturk E., Gokcol C. Fourth Order Compact Formulation of Navier-St.okes Equations and Driven Cavity Flow at. High Reynolds Numbers // Int.. J. Numer. Met.li. Fluids. 2006. V. 50. P. 421 436.
8. Исаев С.А. и др. Моделирование ламинарного циркуляционного течения в квадрат-
ной каверне с подвижной границей при высоких числах Рейнольдса // Ипжеперпо-физ. журп. 2002. Т. 75, Л» 1. С. 55 60.
9. Исаев С.А. и др. Моделирование ламинарного циркуляционного течения в квадратной каверне с подвижной границей при высоких числах Рейнольдса с помощью пакетов VP2/3 и FLUENT // Ипжеперпо-физ. журп. 2005. Т. 78, Л'! 4. С. 163-179.
10. Fortin A. et. al. Localization of Hopf bifurcations in Fluid Flow problems // Int.. J. Numer. Met.li. Fluids. 1997. V. 24. P. 1185 1210.
11. Sahin М., Owens R.G. A novel fully-implicit. finite volume method applied to the lid-driven cavity problem. Part. II. Linear stability analysis // Int.. J. Numer. Met.li. Fluids. 2003. V. 42. P. 79-88.
12. Cazemier W., Verstappen R.W.C.P., Veldman A.E.P. Proper orthogonal decomposition
and low-dimensional models for driven cavity flows // Pliys. Fluids. 1998. V. 10,
No 7. P. 1685 1698.
13. Флетчер К, Вычислительные методы в дшашке жидкостей. Т. 2. М.: Мир. 1991. 552 с.
14. Weinan Е., Liu J.G. Vorticity boundary conditions and related issues for finite difference
schemes // J. Comp. Pliys. 1996. V. 124. P. 368 382.
15. Trottenberg U., Oosterlee C.W., Schuller A. Multigrid. London: Acad. Press, 2001. 631 p.
16. Ruache P.J. Quantification of uncertainty in computational fluid dynamics // Annu. Rev. Fluid. Mech. 1997. V. 29. P. 123 160.
17. Нуриев A.H. Устойчивость плоского стационарного течения в каверне с подвижной крышкой // Труды матем. центра им. Лобачевского. 2008. Т. 37. С. 133 135.
18. Gomes F.M.,Sorensen B.C. ARPACK------------: A С------Implementation of ARPACK
Eigenvalue Package. URL: www.caam.rice.edu/software/ARPACK.
19. Kuznetsov Y.A. Elements of Applied Bifurcation Theory. Berlin: Springer, 1995. 591 p.
20. Johannsen K. On the validity of the Boussinesq approximation for the Elder problem // Comp. Geosci. 2003 V. 7. P. 169 182.
21. Abouhamza A,, Pierre R., A neutral stability curve for incompressible flows in a rectangular driven cavity // Math. Comp. Modelling. 2003. V. 38. P. 141 157.
Поступила в редакцию 03.03.09
Егоров Андрей Геннадьевич доктор физико-математических паук, заведующий кафедрой аэрогидромехапики Казанского государственного университета.
E-mail: Andrey.EgorovQksu.ru
Нуриев Артем Наилевич бакалавр механики, магистрант кафедры аэрогидромехапики Казанского государственного университета.
E-mail: аНет.501 fflist.ru
