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

Вычислительные технологии

Том 18, № 1, 2013

О нестационарных решениях в задачах гидродинамики со стационарными краевыми условиями

Ю.Н. ЗАХАРОВ, К. С. ИВАНОВ Кемеровский государственный университет, Россия e-mail: zyn@kemsu.ru, topspin83@mail.ru

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

Ключевые слова: численное моделирование, вычислительная гидродинамика, уравнения Навье — Стокса, вязкая однородная несжимаемая жидкость.

Введение

В известной статье В.И. Юдовича [1] одна из "великих" проблем современной гидродинамики обозначена как доказательство (или опровержение) глобальной теоремы существования стационарных и вынужденных периодических движений вязкой несжимаемой жидкости. Автор отмечает: "... при увеличении числа Рейнольдса стационарный режим может исчезнуть — уйти на бесконечность, когда Яе ^ Яе*, причем критическое значение Яе* конечно. Однако перед тем как исчезнуть этот стационарный режим теряет устойчивость, порождая автоколебательный периодический режим (возможно, конечно, что сначала происходит ветвление в классе стационарных режимов). Интересно было бы проверить хотя бы в численном эксперименте возможность такого хода событий". Действительно, при решении задач о течении вязкой несжимаемой жидкости во многих случаях используется стационарная постановка исходной системы дифференциальных уравнений Навье — Стокса. Это объясняется тем, что для практических целей зачастую важен не сам процесс развития течения, а его установившийся режим. Для решения стационарных задач гидродинамики многими исследователями в течение последних десятилетий активно разрабатываются различные численные алгоритмы. Например, метод искусственной сжимаемости [2] позволил применить высокоэффективные схемы расщепления [3] для получения стационарных

решений практически важных двух- и трёхмерных задач о течении вязкой несжимаемой жидкости. В действительности же любое течение жидкости является нестационарным и всегда имеет фазы своего начала и развития. Используя стационарную постановку задачи, неявно делается естественное, на первый взгляд, предположение о том, что при использовании стационарных краевых условий течение должно довольно быстро перейти в стационарный режим. В большинстве случаев, как показали численные и натурные эксперименты, это действительно так. Однако стационарная модель может быть применена не к любой задаче о течении вязкой несжимаемой жидкости. Во-первых, нет строгих теорем, которые бы гарантировали установление течения для любых значений чисел Рейнольдса [1]. Во-вторых, существуют известные парадоксы симметрии вязких течений [4] и численные расчёты некоторых нестационарных задач [5], которые демонстрируют отсутствие стационарного режима течения жидкости при стационарных краевых условиях.

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

1. Постановка задачи

Основные уравнения в переменных вихрь — функция тока, описывающие плоское течение однородной вязкой несжимаемой жидкости, имеют следующий вид (все величины обезразмерены):

где V — вектор скорости, И,е — число Рейнольдса, ш — вихрь, ф — функция тока. Будем считать, что течение жидкости происходит в некоторой односвязной области О с границей дО на временном промежутке [0, то]. Для систем (1), (2) ставятся начальные и краевые условия [6, 7].

2. Метод решения

Введём в области О разностную сетку О^. Для численного решения уравнения переноса вихря (1) будем использовать схему стабилизирующей поправки и метод продольно-поперечной прогонки. Аппроксимируя обычном образом на сетке О^ уравнение

(1)

А ф = -ш,

(2)

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

Аи = /,

(3)

где и, / € Ят, А — линейный оператор: Ят ^ Ят. В настоящей работе для решения разностных задач, аппроксимирующих уравнение Пуассона, предлагается использовать итерационный метод, впервые применённый в [8] и затем развитый в [9], который показал высокую эффективность при решении стационарных уравнений Навье — Стокса. Для решения системы (3), следуя [9], построим следующий итерационный процесс:

и

-+1/2 = ип - Тп+1 [Аип - /]

и

п+1

и

п+1/2

и

п+1

ик-1 - У^ <+1 , к = 1

п+1

и

п+1

гeSfc

и-+1, п = 0,1, 2...,

где и0 — произвольный вектор из Ят; тп+1, а-+1 — итерационные параметры; — непересекающиеся подмножества (группы) множества всех индексов такие, что и = Б;

к

;п — векторы из Ят с одной ненулевой г-й компонентой. Оптимальные итерационные параметры тп+1 и а-+1 в (4) определяются из условия минимума норм векторов невязок гп+1/2 и гп+1 соответственно. Как показано в [9], при таком выборе итерационный процесс (4) является сходящимся при любом начальном приближении и произвольном неособенном операторе А. Для определения оптимальных значений ап+1 при каждом к =1, 2,..., I необходимо решить систему линейных алгебраических уравнений вида

(

(А* (А*

V

кРк

А*пх

(А; (А;

'к1

к2

кРк

А*к2 А*к2

> А*к2

кА*п, А*п

\

А*п„,, А;

кРк

( а

х

а.

к1

п+1 \

к2

п+1

Рк

7

\ акрк

кр п+1

( А*;

А;'

п „п+1 к1, 'к

к2, Гк

п+1

А*

рк

гп+1 , 'к

где кг € Бк, г = 1, 2,...,рк, ^ рк = т. Отметим: чем большее количество индексов к

к=1

участвует в группах Бк (т. е. чем больше рк, или меньше /), тем большую эффективность алгоритма (4) следует ожидать. Вместе с тем при увеличении размерности групп Бк увеличивается и размерность системы (5), решение которой в общем случае по степени сложности эквивалентно решению исходной системы (3). Оператор системы (3) не является произвольным, а есть следствие аппроксимации дифференциального оператора Лапласа и, следовательно, его матрица имеет блочно-ленточную структуру. Например, в простейшем случае, когда область расчёта параллелепипед, матрица оператора имеет только семь ненулевых диагоналей. Очевидно, что в силу такой структуры матрицы оператора многие из величин (А;™, А;п) , г,= 1, 2, ...,т, равны нулю. Следовательно, группы Бк всегда можно выбрать так, что для каждого к матрица системы (5) будет иметь диагональный вид и нахождение величин а^+1 не составит труда [9]. Рассмотрим ситуацию, когда имеется всего одна группа индексов, совпадающая с множеством всех индексов Б. Можно показать [9], что в этом случае метод (4) сойдется к точному

к

и

п

решению системы (3) за одну итерацию. Для определения оптимальных итерационных параметров агп+1,% = 1, 2, ...,т, необходимо решить систему (5), размерность которой в данном случае совпадает с размерностью исходной системы. В такой постановке данная задача эквивалентна исходной. Однако в качестве гП могут быть не единичные векторы, а такие, что (А^П, Агп) = 0, г,] = 1, 2, ...,т, % = ], причём матрица системы (5) вновь будет иметь диагональный вид. Иначе говоря, для достижения этой цели необходимо построить ортогональную систему векторов гп, ] = 1, 2, ...,т, в пространстве Я™* а со скалярным произведением (А* Аг, г) , г Е Я™. Задача ортогонализации системы векторов в общем случае эквивалентна по степени сложности задаче решения системы линейных уравнений, однако такой подход в нашем случае обладает рядом значительных преимуществ. Во-первых, в силу упомянутой выше блочно-ленточной структуры матрицы исходного оператора А часть векторов гП уже ортогональны в Я™*а. Во-вторых, если однажды построить указанную ортогональную систему векторов, то возможно её применение и в последующем для решения системы (3) при условии, что оператор А не изменился. В-третьих, можно решить задачу построения векторов, обеспечивающих не диагональный, а, например, трёхдиагональный (не представляющий труда для обращения) вид матрицы системы (5), что значительно экономит вычислительные затраты. Отметим, что при неизменности границы области решения разностный оператор системы уравнений для функции тока не изменяется при переходе от слоя к слою по времени, а также сохраняет структуру своей матрицы при увеличении количества ячеек расчётной области. Это позволяет построить достаточно эффективный алгоритм (4) с помощью оптимального выбора групп Б к, включив в каждую из них максимально возможное число индексов. При проведении расчётов с умеренным количеством узлов разностной сетки очень экономичным оказывается применение ортогональной системы векторов гп, ] = 1, 2, ...,т, обеспечивающей сходимость процесса (4) за одну итерацию к точному решению системы (3). Построив такую ортогональную систему векторов на первом шаге по времени с затратой существенного количества ресурсов, её можно использовать не только для быстрого решения систем уравнений (3) на последующих временных слоях, но и для проведения новых расчётов, например, при других числах Рейнольдса или краевых условиях.

3. Результаты расчётов

Для проверки эффективности предложенного метода решения дискретизированных систем уравнений (4) и обнаружения нестационарных решений были проведены численные расчёты классической модельной задачи о течении вязкой однородной несжимаемой жидкости в квадратной каверне с неравномерно движущейся верхней крышкой и задачи о течении вязкой однородной несжимаемой жидкости в плоском диффузоре. Сначала небходимо было убедиться в работоспособности и устойчивости численного алгоритма на сколь угодно больших интервалах времени, поэтому первая серия расчётов связана с исследованием внутренней задачи со стационарными и периодическими краевыми условиями при умеренных числах Рейнольдса. Область О в этом случае представляет собой квадратную каверну с длиной стороны Ь =1. Верхняя стенка (крышка) каверны движется вдоль оси X с заданной скоростью и (¿), остальные стенки неподвижны. Вся граница дО является непроницаемой. На рис. 1, а приведены линии тока при и(¿) = 1, И,е = 400, к = 0.02, = 0.01 в различные моменты времени, демон-

Рис. 1. Линии тока течения в квадратной каверне при Ив = 400, Н = 0.02, = 0.01: а —

и(*) = 1, б — и(*) = вт(*)

стрирующие фазы начала, развития и установления течения. Полученные результаты свидетельствуют об установлении течения с ростом значений времени и согласуются с данными других исследователей [10].

На рис. 1, б представлены линии тока течения при и(¿) = вт(£) И,е = 400, к = 0.02, = 0.01 в моменты времени, отвечающие периоду движения крышки каверны. Расчёты показали, что течение на достаточно больших промежутках времени (для 20 временных периодов), как и верхняя стенка, сохраняет периодический характер. Этот факт количественно подтверждается графиком распределения по времени вихревой характеристики в произвольно фиксированной точке области. Например, на рис. 2 приведена зависимость функции вихря от времени при X = 0.2, У = 0.8.

0.01 3.01 9.01 15.01 21.01 27 01 33.01 39.01 45.01 48.01 0 -\

Рис. 2. График зависимости функции вихря от времени в задаче о течении в квадратной каверне при X = 0.2, У = 0.8, и(*) = вт(г), Ив = 400, Н = 0.02, Д* = 0.01

Рис. 3. Линии тока течения в квадратной каверне при и(*) = 1, Ив = 6400, Н = 0.02, Д* = 0.01

Вторая серия расчётов связана с поиском нестационарного решения в задаче о внутреннем течении. В этом случае расчёты проводятся при больших числах Рейнольдса, верхняя стенка каверны, как и в задаче на установление, довольно быстро достигает стационарного положения. На рис. 3 приведены линии тока в различные моменты времени, связанные с фазами начала, развития и возникновения неустановившегося режима течения. Расчёты проведены при и(Ь) = 1, И,е = 6400, к = 0.02, ДЬ = 0.01. Полученные результаты демонстрируют отсутствие стационарного режима течения в каверне при стационарном условии движения верхней стенки в случае, когда значения числа Рейнольдса находятся вблизи критического И,е = 6400. На самом деле (хотя это явно не следует из качественной картины течения) характер движения жидкости в данном случае является периодическим, на что также указывают некоторые исследователи [11, 12]. Это подтверждает график распределения по времени вихревой характеристики в произвольно фиксированной точке области (рис. 4).

Третья серия расчётов связана с исследованием характера движения жидкости в задаче протекания при умеренных и больших числах Рейнольдса. Расчётная область в этом случае представляет собой плоский диффузор длиной Ь = 6 и высотой Н = 1. Скорость входного потока имеет симметричный параболический профиль и задается своим максимальным значением с помощью функции и (Ь). На выходе ставятся мягкие условия свободного вытекания. Верхняя и нижняя стенки неподвижны и непроницаемы. На рис. 5, а приведены линии тока при и(Ь) = 1, И,е = 400, к = 0.02, ДЬ = 0.01 в различные моменты времени, демонстрирующие фазы начала, развития и установления течения.

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

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

0 3 6 9 12 15 18 21 24 27 30 33 36 39 42 45 48

■35

Рис. 4. График функции вихря при X = 0.2, У = 0.8 в задаче о течении в квадратной каверне при и (г) = 1, Ив = 6400, Н = 0.02, Дг = 0.01

а

г = 21.02

б

г = 21.02

Рис. 5. Линии тока течения в плоском диффузоре при и (г) = Ив = 400, б - Ив = 1200

1, Н = 0.02, Дг = 0.01: а -

0.01 5.01 15.01 25.01 35.01 45.01 55.01 65.01 75.01 20

15

Рис. 6. График функции вихря при X = 5.5, У = 0.5 в задаче о течении в плоском диффузоре

при Ив = 1200, и (г) = 1, н = 0.02, дг = 0.01

Рис. 7. Линии тока стационарного решения задачи о течении в канале с каверной при Ив = 400,

Н = 0.02, дг = 0.01

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

Из проведённых расчётов следует, что нестационарное решение можно получить не в любой задаче о течении вязкой однородной несжимаемой жидкости. Существенное влияние на его возникновение или отсутствие оказывает геометрия области решения. Например, на рис. 7 представлено стационарное решение известной задчи о течении жидкости в канале с каверной при И,е = 400, к = 0.02, Д£ = 0.01. При этом увеличение числа Рейнольдса не привело к появлению нестационарного решения: течение либо переходит в стационарный режим, либо (по достижении критического значения числа Рейнольдса И,е = 6000) приобретает турбулентный характер.

Заключение

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

переходит в стационарный режим и, как правило, имеет периодический характер. При этом в ряде случаев может наблюдаться отсутствие симметрии движения жидкости при симметричных геометрии области и краевых условиях. При Re < Re* происходит классическое установление течения, при Re > Re* течение теряет свойство ламинарности и становится турбулентным.

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

Список литературы

[1] Юдович В.И. Одиннадцать великих проблем математической гидродинамики // Вестник молодых ученых: Прикладная математика и механика. СПб, 2003. С. 186-192.

[2] Владимирова Н.Н., Кузнецов Б.Г., Яненко Н.Н. Численный расчёт симметричного обтекания пластинки плоским потоком вязкой несжимаемой жидкости // Некоторые вопросы вычисл. и прикл. математики. Новосибирск: Наука, 1966. С. 186-192.

[3] Яненко Н.Н. Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука, 1967.

[4] Гольдштик М.А., Штерн В.Н., Яворский Н.И. Вязкие течения с парадоксальными свойствами. Новосибирск: Наука, 1989.

[5] Кочевский А.Н. Расчёт внутренних течений жидкости в каналах с помощью программного продукта FlowVision.

http://www.tesis.com.ru/infocenter/ downloads/flowvision/fv_sgu_es04.pdf.

[6] Белоцерковский О.М. Численное моделирование в механике сплошных сред. М.: Наука, 1984. 519 с.

[7] Белолипецкий В.М., Костюк В.Ю., Шокин Ю.И. Математическое моделирование течений стратифицированной жидкости. Новосибирск: Наука, 1991. 170 с.

[8] Захаров Ю.Н. Многошаговые схемы с вариационной оптимизацией итерационных параметров. Новосибирск, 1980 (Препр. ИТПМ СО АН СССР).

[9] Захаров Ю.Н. Градиентные итерационные методы решения задач гидродинамики. Новосибирск: Наука, 2004. 238 с.

[10] Атабаев С.Ч., Брайловская B.A., Коган B.P. и др. Течение вязкой жидкости в плоской каверне // Процессы переноса в вынужденных и свободноконвективных течениях. Новосибирск: Ин-т теплофизики СО АН СССР, 1987. С. 168-176.

[11] Пухначёв В.В. Симметрии в уравнениях Навье —Стокса // Успехи механики. 2006. № 6. С. 3—76.

[12] Волков П.К. О природе движения жидкости // Вестник Югорского гос. ун-та. 2011. Вып. 2(21). С. 8-28.

Поступила в 'редакцию 14 ноября 2012 г.