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

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

Том 8, № 5, 2003

ИТЕРАЦИОННЫЕ СХЕМЫ РЕШЕНИЯ СИСТЕМЫ УРАВНЕНИЙ НАВЬЕ - СТОКСА В ПЕРЕМЕННЫХ ФУНКЦИЯ ТОКА, ВИХРЬ

М.Ю. Балаганский, Ю.Н. Захаров

Кемеровский государственный университет, Россия e-mail: [email protected], [email protected]

The steady Navier — Stokes equations in vorticity-stream function formulation are solved numerically using an explicit iterative scheme with the multicomponent optimization. Nonreflecting boundary conditions at the inflow and outflow boundaries are based on a certain approximation of the Navier — Stokes equations. Three test cases are presented.

Введение

Решение стационарной задачи движения вязкой несжимаемой жидкости в переменных завихренности ш, линии тока ф разностными методами имеет длинную историю [1, 2]. Из всех проблем, рассматриваемых в литературе, можно выделить две. Первая — это правильная постановка краевых условий для ш на твердых стенках и краевых условий для ф, ш на “свободной границе” (например, на выходе из канала). Эта проблема достаточно полно отражена в литературе (см. [3] и приведенную там библиографию). Вторая — построение экономичных итерационных методов решения разностных схем, аппроксимирующих систему уравнений Навье — Стокса. В основном используются неявные разностные схемы, аппроксимирующие нестационарную систему уравнений Навье — Стокса с последующим стационированием по времени. Такой подход требует выбирать временной дискретный параметр т не из условия быстрой сходимости нестационарной схемы, а из условия аппроксимации производной по времени и устойчивости разностной схемы. А это обстоятельство накладывает ограничение на скорость сходимости. При этом использование в задачах протекания [4] неявных схем в дробных шагах требует достаточно простых краевых условий на “свободной границе”. Гораздо реже встречаются работы, в которых разностная задача рассматривается как система нелинейных уравнений (см. например [5, 6], так как выбор начального приближения и итерационных параметров, обеспечивающих сходимость итерационных схем, является достаточно трудной задачей.

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

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

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

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

Рассмотрим двухмерное стационарное течение однородной вязкой несжимаемой жидкости. В безразмерном виде стационарная система уравнений Навье — Стокса в декартовой системе координат в переменных вихрь ш, функция тока ф имеет вид [7]

/ д 2ш д2ш \ дфдш дфдш ^ (1)

V \ дх2 + ду2) ду дх + дх ду ’

д 2ф + д2ф = (

дх2 + ду2 Ш, ( )

где (х, у) £ П, П — прямоугольная область с границей Г; V > 0 — коэффициент кинематической вязкости. Для системы уравнений (1), (2) необходимо задать следующие краевые условия:

ф = <Pl(x,y), (х,у) £ Г, (3)

дп = ^2(х,у), (х,у) £ г, (4)

д п

где п — нормаль к Г; функции ^1(х,у), ^2(х,у) определяются заданием на Г вектора скорости w = (и, у) (<^1 — с точностью до аддитивной постоянной); функция тока ф связана с элементами вектора скорости (и, у) равенствами

дф дф

и = т-, у = -т—. ду дх

2. Разностная схема

Введем в области П неравномерную согласованную с границей Г прямоугольную сетку

П = (х* = х— + к7н, у^ = у^-1 + Нщ, г £ I, ] £ 3} с сеточной границей Г) 1, 3 —

множества индексов внутренних точек, С — множество индексов г, ], отвечающих Г). На П аппроксимируем систему уравнений (1), (2) разностной схемой

VД)Шу - ф)шу + Ь)2) ф)шу = 0, (5)

Д)фу = шу, г £ I, ] £ 3. (6)

^ т(1)п \ Т (2) { / \ дф д дф д

Здесь Ь) (фу), Ь) (фу) — некоторые аппроксимации операторов — , соответ-

ду дх дх ду

ственно; Д) — разностный оператор, аппроксимирующий оператор Лапласа со вторым порядком [8]. Краевое условие (3) заменяется точным равенством

фу = ^1(хг,уу),г,3 £ С. (7)

Краевое условие для ш на твердой стенке мы задаем со вторым порядком аппроксимации по формулам, приведенным в [4].

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

Если ввести вектор и = {шу, і Є І, і Є 3, і, і Є О, фу, і Є І, І Є 3} размерности т, то разностная задача (5)-(7) запишется в виде системы билинейных уравнений

А(и, и) = ї, (8)

где А(и, и) = Аі(и, и) + А2и, Аі(и, и) — билинейный оператор, получающийся после ап-

дф дш дф дш

проксимации нелинейной части уравнения (1),------ —-—+ ———; А2 — линейный опера-

ду дх дх ду

тор, получающийся после аппроксимации линейной части задачи (1)-(4); ї — известный вектор, зависящий от функций <^(х, у), і = 1, 2.

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

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

Для решения (8) в вещественном евклидовом пространстве Нт рассмотрим явную итерационную схему

un+l/2 = un - Tn+l[A(un, un) - f],

un+1 = un+l/2 + an+lzn, n = 0, І, 2,.

(9a)

(9b)

где zn 6 Hm — некоторый вектор; u0 — произвольное начальное приближение из области определения оператора A; rn+1 = const; an+1 — итерационный параметр, который может быть как константой, так и матрицей.

Если an+i = const, то норма невязки rn+1 = A(un+1, un+1) — f имеет вид

„n+l і

rn+l/2 її2

где un+1 — решение (9b) и 0ln = 2(rn+l/2, Fln),

03n = 2(F1n, F2n) ,

Fln = A(un+l/2, zn) + Al(zn, un+l/2),

k=l

02n = (Fln, Fln) + 2(rn+l/2, F2n),

04n — (F2n, F2n) ,

F2n = Al(zn, zn).

Параметр ага+1 можно выбирать из условия минимизации ||гга+1||2. Для этого необходимо решить уравнение

nn+l |

да

n+l

5>an+l0

kn =

О.

k=1

2

4

2

Корни этого уравнения находятся по явным формулам Кардано, и из трех возможных корней выбирается тот, который глобально минимизирует ||гп+1||. Следствием такого выбора ап+1, в случае если оператор А и вектор гп таковы, что Е1п = 0, имеет место неравенство ||гп+1|| < ||гп+1/2||. Если и тп+1 выбирается аналогичным образом, то ||гп+1|| < ||гп||, п = 0,1,...

Монотонное убывание нормы невязки в наших расчетах обеспечивал следующий вариант схемы (9).

Пусть ап+1 = (Ои^} — диагональная матрица. Тогда для (9Ь) справедливо равенство

где yn0+1 = un+l/2, z, = (О,... , О, І, О,... , О) и аП+. выбирается последовательно, минимизируя Цг^+.Ц2 = IA(y‘i+1,yni+1) — f||2. Здесь аП+. выбирается по тем же формулам, что и в случае an+l = const. Это означает, что переход от компоненты вектора un+l/2 к компоненте вектора u‘+. происходит оптимальным образом (подробнее см. [9]).

Развивая далее подход вычисления каждой компоненты вектора un+1 со своим оптимальным итерационным параметром, можно находить диагональные элементы матрицы an+l не последовательно, а из условия глобального минимума ||гп+.||.

Имеем

u‘+. = u‘+l/2 + an+lzn + ■ ■ ■ + аП”?1<.

Используя билинейность оператора A, получим гп+. = A(un+1, u‘+.) — f =

= Al(un+l/2, un+l/2) + an+lAl(un+l/2, z‘) + ■ ■ ■ + аП+lAl(un+l/2, z;)+

+ a‘+lAl(zn, u‘+l/2) + (аП^) Al(zn, z,) + ■ ■ ■ + anl))lan+)lAl(zn, zm)+

+ ^їи.ос u‘+l/2) + аП1+іаП+)1Al(zm, zn) + + (an+l) Al(z;m zm)+

+ A2u‘+l/2 + ar^-)l A2zn + ■ ■ ■ + an+)lA2zm — f.

Пусть оператор a. и векторы z, такие, что

A.(z‘, zn) = 0, Vi, j = І,...,т, (10)

тогда

гп+. = rn+l/2 + anl+lG(z‘‘, un+l/2) + ■ ■ ■ + a‘+)lG(zm, un+l/2),

где

G(zn, un+l/2) = A(un+l/2, zn) + a. (zn, un+l/2),i = І, 2,...,m.

Для краткости вместо G(z‘, un+l/2) будем записывать G(z‘). Отсюда

||rn+l||2 = I r‘+l/2 12 +

+ a;,1>1(r"+l/2, G(z?)) + ■ ■ ■ + ai’+)1(r"+l/'2, G(z;))+

+ (G(zi),G(zi)) + ■■ ■ + 2a!>1+la4+i(G(zi),G«)) + ■ ■ ■ +

+ (a,‘:+l)2 (G(z?),G(zm)).

Для минимизации уравнений

к.га+11

(І)

по ап+1, і

, т, необходимо решить систему линейных

Г.П+1 I

да

(і)

п+1

0, і = 1, т,

которая в развернутом виде выглядит следующим образом:

/(0^,0^)) (О^),О^))

(ою,ою) (ою,ою)

(ОК),О^)Л

(О(г5),О(гт))

х

У(°ю,°к)) (ою,ок)) •••

х

(а(1+Л

(2)

V

'сП+1

(т) +1/

^-(гга+1/2, О^))\ -(гп+1/2,о^))

У-(г н^^ощу

(11)

При решении конкретной задачи условие (10) может быть трудновыполнимым, поэтому

(*)

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

Пусть I — множество всех индексов (г : г = 1,т}. Разобьем I на непересекающиеся подмножества 11,12,...,4, для каждого из которых выполняется условие

А^™, zn) = 0, Уі,і Є Іі, І = 1,

к.

Тогда можно записать

уп+1

у(0

уп+1 ,у^ а(і) у(і-1) + 2-^ іЄ/і

п+1

1, 2,

,к,

где У^

и

п+1/2

Параметры а(і+1 при і Є І будем находить из условия глобального минимума

„га+11

Г(і) I

г

п+1/2 '(і-

1/2 + £ а^О^, у",-1))

іЄ/і

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

Кроме того, если для каждой из I выполняется условие

(еде,У3-у,с^,у(1_1))) = 0, г = *

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

Замечание. Если система (8) является сеточной задачей, то матрица А является ленточной и нахождение множеств 11, 12,... ,1^ осуществляется алгоритмом, рассмотренным в [10].

Если система (8) имеет единственное решение, то предложенные выше алгоритмы могут обеспечить практическую сходимость схемы (9) для произвольных начальных данных. Если же (8) имеет не одно решение, то схемой (9) мы будем находить обобщенное решение (8). В сочетании с решением задачи на последовательности сеток предложенная итерационная схема позволяет получать решения, совпадающие с натурными и численными экспериментами других авторов [11-13].

1

2

4. Численные расчеты

Для проверки эффективности предлагаемых алгоритмов была проведена серия численных экспериментов по решению системы (1)-(2), описывающей внутреннее течение и течение в бесконечных каналах:

Задача 1. Течение в квадратной каверне.

0-тг

дп

Гі

дф , д п

Г 2

Здесь и далее п — нормаль к границе.

0

Г

Задача 2. Течение в канале с каверной.

Гз

п Г4 1 1

г2 Г2 Г2

ф дф д п

Гі

= Р1(У), ф

Г2

= 0, ф

Гз

Н - к,

Г2иГз

= 0, ф —> Р1(у) при х —> то.

Задача 3. Течение в канале с уступом.

т4 ф = р1(у), ф = 0, ф

Гі Г2

Н - к,

Гз

дп

Г2иГз

0, ф —> Р2(у) при х —> то.

Здесь Н — ширина канала; к — высота уступа; функции Р1(у) и Р2(у) задают течение Пуазейля в канале шириной Н — к и Н соответственно.

В области П с размерами X х У строится неравномерная, прямоугольная, с числом точек N. х N, сгущающаяся вблизи границ сетка. На этой сетке аппроксимируем систему

(1), (2) со вторым порядком аппроксимации [8]. На твердых стенках из условия д^/дп = ^2(ж,у) задаются значения для завихренности ш со вторым порядком точности [4].

Получающуюся разностную задачу мы рассматривали как систему нелинейных уравнений (8), которую решали итерационной схемой (9) в различных модификациях. Итерации продолжались до выполнения условия ||гп|| < е, где гп = А(ип, ип) — ¥ — невязка.

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

Задача 1 является тестовой, и ее численное и экспериментальное решение хорошо известно. Поэтому решение именно этой задачи рассматриваемыми выше итерационными алгоритмами позволило проверить их работоспособность. Как видно из рис. 1, полученное решение имеет один центральный и два вторичных вихря и оно хорошо согласуется, например, с результатами работ [11, 12].

Хотя существует несколько способов переноса краевых условий с бесконечности на границу конечной области (см., например, [4]), все они не позволяют поставить правую границу Г4 близко к препятствию, где поток жидкости резко меняет направление движения.

Мы, как и в работах [14, 15], предлагаем на границе конечной области Г4 (см. рисунки к задачам 2,3) в качестве краевого условия задавать некоторую аппроксимацию самой системы уравнений (1), (2) внутрь области решения. Тогда разностная задача записывается в виде (8), матрица которой как минимум шестидиагональна. В этом случае использование схем расщепления нецелесообразно и мы решали (8) с помощью итерационной схемы (9).

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

Рис. 1. Течение в каверне: V = 0.002, X х Рис. 2. Течение над каверной: V = 0.002, X х

У = 1 х 1, Хх х Ху = 51 х 51, е = 2 ■ 10-7.

г0У = 29.90, У = 2.3 х 1, Хх х Ху = 93 х 41,

1 ■ 10-4.

= 7.60, е =

0

г

Рис. 3. Течение над каверной: V = 0.002, X х Рис. 4. Течение над каверной: V = 0.002, X х

У = 3 х1, Хх хХу = 121 х 41, ||г01| = 7.85, е = У = 2.3х1,Хх хХу = 185 х81, ||г0|| = 10.05, е =

1 ■ 10-4. 1.4 ■ 10-3.

решений в меньшей области.

Как показали расчеты задачи 2, внешний поток, обтекающий “каверну”, практически не меняется, и в этом случае увеличение области решения вниз по потоку совершенно не влияет на течение в “каверне” (сравним рис. 2, 3). Увеличение числа узлов в меньшей области только уточняет решение (рис. 4).

Трудности решения задачи 3 связаны с тем, что при уменьшении коэффициента вязкости V увеличивается область, занятая вихрем за уступом. Например, для V = 0.05 вихрь

V = 0.05, X х У = 1 х 1, Хх х Ху = 81 х V = 0.01, X х У = 1.5 х 1, Хх х Ху = 31 х 21,

81, ||г0|| = 40.08, е = 3 ■ 10-4 ||г01| = 9.88, е = 1 ■ 10-4

V = 0.0025, X х У = 1.5 х 1, Хх х Ху = 61 х V = 0.0025, X х У = 2 х 1, Хх х Ху = 41 х

41, ||г0|| = 3.13, е = 3 ■ 10-4 41, ||г0|| = 3.11, е = 2 ■ 10-3

V = 0.0025, X х У = 9 х 1, N х N = 361 х 41, ||г0|| = 3.82, е = 4 ■ 10-3

Рис. 5. Течение над уступом.

заполняет область, по длине равную 0.6 высоты уступа, при этом точка отрыва (см. также [16]) находится ниже угловой точки уступа (рис. 5, а). Для v = 0.01 вихрь заполняет уже две высоты уступа (рис. 5, б). Более того, из численных и натурных экспериментов [12, 13] известно, что на верхней стенке канала имеется вихрь, начинающийся за основным вихрем, что делает течение вниз по потоку нетривиальным. В связи с этим обычно используемые краевые условия [4] необходимо ставить достаточно далеко от уступа.

В нашем случае удается получить хорошее решение даже в области, лишь частично захватывающей вихрь. При этом увеличение области решения более чем в три раза не приводит к изменению решения в меньшей области более чем на 1 % и позволяет получить картину течения в большой области (рис. 5, в-d), хорошо согласующуюся с результатами численных и натурных экспериментов, приведенных в работах [12, 13].

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

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

[1] Thom A. An investigation of fluid flow in two-dimensions // Acr. Res. C.R. A M. 1928. N 1194.

[2] WOODS L. Note on the numerical solution of a fourth order differential equation // Aero Quart. 1954. N 5.

[3] Булеев Н.И., ТАрунин Е.Л. Исследование скорости сходимости схемы w, ф при различной структуре условия для вихря у твердой стенки // Численные методы механики сплошной среды: Сб. науч. тр. / АН СССР. Сиб. отд-ние. ИТПМ. 1984. Т. 15, № 6. С. 28-40.

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

[5] Бабский В. Г., Скловский Ю.Б. Применение метода Ньютона — Канторовича к расчету течения вязкой жидкости между вращающимися эксцентричными цилиндрами // Тр. IV Всесоюз. сем. по численным методам механики вязкой жидкости. Новосибирск: ВЦ СО АН СССР, 1973. С. 63-67.

[6] Venkatakrishnan V. Newton solution of inviscid and viscous problems // AIAA J. 1989. Vol. 27, N 7. P. 885-891.

[7] Пухначев В. В. Лекции по динамике вязкой несжимаемой жидкости. Новосибирск: Изд-во НГУ, 1969. Ч. 1. 198 с.

[8] Самарский А. А. Введение в теорию разностных схем. М.: Наука, 1971. 552 с.

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

[10] БАлАгАнский М.Ю., Захаров Ю.Н., Ханевт В. Использование итерационных методов в решении нестационарных задач движения стратифицированной жидкости // Вычисл. технологии. 2002. Т. 7, № 5. С. 3-10.

[11] Ghia U., Ghia K.N., Shin C.T. High-re solutions for incompressible flow using the Navier — Stokes equations and a multigrid method // J. Comput. Phys. 1982. Vol. 48. P. 387.

[12] Rogers S.E., Kwak D. An upwind differencing scheme for the incompressible Navier — Stokes equations // Appl. Numer. Math. 1991. Vol. 8. P. 43-64.

[13] Armaly B.F., Durst F., Periera J., Schonung B. Experimental and theoretical investigation of backward-facing step flow // J. Fluid Mech. 1983. Vol. 127. P. 473.

[14] ЗАхАров Ю. Н. Об одном методе решения уравнений с краевыми условиями на бесконечности // Вычисл. технологии. 1993. Т. 2, № 7. С. 55-68.

[15] Захаров Ю.Н., Кривушин С. А. Метод минимальных невязок решения системы уравнений Навье — Стокса // Вест. Кем. гос. ун-та. Математика. Вып. 4. Кемерово, 2000. С. 108-113.

[16] Тейлор Т.д., Нфедо Э. Расчет течения вязкой жидкости в канале при помощи метода расщепления // Численные методы в механике жидкостей. М.: Мир, 1973. С. 218-229.

Поступила в редакцию 7 апреля 2002 г., в переработанном виде — 28 июля 2003 г.