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

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

Том 14, № 2, 2009

Метод факторизации для численного решения уравнений вязкой несжимаемой жидкости

А. В. Блзовкин, В.М. Ковеня Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: [email protected], [email protected]

О. М. Вавилова Новосибирский государственный университет, Россия

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

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

Введение

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

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 08-01-00264a), интеграционных проектов СО РАН № 44, № 103.

© ИВТ СО РАН, 2009.

Для численного решения уравнений Навье—Стокса в естественных переменных часто применяются разностные методы и методы конечных элементов [2-11]. Большая часть этих методов использует предположение, что на промежуточном слое течение со-леноидальное, а на полном шаге решается уравнение для давления, т.е. в алгоритме используется идеология расщепления. Альтернативой этим подходам выглядит метод искусственной сжимаемости [8], состоящий во введении в уравнение неразрывности производной по времени с малым коэффициентом. Получающаяся система уравнений является классической — параболической для уравнений движения и гиперболической для уравнения неразрывности, но для нестационарных задач необходимо обеспечить сходимость решения модифицированных уравнений с решением исходных уравнений [8].

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

Тестирование алгоритма проведено на примере решения двумерных и пространственных задач. Важно было понять работоспособность алгоритма и в криволинейных координатах, поскольку при решении реальных задач применяется, как правило, преобразование исходных координат. В качестве двумерных задач рассмотрены плоские и осесимметричные течения около уступа и течение вязкой несжимаемой жидкости в расширяющемся канале. Решение этих задач находилось в приближении уравнений Навье—Стокса. Для повышения точности расчетов вводилось преобразование координат, сгущающее сетку в области границ. Работоспособность алгоритма в трехмерной постановке проверена на задачах о тепловой конвекции в кубе (с подогревом одной из сторон). В этом случае в качестве исходной модели уравнения Навье—Стокса дополнялись уравнением теплопроводности. Разностная схема для его решения была основана на методе приближенной факторизации. Рассмотрены режимы течения, при которых влиянием конвективных членов можно было пренебречь. Этот класс задач считается одним из основных в геофизике при исследовании конвекции в недрах Земли [3, 10, 12, 13].

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

В качестве исходной модели будем рассматривать уравнения Навье—Стокса вязкой несжимаемой жидкости в приближении Буссинеска в безразмерном виде [1]:

V = 0,

д V

— + а1у (V ■ V - с) + Ур = г, (1) дТ

— + (у- v)т = (кут).

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

д U

M— = -W, dt

д

W = Е дцW- ф

1=1

(2)

где

U = H

(Р \

V1 V2 \ V3 J

W1 = H

v

Vi

2 + p - Til V1V2 - T12

\ V1V3 - T13 /

W2 = H

V2

V1V2 - T12 V22 + Р - T22 \ V2V3 - T23 )

W3

V3

V1V3 V2V3

T13 T23

\ V2 + p - T33 /

Ф

0 0

V2 + p - T33 V -V2V3 + T23

M

0000 0 10 0 0 0 10 V 0 0 0 1 /

dV1 dV2 2^dV3 . N

Tn = 2^д^, T22 = 2^д^, T33 = H d^ + j V2

dx1 dX2 H \ dX3 ,

T12 = V

' dV1 dV2

дж2 dx1

T13 = V

dV3 1 dV1

3 + 1

dx1

H дж3

/ dVs 1 dV2 .V3N

T23 = Чд*2 + Hd^ - jH

Здесь x, v , (/ = 1, 2, 3) — оси системы координат и компоненты вектора скорости v, П — коэффициент вязкости, H = rj — коэффициент Ламе. Значение параметра j = 0 соответствует декартовым координатам, где H =1, ж1 = ж, ж2 = у, ж3 = z и V1 = u, v2 = V, v3 = w; j = 1 — цилиндрическим, где H = r, x1 = z, x2 = r, ж3 = ^ и V1 = vz, v2 = vr, v3 = vv соответственно. Уравнения (2) записаны в безразмерном виде. При переходе к безразмерному виду величины нормировались на их невозмущенные значения, а в коэффициент вязкости внесено число Рейнольдса, т.е. принято п = По/Re. Отметим, что матрица M в системе уравнений (2) вырожденная. Уравнения (2) дополняются граничными и начальными условиями. Их постановку проиллюстрируем на задаче об осесимметричном течении за торцом цилиндра (рис. 1).

На поверхности тела заданы условия прилипания vz |г = vr |г = 0.

На входе, при z = 0, заданы условия набегающего потока vz = vz(r), vr =0, а на выходе при z = L задаются мягкие условия дf/dz = 0, где f = (p, vz,vr,v^)T.

Если течение симметричное, то во всей расчетной области vv = 0, д/д^ = 0, а на оси симметрии при r = 0 справедливы соотношения дvz^r = vr = vv = 0.

Верхняя граница задается достаточно далеко от обтекаемого тела и на ней ставятся условия набегающего потока: vz = 1, vr = 0,v^ = 0.

В начальный момент времени при t = 0 во всей области, кроме границ, заданы нулевое значение скорости, или vz = vz(r), и давление p = const.

3

j

Рис. 1. Расчетная область

Наряду с дивергентной формой уравнений Навье—Стокса (2) при построении алгоритма будет использован их недивергентный вид:

Л-Р

м— + ы = И, д^

(3)

где В = £ Бг, /=0

f =

( Р \

VI "2 V "з )

В0 =

0

д

д дХ2

\ ЯдХЗ

д 1 д

дх1 0

0

Я дх2 0

я

0

1 д \

Я дх3 0

0

И

0 /1

+3Я

."2 "з

+ 3

Я /

В1 =

00

д д д 0 --

дх1 дх1 дх 00

00

"1

д

дд п-

дх1 дх1 дх 0

"1

дх1

0 0

дд -п

дх1 дх1 /

0

В2 =

0 "2

0

0

А

дх2

0

1 д Яп д

Я дх2 0

0

дх2

"2

А

дх

0 0

А А

Я дх2 0

Яп

А

дх2

"2

дх

0 0

1А

Я дх

Яп

д

дх2 )

з

2

0

0

0

0

0

0

1

0

Вз

( 0 0

0 0

V3 _д_

H дх3

0

H2 дх3 0

0

V3

H дх3

0 1

д д -п-

H2 dx3 дх3 0

V3 д

0 0

2 д д

H dx3 H2 дх3 дх3 )

При записи уравнения (3) матричный оператор В представлен в виде расщепления по физическим процессам и пространственным направлениям. Матричные операторы Вг (/ = 1,..., 3) учитывают конвективные и вязкие члены (повторные производные) в уравнениях движения, оператор В0 содержит уравнение неразрывности и члены с давлением в уравнениях движения, а компоненты /1; /2, /3 вектора Ф — смешанные производные в уравнениях движения. Отметим эквивалентность представления уравнений в дивергентном (2) и недивергентном (3) видах, т. е. выполнения соотношений

Вf - R = W/H. Этот факт будет использован при построении разностных схем.

(4)

2. Численный алгоритм

Введем в расчетной области равномерную сетку с шагами Л-1,Л,2, Л,3 соответственно. Определим в узлах сетки г, к, I сеточные функции /п = /"¿-г и аппроксимируем первые и вторые производные симметричными операторами со вторым порядком по формулам

1

Л1/г , k ,Z = , k , Z — /г—1, k , Z )j

2hi

Л1^Л1/г, k,Z = Tö

ЛlWЛ2/г,k,í =

^i+1/2,k,Z (/i+1,k,1 — Ak,Z) — ^г— 1 /2,k,1 (/i,k,1 — Уг—1,k,1)

(5)

1

[^i+1,k,Z(/i+1,k+1,Z — /i+1,k—1,Z) — ^г— 1,k,Z(/i— 1,k+1,1 — /"г—1,k—1,Z)] •

4^2

Производные по другим направлениям аппроксимируются подобным образом. Аппроксимируем матричный оператор В и вектор W разностными операторами В^ и Wh со вторым порядком согласно (5). Для численного решения уравнений (2) рассмотрим разностную схему:

(3 \ .п+1 _ .п 1

М + та 1=0 ВП^-т-= _ н Wh, (6)

т

где стабилизирующий оператор в левой части схемы основан на аппроксимации уравнений Навье—Стокса в недивергентной форме, а в правой — в дивергентной. Это возможно в силу эквивалентности форм записи уравнений (4). Схема (6) линейна относительно (п + 1)-го слоя и аппроксимирует уравнения (2) с порядком 0(тт + ^2), где Л, = шах(^1, Л.2;Л,3), а т = 2 при а = 0.5.

Хотя матрица М в (6) и является вырожденной, справедливо приближенное равенство

3

М + таВ^ = (М + таВоЛ) П (I + таВ№) + 0(т2),

г=1

0

0

и тогда вместо схемы (6) может быть рассмотрена схема приближенной факторизации

3 .рга+1 _ гга 1

(М + таД*) п (I + таДь)-= _ 77Wh,

1=1 т 7

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

еп=_ -

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

(/ + таВ1^)ега+4/6 = еп+1/2, (7)

(/ + та^еп+5/6 = еп+4/6,

(/+таВз^)ега+1 = еп+5/6,

. п+1 = + теп+1

Остановимся на ее реализации. На нулевом шаге значение вектора невязок е = (ер, е1, е2,е3)Т во внутренних узлах области вычисляется явно из первого уравнения (7), где компоненты вектора Wh и матриц являются сеточными аппроксимациями соответствующих компонент вектора W и матриц Д.

На первом дробном шаге решается система разностных уравнений:

та (Л1еп+1/2 + -1л2яеп+1/2 + -1лзеп+1/^ = е

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

¡1 -й^шчГ

еп+1/2 = еп _ таЛ2еРп+1/2, еп+1/2 = еп _ - ЛзеРп+1/2-

Исключая значения еП+1/2 (^ = 1 2, 3) из первого уравнения, получим уравнение

относительно невязки для давления

тад^еп+1/2 = /;, (8)

/п = (л^ + -^-еп + —Лзезга) _ ¿еп, = Л1Л1 + -1Л2ЯЛ2 + —Л3Л3,

где Д^ — сеточный аналог оператора Лапласа. Его решение найдем как предельное решение нестационарного уравнения методом установления, например, по схеме приближенной факторизации:

(I _ таЛ1Л1) (/ _ -Л2-Л^ (/ _ —Л3Л^ Ср+\ ер = д^е; _ /П,

или эквивалентной ей схеме в дробных шагах, реализуемой скалярными прогонками по каждому пространственному направлению:

с * = дер _ /;,

(I _ таЛ1 Л1Х"+1/3 = С",

(/ _ —Л2ЯЛ2) с^+2/3 = с^+1/3, (/ _ —Л3Л3) С"+1 = С^+2/3,

ег1 = ер + тс "+1.

Итерации продолжаются до сходимости, т. е. до выполнения условия |е^+1 _ еУI — т8, где параметр 8 ~ О (Л2). Наряду с методом установления разностное уравнение (8) может быть решено и по какой-либо итерационной схеме. На втором дробном шаге система разностных уравнений

еп+4/6 = е п+1/2

[/ + та(^1 _ 2ЛщЛ1)]еп+4/6 = еп+1/2, [/ + таы Л1 _ лщЛ1)]еп+4/6 = еп+1/2, [/ + та(^1 _ лщЛ1)]еп+4/6 = еп+1/2

реализуется скалярными прогонками для каждой компоненты вектора скорости, а значение ерп+4/6 переносится с предыдущего шага. Аналогично решаются уравнения по второму и третьему направлениям. Наконец, новые значения функций вычисляются явно по формулам

<+1 = < + тег1, = + те2п+1, = ^^ + теп+1, рп+1=+ те;+1

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

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

двумерных и пространственных задач

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

в [9].

3.1. Моделирование двумерных течений

В первой серии расчетов в рамках полных уравнений Навье—Стокса в двумерном приближении исследовано течение вязкой несжимаемой жидкости в канале за уступом

(7 = 0) и за торцом цилиндра (7 = 1) при числах Ие = 100, 500,1000. Расчеты проведены для расчетной области, представленной на рис. 1, при следующих параметрах: высота расчетной области Н = 1, высота Я = 0.25 и длина Ь0 = 0.25 уступа в плоском случае или цилиндра в осесимметричном. Длина расчетной области при обтекании уступа задавалась равной Ь = 2.5, а при обтекании цилиндра выбиралось Ь = 1.0. Распределение скорости во входном сечении задавалось из решения Пуазейля [2] при различных числах Ие. Для оценки точности полученных решений по предложенному алгоритму были проведены расчеты на различных расчетных сетках. На рис. 2 приведены распределения продольной компоненты скорости ух при течении около торца цилиндра для Ие = 500 на сетках 40 х 80 (цифра 1), 60 х 120 (2) и 80 х 160 (3) в центральном сечении хх = 2Я = 0.5. Отметим, что распределения скоростей на сетках 60 х 120 и 80 х 160 совпадают с высокой точностью. Поэтому дальнейшие расчеты проводились на сетке 60 х 120. Типичная картина течения за торцом цилиндра для Ие = 500 приведена на рис. 3. Известно [1, 2], что длина возвратной зоны за телом существенно зависит от формы обтекаемого тела. Она меньше за цилиндром, чем за уступом, что объясняется растеканием потока за торцом цилиндра в азимутальном направлении. Для оценки длины области отрывных зон была проведена серия расчетов течений около торца цилиндра и уступа для различных чисел Ие. Результаты этих расчетов представлены на рис. 4, где кривые 1 и 2 соответственно отвечают плоскому и осесимметричному случаям. Из расчетов следует, что длина возвратной зоны за торцом цилиндра не меняется уже при числах Ие > 700.

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

(г*1 при 0 < г < гх, (г - гх)(гх - г2)/(*! - г2) + гх, г2 при г2 < г < 1.

Для повышения точности расчетов было введено преобразование координат q1 = q1 (г,г), q2 = q2(г,г), отображающее исходную расчетную область на единичный квадрат. Уравнения Навье—Стокса в преобразованных криволинейных координатах реша-

Рис. 2. Распределение скорости на различных сетках: 40 х 80 (1), 60 х 120 (2), 80 х 160 (3)

Рис. 3. Распределение вектора скорости за торцом цилиндра для Ив = 500

Рис. 4. Зависимость длины возвратной области за торцом пластины (!) и цилиндра (2)

Рис. 5. Расширяющийся канал

Рис. 6. Распределение вектора скорости в расширяющемся канале для ^ = arctg0.7 (а) и

^ = аг^ 2 (б)

а

лись по схеме приближенной факторизации. Вид уравнений в криволинейных координатах и разностная схема их решения приведены в Приложении. При проведении расчетов использовалось преобразование координат ^ = г, q2 = [г/г0(г)]Ь, что позволяло сгущать узлы сетки в окрестности пограничного слоя путем варьирования параметра Ь. Решение задач находилось при различных значениях итерационного параметра т, причем устойчивость алгоритма сохранялась в широком диапазоне его изменений. Расчеты проведены на сетке 100 х 100 узлов. Типичное число шагов до установления составляло « 103 ... 104 в зависимости от параметров течения и начального поля течения. При выборе в качестве начального распределения параметров потока (например, при другом числе Re) число шагов до установления уменьшалось в несколько раз. Ниже приведены некоторые результаты расчетов течений в расширяющемся канале. Для чисел Re = 102 вплоть до угла раствора канала ^ = arctg0.75 течение было безотрывным, а с увеличением чисел Рейнольдса и угла наклона канала длина области возвратного течения начинала возрастать. На рис. 6 приведены распределения поля скоростей в расширяющемся канале при угле наклона стенки ^ = аг^0.7 и ^ = аг^2.0 для Re = 500 и степени сгущения узлов Ь = 2.0.

Опытным путем показано, что оптимальное сгущение узлов сетки достигается в диапазоне 1.3 < Ь < 2.0, для него требуется максимальная точность при минимальном числе узлов сетки 100 х 100 для чисел 102 < Ие < 103.

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

3.2. Моделирование пространственных течений

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

Л 9У1 Т —- = 0,

1=1 ОХ1

dvj j dp dri

+ jx = £ dj + ^3Ra ■T j = 1' 2 3)' (9)

dt dxi = dx.

дТ дт д2т\

Ж + 1=1 V1 дХ - ~дХ!) =0'

где вид элементов тензора ту приведен в (2). При переходе к безразмерному виду параметры потока в исходных уравнениях нормировались, как это принято в геофизике (см., например, [12, 13]), на следующие величины:

/___ч (х у г \ г ъ П - Рь2

L'L'L/' KoL2 по' По ко'

- Т _ _ (uL vL wL \

T =-г^, (u,v,w)=[—,—,- ,

AT V ко ко ко /

где L = 1 — размер каверны, щ ,ко — коэффициенты вязкости и температуропроводности, AT — разность температур на границах области, p = ро — pgz+ const — отклонение давления от гидростатического, Ra = вродхL3AT/(поко) — число Рэлея, в — коэффициент теплового расширения, gz — компонента силы тяжести, направленная противоположно направлению оси z. Решение уравнений (9) отыскивалось в единичном кубе (рис. 7).

Полагалось, что через верхнюю и нижнюю поверхности куба жидкость не протекает, и на них задавались условия прилипания для скоростей и значения температуры, т. е. на поверхностях z = 0 и z =1 (0 < x < 1, 0 < y < 1) выполнялись условия u = v = w = 0, T(z = 0) = 1, T(z = 0) = 0.

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

— на поверхностях x = 0 и x = X (0 < y < 1, 0 < z < 1) u = dv/dx = dw/dx = 0;

— на поверхностях y = 0 и y = Y (0 < x < 1, 0 < z < 1) и = du/dy = dw/dy = 0.

Рис. 7. Расчетная область и постановка граничных условий

Краевые условия для температуры приведены на рис. 7. В начальный момент времени ставились нулевые значения скорости u = v = w = 0, значение давления полагалось постоянным и равным p0 = 0.4, а температура T0 (x,y, z) задавалась в виде (см. [13])

T0 (x, y, z) = 1 — z + a (cos (nx) + cos (ny)) sin (nz). (1)

Численное решение упрощенных уравнений неразрывности и движения в (9) находилось по разностной схеме (7), а для решения уравнения теплопроводности использовалась схема приближенной факторизации

3 Tn+1_тп 3

П [I + ™(viAi _ Л)]-= _ E [viAl _ All] Tn (2)

j=i T 1=1

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

СТ = - D [viAi - Л11 ] Tn, i=i

[I + ra(viЛ1 - Лц)] СТ+1/3 = СТ, [I + ra(v2Л2 - Л22)] СТ+2/3 = СТ+1/3, (3)

[I + ra(v3Л3 - Л33)] СТ+1 = СТ+2/3, Tn+1 = Tn + r^T+1.

Разностные схемы (11) или (12) аппроксимируют исходное уравнение при a = 0.5 + 0(т) со вторым порядком по всем переменным и реализуются на дробных шагах скалярными прогонками.

По предложенной модели и алгоритмам проведены многочисленные расчеты нестационарных течений жидкости в квадратной каверне с постоянной и переменной вязкостью (двумерный случай) и течения в кубической каверне при подогреве снизу. Стационарные решения данной задачи (при отсутствии конвективных членов) соответствуют равновесному неподвижному состоянию в широком диапазоне параметров — чисел Рэлея и коэффициентов вязкости (1 < Ra < 103, 10-3 < rj < 1). При этом картины течения в процессе установления в зависимости от вида начального распределения температуры (10) существенно различались и для двумерного, и для трехмерного случаев.

3

В двумерном случае при И,а = 1, п = 10-3, на момент времени I = 50 при начальном задании параметра а = 0 в (10) движение жидкости было двухвихревым (рис. 8, а), а при задании а = 0.2 содержало один вихрь с круговым движением вокруг центра полости (рис. 8, б). Подобный режим установления отмечен в работе [3]. При установлении оба режима выходили на неподвижное состояние с почти линейным распределением температуры. Для примера на рис. 9 показано распределение давления. С высокой точностью оно удовлетворяет уравнению для давления др/дг = И,а д2 (1 — z).

а б

Рис. 8. Функция тока: а — при а = 0; б — при а = 0.2

Рис. 9. Распределение давления в состоянии Рис. 10. Зависимость вектора скорости от чис-равновесия ла итераций

При проведении расчетов в трехмерной каверне расчетная сетка содержала, как правило, 20 узлов по каждому направлению, а временной шаг т менялся в зависимости от параметров потока (чисел Рэлея и вязкости). Расчеты течений были проведены для различных значений параметров потока. Для иллюстрации приведем результаты расчетов при И,а = 3 • 104, п =1 и начальном распределении температуры в (8) с параметром а = 0.2. График выхода на стационарный режим, соответствующий неподвижному состоянию, представлен на рис. 10. Отметим, что выход на стационарный режим с погрешностью О (К2) достигался за к 104 итераций, хотя для его подтверждения проводилось 2 • 104 итераций. Одна из важнейших характеристик алгоритма — его способность обеспечить сходимость итерационного процесса при численном интегрировании уравнения Пуассона для невязки давления (8). В начальный момент времени для установления уравнения Пуассона требовалось до 200 итераций, а затем число итераций до сходимости уменьшалось до 2-5. Отметим, что в процессе установления наблюдалась более сложная, по сравнению с двумерными задачами, конфигурация течения. На рис. 11 и 12, а приведены поля вектора скорости для трех сечений куба г = 0.5, у = 0.5 и х = 0.5 на момент времени I = 50. Отметим, что в плоскости куба при г = 0.5 течение жидкости стремится к центру (см. рис. 11, а), в плоскости у = 0.5 (рис. 11, б) возникает вихревое течение с круговым вращением в центре области и застойными зонами около верхней и нижней плоскостей г = 0 и г =1.

Еще более сложная картина течения имеет место в сечении х = 0.5 (см. рис. 12). В центре плоскости образуется центральный вытянутый вихрь большой интенсивности, а у верхней и нижней границ области — вихревые зоны малой интенсивности с вращением в противоположном к центральному вихрю направлении. Нормальная составляющая скорости представлена на рис. 12, б. Видно, что в поперечном направлении к плоскости х = 0.5 поток жидкости движется в разные стороны симметрично оси г = 0.5. На рис. 13 показаны изолинии отклонения температуры = Т — (1 — г) от линейного

Рис. 11. Поля вектора скорости: а — сечение г = 0.5, тах |V| тах |V| к 3.64

0.19; б — сечение у = 0.5,

а

Рис. 12. Сечение х = 0.5: а — поле скорости, тах |V| ~ 1.91; б — изолинии и

Рис. 13. Изолинии отклонения температуры

по высоте распределения, отвечающего равновесному неподвижному состоянию, в срединном сечении г = 0.5. В данном случае в правом нижнем углу сечения происходит подъем нагретой снизу жидкости, а в левом верхнем углу — опущение жидкости охлажденной. При установлении решение соответствует неподвижному состоянию, как и в двумерном случае. Отметим, что распределения температуры, полученные из решения задачи о течении жидкости в квадрате и в кубе, практически совпадают, что позволяет сделать вывод о достоверности полученных результатов.

Приложение

Уравнения Навье—Стокса в преобразованных координатах

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

В цилиндрических и декартовых координатах уравнения Навье—Стокса (2) в двумерном случае представляются в виде

М^ = -ТО", W 5*

V - Р

^ 5т- '

(13)

где

/

и = н

р

VI \

Wl = н

VI

V,2 + р - Гц VlV2 - Т12

= Н

V2

VlV2 - Т12 V V2 + р - Т22

Т11 =

5v1 5т 1

Е

М

/000 0 1 0 V 0 0 1

5vЛ ./ «2 \

ЧдТ2 + 5X1,1' 6 = Пр-

Введем невырожденное преобразование координат, переводящее расчетную область в единичный квадрат = 51(т1,т2), д2 = 52(т1,т2). Тогда существует обратное преобразование т1 = т1(д1,д2), т2 = т2(д1,д2), а исходные и новые переменные связаны соотношениями

5 5 5^1 5 дф _ 5 5

т;— = й— п--' п— п— = п--г ¿21 ^— = 1, 2).

5x1 5д15x1 5д2 5т 5д1 5д2

В преобразованных переменных уравнения (13) вновь могут быть записаны в дивергентном виде:

5И _

М— = -W, 5*

(14)

где

_ И _ 5W1 5 W2 — ___

И = -, W = —- + —^ - Е, Wl = (¿1^1 + ¿^2^0, = + ,

3 5^1 5^2

30 = ( дт-5X2 - 5X15X2 ) ,3 = "г--якобиан преобразования, W = (Ж1, Ж2, Ж3)т,

\ 5^1 5^2 5^2 5^ / 3о

- Н д д

Е = Е3о, 3 = -, Ж = —(3^1) + — (3^), 3 д^1 д^2

дд ^2 = — ЗД^ + ¿ир - СТц) + — 71(^2^1 + ¿21Р - СТ12),

дд Жз = — 31(^1^2 + ¿12Р - СТ21) + 31(^2^2 + ¿22Р - ^22) - ^'То,

VI = ¿11V + ¿^ (1 = 1, 2), 3 = р - 2^Н> °"11 = ¿11Т11 + ¿12Т12, аи = ¿21 ти + ¿22Т12,

п I дь1

021 = ¿11Т12 + ¿12Т22, ^22 = ¿21Т12 + ¿22Т22, ТИ = 2^ ¿ц---+ ¿21 —

V %

д^2

( д^о д^Л / дг>1 д^о д^' Т22 = 2^ ¿ю"--+ ¿22— , Т12 = И ¿12^--+ ¿22^--+ ¿ц---+ ¿21 —

V д<21 дq2J V д^1 д^ д^1 д^2.

В недивергентном виде уравнения (14) примут вид

где

В

д

г\.о

М— + В£ = о,

2 д

Е ¿11т;— 1=1 дqz

д

(15)

£¿11^ - Ь2)

1=1 дq1 1=1 дq1

2д

Е ¿12^ 0

\ 7=1 д®

Л ¿12 д тт \

1=1 Н д® 0

2д Е (Ит^ - Ь3) 1=1 дq1 3' /

д

д , ¿12 д

д

Ь2 = 2^1^— ^.¿ц — + — — —,

дq1 дq1 Н д®

дq1;

д

Ьз = ¿11^— ^¿11^- +

д 2¿12 д

д

• ^

Н8, ^ - <' 2)

£п+1 _ £ п

т

- Н ^

(16)

. гл гл

дq1 дq1

Оператор В содержит первые и повторные вторые производные, а вектор О — смешанные. Соответствующим образом в новых координатах переписываются и граничные условия. Подобно (2), (7) для численного решения системы уравнений (14), (15) рассмотрим разностную схему

(М + таВоь,)(1 + таВш)(/ + таВ2^) или эквивалентную ей схему в дробных шагах:

£п = ^3/Н, (М + таВо^)С+1/3 = С, (I + таВш)£п+2/3 = £п+1/3, (17)

(I + таВ2^)С+1 = С+2/3,

£ П+1 = £п +

аппроксимирующие исходные уравнения с тем же порядком 0(тт + Л2), что и в (2). Здесь

0 ¿ПЛ1+¿21Л2 нЛ1Н+нЛ2Н ^

¿11Л1 + ¿21Л2 0 0

\ ¿12Л1 + ¿22Л2 0 0 )

(

В

о^

0

2

( 0

B

lh

0 V^i - b1

0

0

0

2 0 v^ - b3

0

B

2h —

0

0

о - b2 о V о о У2Л2 - b3

а коэффициенты матричных операторов Bj заменены их разностными аналогами Bjh.

Остановимся на особенностях реализации схемы (17). Значение вектора находится явно из первого уравнения схемы, где компоненты вектора Wh являются сеточными аппроксимациями компонент вектора W:

Wih — Л1 ti +Л212, W2h — Л^^) + Л1 Ji(ziip - ail) + Л2(¿2Vl) + Л2 Jl(z2ip - CT12),

W3h — Л1 (tlV2) + Л1 Jl(zi2p - 021) + Л2(¿2^2) + Л2 Jl(z22P - ^22) - j Jo (p - '

ti — Ji Vl, а коэффициенты преобразования координат заменены их разностными ана-

логами:

Zii — J Л2Ж2, Z12 — - JЛ2Xl, Z21 — - JЛlX2, Z22 — JЛlXl.

На первом дробном шаге решается система уравнений:

*n+l/3

га(гпЛ1 + ^21Л2)СП+1/3 + та (fЛ1Я + f Л2Я) £Г1/3 — Cpn,

Я Я

£П+1/3 — СП - та(^пЛ1 + Z2^2)£pn+l/3, еП+1/3 — C2n - Та(^12Л1 + Z22Л2)en+1/3•

(18)

Исключая значения ^п+1/3 из первого уравнения (18), получим сеточный ана-

лог уравнения Пуассона:

я^+1/2 = /;, (19)

где

Rph — та

2 ' Z12

= ^иЛг^иЛг + f Лг Я^Лг J + zn Л^Л2+

+^21Л2^ПЛ1 + f ЛlЯz22Л2 + f Л2ЯZl2Лl

/n — (zi^ien+z^e? + Я л1я^2п+f Л2я^2п) - таС

n ep •

Его решение будем находить как предельное решение нестационарного уравнения, например, по схеме приближенной факторизации методом установления:

п

i=l

1 Zi2

I - та (^Лг^Лг + fЛгЯ^Лг

ev+i - e S>p S>i

R £v _ fn Rph?p / ,

или по эквивалентной ей схеме в дробных шагах, реализуемой скалярными прогонками по каждому направлению:

/V _ р tv _ in

z Rph?p / ,

zl2

I - та f zi^izi^i + fЛlЯzl2Лl Z22

I - та(z21 Л2Z2lЛ2 + f Л2ЯZ22Л2)

epv+1 — epv + тс

c v+i/2 — c v, Z v+l — Z v+1/2

*v+l

0

Итерации продолжаются до сходимости, т. е. до выполнения условия |£^+1 - £у | < т8, где параметр 8 ~ 10-3. Наряду с методом установления разностное уравнение (19) может быть решено и по итерационной схеме. На втором дробном шаге система разностных уравнений

1 + та (V^i - - ЯЛ^^Ль)

1 + та f V^i - - 2—Л^^Ль

¿n+2/3 _ ¿n+1/3

51 — ?1 >

¿n+2/3 _ ¿n+1/3

52 — S2

^п+2/3

реализуется скалярными прогонками для каждой компоненты вектора £1 . Аналогично решаются уравнения и на третьем шаге:

I + та ( V^2 - 2Z21 Л2^21Л2 - Z22Л2Я^^21Л2

I + та ( У2Л2 - ¿21Л2^^21Л2 - 2ЯЛ2Я^^22Л2 + 2jH

¿n+2/3 — S1 >

tn+1 _ ¿n+2/3 S2 — S2 •

Наконец, новые значения функций вычисляются явно по формулам

n+1

,.n+1

+ т^

■n+1

1,

,.n+1 _ n

— vn + т^

n+1 2,

p

n+1 n

— Pn + т£

n+1 p

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

1

1

2

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

[1] Ландау Л.Д., Лифшиц Е.М. Теоретическая физика. М.: Наука, 1988. Т. 6. Гидродинамика. 736 с.

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

[3] Тарунин Е.Л. Нелинейные задачи тепловой конвекции. Избранные труды. Пермь: Изд-во ПГУ, ПСИ, ПСИ МОСУ, ПССГК, 2002. 216 с.

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

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

[6] Толстых А.И. Компактные разностные схемы и их применение в задачах гидродинамики. М.: Наука, 1996. 230 с.

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

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

[9] Ковеня В.М. Об одном алгоритме решения уравнений Навье—Стокса вязкой несжимаемой жидкости // Вычисл. технологии. 2006. Т. 11, № 2. С. 39-51.

[10] Трубицын В.П., Рыков В.В., Трубицын А.П. Роль конвективных процессов при образовании высоковязкой континентальной литосферы // Вестник ОГГГГН РАН. 2001. № 4(19). http://www.scgis.ru/russian/cp1251 /h_dgggms/4-2001 /trubitsyn.htm

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

[12] Червов В.В. Моделирование трехмерной конвекции в мантии Земли с применением неявного метода расщепления по физическим процессам // Вычисл. технологии. 2006. Т. 11, № 4. С. 73-86.

[13] BussE F.H., Christensen U., Clever R. et al. 3D Convection at infinite prandtl number in cartesian geometry — a benchmark comparision // Geophys. Astrophys. Fluid Dynamics. 1993. Vol. 75. P. 39-59.

Поступила в редакцию 15 сентября 2008 г., в переработанном виде — 6 ноября 2008 г.