Численное решение нестационарных уравнений Навье-Стокса Текст научной статьи по специальности «Математика»

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

Том 13, Специальный выпуск 4, 2008

Численное решение нестационарных уравнений

Навье—Стокса

К. С. Иванов Кемеровский государственный университет, Россия е-таП^орвр1п83@та11. ги

In this paper we propose a method for solution of the non-stationary Navier—Stokes equations describing a viscous incompressible flow. Our method relies on the minimal residuals method with a multiparametric optimization based on a componentwise minimization for the residual norm of the approximate solution.

Введение

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

например, обзор в [2]), Менее популярны методы решения системы уравнений Навье—

ф

подхода является отсутствие каких-либо существенных проблем постановки краевых

ф

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

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

Система уравнений Навье—Стокса, записанная в постановке ф — и, имеет вид

дш д(шо) д(ьи) д дЬ дх ду

Дф = и, (2)

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.

начальные условия

ф\г=о = Ф(х,у),х,у е й; (3)

краевые условия

Ф1дс = Фг(х,у,г),г е [0;т], (4)

^\дс = ф2(х,у,1),1е[0]Т]. (5)

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

Ф

дг дх ду

ф |4=о = ф(x,y), х,у е С; (7)

Ф 1дс = Фг(х,у,1), г е [0;т], (8)

^\дс = ф2(х,у,г), ¿€[0;Т]. (9)

В системах уравнений (1)-(5) и (6)-(9) V > 0 — коэффициент вязкости; ф(х,у), фг(х, у, г), ф2(х, у, г) — заданные функции своих аргументов; С — выпуклая одноевязная область решения; дС — гладкая граница области С, Будем считать, что задачи (1)-(5) и (6)-(9) имеют единственное решение [4],

2. Численное моделирование

Введем в области С неравномерную по г, х, у, согласованную с границей дС сетку С^, Аппроксимируя задачи (1)-(5) и (6)-(9)

па сетке С^, некоторыми разностными схемами,

получим разностные задачи:

начальные условия

краевые условия

п+1/2 Ми - Ш,

;п

Н , тп+1 п+1/2 Тп+1 п _ гп+1/2. /1Пч

т/2 ,п+1 — , ,п+1/2

;^ ^ ^нх шн т шн — ¿ь, >

+ ЬП++1^П+1/2 + Ь1+1МП+1 = /П+1; (11)

и

ЛнФП+1 = МП+1; (12)

мП+1|дс = ИнфП; (13)

ФП+1|дС = ФП+1 (14)

- АпФ1 + ь1+ж+1 = (15)

т

К1гФП+1 = 9П+1. (16)

В задачах (10)—(14) и (15), (16) Л}г есть некоторая аппроксимация оператора Лапласа; ЬП+\ ¿П+1 и ЬП+1 — некоторые аппроксимации конвективных слагаемых в уравнениях движения (1) и (6); М^ — какой-либо вариант условий Тома [2], а уравнение (16)

есть разностный аналог краевых условий (8) и (9), Заметим, что операторы Ь^1, Ь^^1 и могут быть как линейными, так и нелинейными в зависимости от способа аппроксимации конвективных слагаемых в (1) и (6), В случае, если для аппроксимации этих величин используются значения на верхнем временном слое, то указанные операторы являются нелинейными, в случае же использования для аппроксимации конвективных слагаемых значений с нижнего временного слоя данные операторы станут линейными и мы получим линеаризованные разностные схемы.

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

Независимо от способа аппроксимации конвективных слагаемых разностные задачи (10)—(14) и (15), (16) для каждого дискретного момента времени можно записать как систему алгебраических уравнений вида

Л(п,п) = /, (17)

где п, / — векторы размерности т (число узлов сетки); Л(п, у) = Л1(п, у) + Л2у Л2 —

Л1

Л1(еш(1) + 6п(2),П1У(1) + П2У(2)) = е1П1Л1(п(1),у(1)) + +е1П2Л1(п(1),у(2)) + С2П1Л1(п(2),у(1)) + ^2П2Л1(п(2),у(2)), (18)

п(г), у(г) — произвольные векторы раз мерности т; £г,пг _ произвольные постоянные, г = 1, 2.

Отметим, что в случае линеаризованных разностных схем на каждом шаге по времени мы имеем систему линейных алгебраических уравнений (Л1 = 0), матрица которой, однако, зависит от временного слоя. При этом достаточно сложно установить некоторые свойства матрицы получившейся системы линейных алгебраических уравнений (например, неособенность и знакоопределенность), которые позволили бы применять богатый арсенал методов для ее решения [5]. В случае же нелинейной (билинейной) системы

Л1 = 0

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

Независимо от того, является ли система (17) линейной или нелинейной, для ее решения построим единообразный итерационный процесс [6]

пп+1/2 = пп - гга+1[Л(пга,пга) - /]; (19)

пп+1 = пп+1/2 + ага+1Жга, п =1, 2 ..., (20)

где хп — некоторый вектор размерности т; п0 — произвольное начальное приближение из области определения оператора А; тп+1,ап+1 — итерационные параметры.

Пусть ап+1 — квадратная матрица с т ненулевыми элементами аП+\ г,^ = 1...т, к — произвольное целое число от 1 до т. Перепишем (20) в виде

т

пп+1 = Урп_+11/2 + £ 1е'г, Р, 3 = 1, 2 ... т, (21)

г=р

где У;+11/2 = ип+1/2 + аЦ V1 + .. ненулевой Аг-й компонентой. Введем обозначение:

+ е^1, уП+1/2

и

га+1/^ _ вект0р с ОДНОЙ

Г

п+1/2

'(О

А

У—11/2 + «й^, уП-+11/2 + ^

- /, г = 1, 2 ...т.

(22)

Очевидно, что Гп+1 = Гт+1/2 = А(ип+1, ип+1) — / и гп+1/2 = гп = А(ип, ип) — / — невязки схемы (20), Переписывая (21) относительно нормы невязки и выбирая аП+1 из условия

минимума ||ггп+1/2||2 [6], можно получить

„п+1/2 II <

п+1/2 г—1

, г = 1, 2 ...т, п =1, 2 ...

(23)

Неравенство (23) означает, что на каждом итерационном шаге норма вектора невязки не возрастает. Необходимо отметить, что в случае линейной системы уравнений (А1 = 0) можно показать [7], что ||гп|| ^ 0, при п ^ то, т.е. итерационный процесс (19), (20) сходится при любом начальном приближении.

Приведенный алгоритм означает, что элементы матрицы ап+1 выбираются последовательно, исходя из условия минимума соответствующей невязки, В ряде случаев удается использовать не последовательную, а многопараметрическую оптимизацию [7]. В этом случае

иГ1/2 = и—11/2 + <4 + ... + <4, (24)

гГ+1/2 = ЖиГ1/2,иГ1/2) — /, г = 1, 2 ..., (25)

и итерационные параметры «П+1... аПг+1, составляющие группу, выбираются из условия

глобального минимума нормы вектора г™+1/2. Для их определения па каждом итерационном шаге необходимо решить систему алгебраических уравнений (линейных или нелинейных в зависимости от исходного оператора), размерность которой равна числу итерационных параметров в группе. Например, если оператор А является линейным, то система линейных алгебраических уравнений для определения итерационных пара-

(п+1) (п+1) " „

метров а ... а имеет следующий вид:

/ЧАяП ^) (А^П ^) ... (А^ (АгП (А*, Azn2) ... (А*

Aznг) \ ^Пг)

X

V (-, Aznl) (Azn2, Aznг) ... (Aznг, Aznг) /

п+1

П+1 к2

V /

(

\

((А (А^

zn1, гг+1/2) n

п к2 ,

г

п+1/2

v (aznг ,гП+1/2)

Отметим, что в большинстве случаев матрица исходного оператора имеет блочно-диагональную структуру. Например, при решении разностной задачи Дирихле для уравнения Пуассона матрица оператора имеет вид

\

А

( ) а1,1 ( ) а1,2 0 0 0 ( ) а1,в 0 0

( ) (п) а2,2 ( ) 2,3 0 0 0 ( ) 2,7 0

0 (п) а3,2 ( ) 3,3 ( ) 3,4 0 0 0 1 а.

0 0 ( ) а4,3 ( ) 4,4 ( ) 4,5 0 0 0

0 0 0 ( ) 5,4 ( ) <5 0 0 0

( ) 0 0 0 0 ( ) аб,в ( ) 6,7 0

0 ( ) 7,2 0 0 0 ( ) а7,6 ( ) 7,7 1 а

0 0 ( ) 8,3 0 0 0 ( ) 8,7 1 а.

(п) 3,8

(п) 7,8 (п) '8,8

V

(27)

/

поэтому, если итерационные параметры сгруппировать следующим образом:

{а

п+1 п+1 п+1 п+1

а

а

12

а

15

...}, {а

п+1 ап+1 а

а?з+1, ап+1 а16 , . .. }

а?8+1, ап+1 а21 , . .. }

ап+1 а20 , ап+1 а23 , . .. }

(28)

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

Особо отметим случай, когда оптимизация проводится по всем итерационным пара-

. а^+1 одновременно. Для определения данных параметров необходимо

метрам аП+1

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

п+1

... а

п+1

имеет вид

^^п^п) (А7п,А7п) (А7п,А7п) (А7п,А7п)

V (АС,А^п) (АС,А^п)

(Агп ,А^) \ (Агп ,А^)

(Ас, /

а1 +1

X

+1

а

а +1 \ т

( (Azn ,гп+1/2) \ (Azn ,7

/

(29)

V ^т,гп+1/2) )

В общем случае решение системы представляет не меньше трудностей, чем решение исходной системы уравнений. Однако, если выбрать векторы ... ¿т так, чтобы

(Аг-,Аг-) = 0, г = то матрица системы уравнений для определения ап+1... а^^1 будет иметь диагональный вид и определение точного решения данной системы уже не представляет труда. Более того, можно показать [7], что в линейном случае при таком выборе итерационных параметров схема (19), (20) будет сходиться к точному решению за одну итерацию.

Алгоритм групповой оптимизации по всем параметрам одновременно оказывается весьма полезным при проведении серийных расчетов нестационарных задач в линейном случае. Например, если использовать разностную схему (10)—(14) для решения нестационарной системы уравнений Навье—Стокса, то можно заметить, что на каждом дискретном временном шаге возникает необходимость в решении разностного уравнения Пуассона, где разностный оператор не зависит от временного шага. Поэтому для решения данной задачи можно использовать алгоритм с полной групповой оптимизацией, затратив на первом дискретном временном шаге некоторое время на построение

1

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

Если система (17) нелинейная, то в случае плохой сходимости схемы (19), (20) для нее аналогично линейному случаю [8] можно построить процедуру ускорения сходимости, суть которой заключается в комбинации приближений схемы (19), (20) на n-м и (n + 2)-м итерационных шагах:

xn+2 = (1 + un)un+2 - unun, (30)

где un+2, un — приближения схемы (19), (20), a un выбирается из условия [7]

min 11 г™+2 || = тт\\А(хп+2,хп+2) - /||. (31)

Заключение

Для оценки предложенного метода решения задач (1)-(5) и (6)-(9) проведены численные расчеты классической модельной задачи о течении вязкой однородной несжимаемой жидкости в квадратной каверне с неравномерно движущейся верхней крышкой и задачи об обтекании вязкой однородной несжимаемой жидкостью обратного уступа. Расчеты проводились при различных значениях коэффициента вязкости и при различных функциях скорости движения крышки и входного потока. Проведенные расчеты показали эффективность предложенных алгоритмов. При достаточно больших значениях шага по времени получены устойчивые результаты. Также обнаружены нестационарные решения данных задач при стационарных краевых условиях.

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

[1] Лойцянский Л.Г. Механика жидкости и газа. М.: Наука, 1987. 840 с.

[2] Роуч П. Вычислительная гидродинамика: Пер с англ. М.: Мир, 1980. 616 с.

[3] Beam R.M. Newton's methods for the Xavicr Stokes equations // Comput. Mech. 1988. Vol. 2. P. 51.II.1-51.II.4.

[4] Ладыженская А.О. Математические вопросы динамики вязкой несжимаемой жидкости. М.: Наука, 1970. 340 с.

[5] Самарский A.A., Николаев Е.С. Методы решения сеточных уравнений: Учебник для вузов М.: Наука, 1978. 592 с.

[6] Захаров Ю.Н., Егорова Е.Ф., Толстых М.А., Шокин Ю.И. Метод минимальных невязок решения одного класса нелинейных уравнений. Красноярск, 1991 (Препр. № 9).

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

[8] Николаев Е.С. Нелинейное ускорение двухслойных итерационных методов вариационного типа // Журн. вычисл. математики и мат. физики. 1976. № 6. С. 1387.

Поступила в редакцию 20 февраля 2008 г.