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

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

Том 11, № 2, 2006

ОБ ОДНОМ АЛГОРИТМЕ РЕШЕНИЯ УРАВНЕНИЙ НАВЬЕ - СТОКСА ВЯЗКОЙ НЕСЖИМАЕМОЙ ЖИДКОСТИ*

В. М. КОВЕНЯ

Институт вычислительных технологий СО РАН, Новосибирск, Россия

e-mail: [email protected]

An economical difference scheme is proposed for solving the Navier — Stokes equations for a viscous incompressible fluid in physical variables. The scheme is based on splitting the equations on physical processes. The scheme is implemented using scalar sweeps. The algorithm is tested on the problem of the convective flow in a cavity and on the problem of a flow in the closed domain with one heated wall.

Введение

Уравнения Навье — Стокса вязкой несжимаемой жидкости не являются системой типа Коши — Ковалевской, и поэтому для их численного решения не удается напрямую применять численные методы, разработанные для решения гиперболических или параболических уравнений [1-3]. Вторая особенность этих уравнений заключается в отсутствии в постановке задач краевых условий для давления. Поэтому традиционно задачи о течениях вязкой несжимаемой жидкости решаются в формулировке "функция тока — вихрь". В этих переменных система уравнений сводится к параболической системе уравнений для компонент вихря и эллиптической системе для векторного потенциала (в двумерном случае — скалярных уравнений для вихря и функции тока). При необходимости давление восстанавливается из уравнения Пуассона, полученного из уравнений движения. Отметим, что для его решения задаются граничные условия второго рода (условия Неймана), что вызывает определенные трудности при его реализации [2-9]. При решении многих, особенно пространственных, задач использование физических переменных может оказаться предпочтительным как из-за уменьшения числа решаемых уравнений (в трехмерном случае), так и вследствие практической необходимости вычисления поля давления. Для численного решения уравнений в физических переменных наиболее часто применяются разностные методы и методы конечных элементов [1-15]. Большая часть из конечно-разностных методов использует предположение, что на промежуточном слое течение является соленоидальным, т.е. полагается, что градиент давления равен нулю, а на полном шаге решается уравнение для давления, т.е. фактически используется идеология расщепления. Отметим также метод искусственной сжимаемости, состоящий в замене уравнений

* Работа выполнена при частичной финансовой поддержке Российского фонда фундаментальных исследований (грант № 05-01-00146) и интеграционного проекта СО РАН № 116.

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

Навье — Стокса несжимаемой жидкости моделью слабо сжимаемой жидкости путем введения в уравнение неразрывности дополнительного члена — производной по времени от давления с малым параметром [5, 9]. Заметим, что в рамках этой модели необходимо обеспечивать сходимость решения модифицированного уравнения неразрывности к исходному уравнению при стремлении малого параметра к нулю.

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

1. Исходные уравнения и постановки задач

В качестве исходной модели рассмотрим уравнения Навье — Стокса вязкой несжимаемой жидкости в приближении Буссинеска в безразмерном виде (см., например, [1, 2, 4, 7-9, 14, 15]):

а^у = о,

дУ + (V ■ У)У = -Ур + ду/Ив - Ragт, (1)

дТ

— + (V ■ у)Т = ДТ/(ИвРг).

Здесь У = (V ,У2) — вектор скорости; р, Т — давление и температура, отнесенные к их невозмущенным значениям Ц, р0,Т0; Ив = ЬЦ0/V — число Рейнольдса; Иа = аgЬT0/ Ц0 — число Рэлея; Рг = V/ к — число Прандтля; V, к — коэффициенты кинематической вязкости и теплопроводности, значения которых полагались постоянными; g — ускорение свободного падения, g = (0,1); а — коэффициент объемного расширения. При переходе к безразмерному виду компоненты скорости нормировались на значение и0, давление — на Цц, а температура — на Т0. Плотность полагалось постоянной, она исключена из уравнений (1). В качестве характерного размера задавалось значение Ь = 1.

Предлагаемый ниже численный метод апробировался на решении задач в области Q с границей Г:

1. Течение в каверне, верхняя стенка которой двигалась со скоростью V = 1. В этом случае уравнение для температуры не решается (Иа = 0). На боковых и нижней стенках области задавались условия прилипания У = 0.

2. Течение теплопроводного газа в замкнутой области Q. На ее границах задавались условия прилипания. На отдельных границах задавалась различная температура Т (возможно задание на части границ условий тепловой изоляции дТ/ дп = 0).

Для указанных задач имеются многочисленные расчеты [3-5, 7-12], что позволяет провести сравнения по точности предлагаемого метода и дать оценки его эффективности.

Запишем исходные уравнения (1) в консервативной форме в виде

2

м" + £ А ш,

,=1

дх,

0,

(

где f

('Р \

и>1

\Т )

и,

1

V

г^и, + Ь,1' - ■ д^1/ дх, и2и, + <у2'р - И ■ ди2/ дх, + Иа ■ Т и,Т - „ ^ ■ дТ/ дх,

М

(2)

0000 0100 0010 \ 0 0 0 1

, Ие ■ Рг ' ,

Заметим, что уравнения Навье — Стокса не являются системой типа Коши — Ковалевской, так как матрица М вырожденная [2, 9].

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

Введем в расчетной области Q равномерную сетку. Пусть т, к2 — шаги сетки по времени £ и пространству х1, х2 соответственно, Л,, Л„ — симметричные разностные производные, аппроксимирующие производные д/ дх, и д2/ дх, = 1, 2) на трехточечном шаблоне со вторым порядком, fп = ¡2 — сеточные функции, заданные в узлах сетки п, ¿1, 12. Аппроксимируем систему уравнений (2) разностной схемой с весами

^+1_х?п 2

М-+ ^ Л, (ашп+1 + ршп) = 0, в = 1 - а. (3)

т ,=1

Схема (3) нелинейна и аппроксимирует исходные уравнения (2) с порядком О(тт + к2 ), где т =2, если а = 0.5 + О (т), и т =1 — в противном случае.

Линеаризируем векторы ШП+1 относительно вектора Р+1 по формуле

Шп+1 = Шп + т (дШ,/дf)п ■ (дf/д£)п + О(т2) = Шп + ДпАГ + О(т2),

где

ДП

( 0 Ь,1 Ь,2 0 \

1

, V-1- 1 , и, 1 ^ V"-

, V- - , и2 (1 + ,

V 0 0 0 и, - Л,/(ИеРг) У

Ь1 (1 + Ь1) и, - Л, / Ие (1- , и 0

Ь2 (1 - Ь,2) и (1 + Ь,2) и, - Л, / Ие 0

/п _ / п+1 _ /

С учетом линеаризации разностная схема (3) примет вид

(М + таДп) Afn = -т ^ Л,Шп. (4)

,=1

Здесь Д = Л1Д1 + Л2Д2, А^ = Л11 + Л22 — разностный оператор Лапласа. В силу эквивалентности дивергентной и недивергентной форм записи уравнений в дифференциальной форме в разностном виде будем иметь

Д = В + О(к2), (5)

2

0

где В =

Л?

V

Л1

Л1 Е V]Л3 - Лл/Ив

3=1

0

Л2 0

2

Е V]Л] - Лл/ Ив

3=1

0 0 0

Тогда разностная схема (4) примет вид

0 0

Иа

ЕиЛ] - Лл/(ИвРг) 3=1

(М + таВ) ДГ = Л3

(6)

3=1

Она аппроксимирует исходные уравнения с тем же порядком, что и базовая схема (3). Применяя расщепление по физическим процессам [1], приближенно факторизуем стабилизирующий оператор М + таВ оператором С со вторым порядком по формуле

Здесь

М + таВ « С = (I + таВ1) (М + таВ21) = I + таВ + т2а2ВВ

0

В1

0

0 Е иЛ,- - Лл/Ив

3=1

V

0 0

0 0

0 0

2

Е и Л3 - лл/ Ив

3=1

0

0 0

Иа

2

Л3 - Лл/ИвРг

3=1

В2

0 Л1 Л2 0 Л1 0 0 0 Л2 0 0 0

\ 0 0 0 0 у

С учетом факторизации (7) разностная схема (6) перепишется в виде

(7)

С ДГ = -т ^ Л3 ^

3=1

или в виде эквивалентной ей схемы в дробных шагах

2

Г = -т е Л3 wn,

3=1

(I + таВ1) |га+1/2 = Г, (М + таВ2) |га+1 = |га+1/2,

£га+1 = р. + т

(8)

(9)

Здесь ^ = (£р,£ь£2,£т) — вектор невязки. Как и схема (3), схемы (8) или (9) аппроксимируют уравнения (2) с порядком 0(тт + к2).

2

2

Остановимся на реализации схемы (9). Значения вычисляются явно из первого уравнения схемы, причем Л3-ип = 0, т. е. уравнение неразрывности удовлетворяется тождественно. На первом дробном шаге решается система уравнений

£П+1/2 = с = 0,

I + та ( ЕиЛ3 - Лл/Ив 3=1

I + та ( Е иЛ3 - Лл/Ив 3=1

¿п+1/2 _ ¿п ?1 = ?1 ,

еп+1/2 = еп + тае п+1/2,

(10)

I + та ( ЕиЛ3 - Лл/(ИвРг) 3=1

еп+1/2 _ еп

ет = ?т •

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

(6)

п+1/2, V

е3,г> Ли 3,г быть записана в виде

п+1/2, V+1

е

V+1

3,г

е3',г. Тогда итерационная схема для него может

е3,» + а1 (е3+1,г - е3-1,») - Ь1 (€3+1,1 - + €3 —1,г) +

+«2 (е^ - о,—о - Ь2 (е^ - 23+€3,^—1) = е

п,

где « = таиг/(2кг), Ь = та/ (кг2Ив). Это уравнение может быть решено трехточечными скалярными прогонками до сходимости итераций, т. е. до выполнения условия

тах |3 - < 0(тк2)

(11)

для всех г и ^. При сходимости итераций получим значение функции £п+1/2. Наряду с итерационными методами можно рассмотреть решение уравнений (10) и по схеме приближенной факторизации вида

П [I + та (и-Л3 - Л33/Ив)]еп+1/2 = еп

3=1

.п+1/2

или эквивалентной ей схеме в дробных шагах

[I + та (и1Л1 - Лп/Ив)] п = еп, [I + та (и2Л2 - Л22/Ив)] еп+1/2 = П,

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

та (Л1еп+1+Л2еп+1) = еп+1/2, еп+1 + таЛ1 еРп+1 = еп+1/2, еп+1 + таЛ2 еРп+1 = еп+1/2,

еп+1 = еп+1/2 еТ = ?Т

2

2

решается в такой последовательности: исключая £п+1 и Сп+1 из первого уравнения системы (уравнения неразрывности), получим для ¿р^1 уравнение

т2а2А^еГ1 + г1/2 - та (л^ + Л2СП+1/2) = 5 = 0, (13)

так как А^ = Л1Л1+ Л2Л2. Его решение может быть получено подобно итерационной схеме вида (11) или по схеме приближенной факторизации

¿V + 1 _ ¿V

(I - таЛи) (I - таЛ22) ^-^ = (14)

т

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

(I - таЛп) е1 = 5,

(/ - таЛ22) е2 = е1,

ег1 = ¿V + те2

трехточечными скалярными прогонками. При сходимости итераций вычисляются значения ерп+1, после чего по схеме (12) явно вычисляются значения ¿п+1 и ¿п+1 . Наконец, значения функций fn+1 определяются явно из последнего уравнения схемы (9). На этом один временной шаг заканчивается. Отметим, что уравнение неразрывности удовлетворяется тождественно на каждом дробном и полном шаге. При нахождении стационарного решения вычисления повторяются до установления, т. е. до выполнения условия

тах/п+1- /п| < О (тк2) . (15)

Остановимся на реализации краевых условий при реализации схемы (9).

При выполнении условий прилипания краевые условия для невязок скорости полагаются равными нулю

еп1г = еп+1/2г = ег+1 = 0, г = 1,2. (16)

На части границы, где задана температура, невязки для нее полагаются равными нулю подобно (16), а для границ, где задано условие тепловой изоляции дТ/дп|=0 = 0, невязка

для температуры ¿п задается с первым порядком по формуле еп+1/2(0) = еп+1/2(1), а при пересчете температуры на полном шаге — со вторым порядком односторонними разностями внутрь расчетной области по формуле [3Т(0) - 4Т(1) + Т(2)]/(2к)|г = 0.

Как отмечалось выше, краевое условие для давления отсутствует, но его задание необходимо для реализации схемы (13) на дробном шаге. Воспользуемся следствием из уравнений Навье — Стокса (1), полагая их справедливыми и на границе Г:

др/ дп = (д 2ип/дп2) / Ие.

Аппроксимация этого условия со вторым порядком даст требуемое значение давления на границе из формулы

3р(0) - 4р(1) + р(2) -10^(1) + 8^(2) - 2^(3)

к1 к1Де

Таким образом, краевые условия задачи, как и разностная схема, задаются со вторым порядком аппроксимации.

Предложенная разностная схема (8) или (9) для уравнений (2) в линейном приближении при а > 0.5 + О(т) безусловно устойчива.

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

В первой серии расчетов исследовалось течение в каверне с движущейся крышкой со скоростью и = 1 (рис. 1) в рамках модели Навье — Стокса (2) без уравнения теплопроводности. На неподвижных границах задавались условия прилипания и = и2 = 0. В начальный момент £ = 0 в расчетной области (кроме верхней стенки) жидкость полагалась неподвижной. Расчеты проведены при различных числах Рейнольдса. Для оценки точности разностной схемы и ее эффективности расчеты были проведены на последовательности сеток с числом узлов 21 х 21, 41 х 41, 81 х 81. На рис. 1 представлено распределение линий тока для Ие = 300 (а) и Ие = 1000 (б). С возрастанием чисел Рейнольдса размер и интенсивность вторичных вихрей в углах возрастают, как и следует из теоретических

Рис. 1.

1

I/ 1-Ы= 81

к/ 2- N=41

3- N=21

а

а

-О .4 -0.3 -0.2 -0.1 О ОН 02 03 О А 0,5 0^6 01 08 0 9 1 0 0,1 0,2 03 04 03 Об 07 08 09 1

а б

Рис. 3.

оценок и расчетов других авторов. На рис. 2 приведены распределения продольной (а) и поперечной (б) скоростей в центре каверны в зависимости от числа узлов для Яе = 300. Можно говорить о достаточной точности схемы для сеток с числом узлов 41 х 41 и 81 х 81. Подтверждением вышесказанного служит рис. 3, где приведены распределения продольной (а) и поперечной (б) скоростей в различных сечениях каверны при Яе = 1000 на сетке с количеством узлов 81 х 81. Сплошной линией отмечен результат текущих расчетов, точками — расчет из работы [1]. Все дальнейшие расчеты проводились на сетках 81 х 81.

В табл. 1 и 2 приведены значения функции тока ф и вихря и, полученные в настоящей работе и в работах [11, 12] для различных чисел Рейнольдса. Отметим хорошее соответствие результатов и для больших чисел Рейнольдса (относительная погрешность не превышает 1.5% ), хотя в [12] были использованы разностные схемы 5-го порядка, а в [11] — большее число узлов.

Результаты исследований скорости сходимости схемы к стационарному решению (удовлетворению условия (14)) при различных числах узлов N для Яе = 300 приведены в табл. 3, а числа Рейнольдса — в табл. 4. Как следует из таблиц, число итераций возра-

Таблица 1. Значения функции тока ф в центре основного вихря

Бе Текущие расчеты И^егв, К--ак [12] СЫа [11]

100 -0.1006 (81 х 81) -0.1030 (81 х 81) -0.1034 (129 х 129)

400 -0.1084 (81 х 81) -0.1131 (81 х 81) -0.1139 (129 х 129)

1000 -0.1162 (81 х 81) -0.1171 (81 х 81) -0.1179 (129 х 129)

Таблица 2. Значения вихря и в центре основного вихря

Бе Текущие расчеты Б^еге, К-ак [12] СЫа [11]

100 3.213 (81 х 81) 3.104 (81 х 81) 3.166 (129 х 129)

400 2.255 (81 х 81) 2.296 (81 х 81) 2.296 (129 х 129)

1000 2.034 (81 х 81) 2.044 (81 х 81) 2.050 (129 х 129)

Таблица 3. Число итераций для различных N Таблица 4. Число итераций для различных Яе

N Число итераций Re Число итераций

21 1231 100 1111

41 1921 400 2731

81 2521 1000 3651

стает по линейному закону с ростом числа узлов и также растет при увеличении чисел Рейнольдса.

При решении уравнения Пуассона (13) до сходимости итераций (выполнения условия типа (11), как правило, требовалось одна-три итерации, кроме первых шагов расчета, когда начальное поле данных было сильно несогласованным (значение скорости во всей области расчета равно нулю, а на пластине — 1). В этом случае требовалось порядка 100-300 итераций.

Следующая серия расчетов посвящена исследованию конвективных течений с посто-

▲

Т=1-х1

Т= 1 Т= 0

-►

Т= 1-jci xi

Рис. 4.

янной вязкостью в замкнутой плоской каверне при подогреве одной из сторон (рис. 4) в приближении Обербека — Буссинеска (1) или (2).

На границах области Г заданы условия прилипания V = 0. Температура Т подобно [8] на боковых границах задавалась постоянной и равной Т(£, 0, х2) = 1, Т(£, N , х2) = 0, а для горизонтальных участков границ — по формулам Т(£, х\, 0) = 1 — х\, Т(£, х\, N) = 1 — х\. Были проведены расчеты для различных чисел Рейнольдса, Прандтля и Рэлея. На рис. 5 представлены распределения скорости, а на рис. 6 — распределения температуры для Рг = 1, Иа = 0.1 в продольном (а) и поперечном (б) сечениях при числах Ие = 100 (кривая 1), 400 (кривая 2), 1000 (кривая 3). Возрастание чисел Рейнольдса приводит к увеличе-

Рис. 8.

Рис. 9.

нию амплитуд скорости и температуры. Изменение чисел Прандтля для фиксированных чисел Рейнольдса и Рэлея также приводит к существенной перестройке картины течения (рис. 7), а изменение чисел И,а в широком диапазоне параметров слабо влияет на конвекцию (рис. 8), как и следует из расчетов других авторов [3-10]. На рис. 9 представлены линии тока при И,е = 400 (а) и И,е = 1000 (б). Выделяются обширная центральная область циркуляционного течения ("центральный вихрь") и вторичные "угловые вихри" в боковой части каверны. Отметим возрастание скорости вихрей при возрастании чисел Рейнольдса.

Выводы

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

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

Автор благодарит Г. Григорьеву за помощь в проведении расчетов.

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

[1] КовЕня В.М. Разностные методы решения многомерных задач: Курс лекций. Новосибирск: Изд-во Новосиб. гос. ун-та, 2004. 146 с.

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

[3] Bruneau Ch.-H., Jouron C. An efficient scheme for solving steady incompressible Navier — Stokes equations // J. Comput. Phys. 1990. Vol. 89, N 2.

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

[5] ПАсконов В.М., Полежаев В.И., Чудов Л.А. Численное моделирование процессов тепло- и массообмена. М.: Наука, 1984. 288 с.

[6] Самарский А.А., Вавишевич П.Н., Матус П.П. Разностные схемы с операторными множителями. Минск: Ин-т мат. моделирования; Ин-т математики АН СССР, 1988.

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

[8] Тарунин Е.Л. Вычислительный эксперимент в задачах свободной конвекции. Иркутск: Изд-во Иркут. ун-та, 1990. 228 с.

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

[10] Kayamura T., Takami H., Kuwahara K. New higher-order upwind scheme for incompressible Navier — Stokes equations // Lect. Notes Phys. 1985. Vol. 218.

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

[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] Вабишевич П.Н., Самарский А.А. Об устойчивости разностных схем для задач конвекции-диффузии // Журн. вычисл. математики и мат. физики. 1997. Т. 37, № 2. С. 188-192.

[14] Васильев О.Ф., Воеводин А.Ф., Никифоровская В.С. Численное моделирование двумерных в вертикальной плоскости ветровых течений в стратифицированных водоемах методом расщепления // Вычисл. технологии. 2005. Т. 10, № 5. С. 19-28.

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

Поступила в редакцию 29 декабря 2005 г.