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

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

Том 14, № 1, 2009

Реализация метода расщепления

по физическим процессам для численного решения

*

трехмерных задач конвекции A. Ф. Воеводин

Институт гидродинамики им. М.А. Лаврентьева СО РАН, Новосибирск, Россия

e-mail: [email protected]

O. Н. Гончарова Алтайский государственный университет, Барнаул, Россия e-mail: [email protected]

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

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

1. Введение в метод расщепления по физическим процессам для задач конвекции

Классические уравнения конвекции с постоянной вязкостью, или уравнения Обербека— Буссинеска, применяются для изучения тепловой гравитационной конвекции в замкнутых областях с твердой непроницаемой границей (см., например, [1]). Начально-краевая задача для системы уравнений Обербека—Буссинеска состоит в определении скорости жидкости, давления, температуры, удовлетворяющих в области П c границей Е системе уравнений, которая в безразмерном виде выглядит следующим образом:

* Работа выполнена при финансовой поддержке гранта президента Российской Федерации для ведущих научных школ РФ № НШ 2260.2008.1.

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

vt + v ■ Vv = -Vp' + — Av + F,

1

(1)

Re

V = 0, (2)

Т* + V ■ УТ = (3)

* ИеРг ^ '

Здесь V — безразмерная скорость, р — безразмерное модифицированное давление

(отклонение от гидростатического давления), Т — безразмерная температура, Е =

— (Ог/И,е2^0Т — сила плавучести, go — единичный вектор ускорения силы тяжести,

Ог, Ие, Рг — числа Грасгофа, Рейнольдса и Прандтля соответственно.

Для данной системы уравнений рассматривается начально-краевая задача. Начальное состояние при Ь = 0 определяется состоянием покоя нагретой по определенному закону жидкости: V = 0, Т = Т0(х) (х £ П). Граничные условия для вектора скорости выражают условия прилипания: V = 0 (х £ Е,Ь £ [0,Ь*]). Для температуры может быть рассмотрено на границе одно из условий первого, второго или третьего рода, когда на границе области течения известна температура или тепловой поток либо определен теплообмен с внешней средой:

дТ дТ

Т = Тв(х,Ь), — = Тв(х,Ь), — + а(Т — Тв) = 0, х £ Е, Ь £ [0,Ь*]. дп дп

Здесь а — коэффициент, характеризующий теплообмен с внешней средой, Тв — температура внешней среды или заданный тепловой поток на границе.

Так называемые конвективный и диффузионный переносы естественным образом выделяются как в уравнениях движения (1), так и в уравнении переноса тепла (3). Тогда вполне обосновано применение расщепления по физическим процессам в уравнениях на каждом временном слое. Идея метода расщепления по физическим процессам базируется на методе слабой аппроксимации [2] и аддитивности этих процессов для достаточно малых шагов по времени. С математической точки зрения расщепление разностного уравнения на составляющие, а затем обоснование аддитивности процессов, описываемых отдельными частями, рассматриваются в [3], где проводится доказательство суммарной аппроксимации исходного уравнения вследствие расщепления. Общая теория расщепления наиболее полно изложена в [4]. Явная схема расщепления по физическим факторам в задачах гидродинамики используется в работах [5-7] и заключается в трехшаговой реализации, включающей расчет давления.

Особенности разрабатываемого авторами метода расщепления для двумерных задач конвекции излагаются в [8-10]. Следует отметить, что в отличие, например, от задач протекания, в задачах конвекции в замкнутых областях и при длительных процессах во времени необходимо точно выполнять условия прилипания и непротекания на твердых стенках, гарантировать сохранение свойства соленоидальности поля скоростей и его энергетическую нейтральность (свойство сохранения среднеквадратичной нормы скорости при переходе со слоя на слой) [9, 10]. Только тогда представляется возможным обеспечить сохранение объема и контролировать выполнение закона сохранения энергии. Схема расщепления, предлагаемая в работе, является наиболее физически оправданной и отвечает указанным требованиям. Очевидно, что выделение этапа конвекции в исходных уравнениях позволяет обеспечить корректность расщепления с точки зрения удовлетворения граничным условиям. Из условий прилипания и гиперболичности системы уравнений следует, что граничные условия на этапе конвекции — следствие уравнений. Достаточной аргументацией выделения этапа конвекции является также и последующая реализация этапа диффузии в переменных "вихрь—функция тока" для двумерных задач и "ротор скорости—векторный потенциал скорости" для трехмерных

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

В данной работе численный метод, основанный на идее расщепления по физическим процессам, предлагается для исследования конвективных движений жидкости в трехмерных областях с твердыми непроницаемыми границами. Расщепление на конвективный и диффузионный переносы проводится для уравнений конвекции, записанных в физических переменных. Этап конвекции реализуется для компонент вспомогательной функции, условно названной конвективной скоростью. На этапе диффузии осуществляется переход к новым искомым функциям: ротору скорости и векторному потенциалу скорости. Применением оператора "ротор" на этапе диффузии исключается расчет давления. Поскольку вычисляемыми функциями являются как компоненты скорости, так и новые искомые функции, принципиальный момент в организации расчета — введение смещенных сеток (см., например, [11]). Смещенные сетки вида (п',т,1), (п,т',1), (п, т, I') вводятся для расчета соответствующих компонент векторного потенциала и ротора скорости. Сетки (п,т',1'), (п',т,1') и (п',т',1) используются для расчета компонент вспомогательной конвективной скорости, а (п, т, I) — для расчета температуры. Смещенные сетки позволяют проследить за автоматическим выполнением условия несжимаемости поля скоростей и организовать восстановление скорости через векторный потенциал. Для реализации этапа диффузии строится алгоритм прогонки с параметрами [9, 12, 13]. При написании граничных условий для новых искомых функций будем исходить из выполнения условий прилипания на твердых непроницаемых границах. Эти условия используются непосредственно для компонент ротора скорости, а для векторного потенциала используются условия, выражающие равенство нулю касательных составляющих векторного потенциала и производной по нормали его нормальной составляющей [14-16]. Граничные условия прилипания для физической скорости и их следствия выполняются на этапе диффузии с учетом, однако, смещенности некоторых границ.

Метод расщепления по физическим процессам в трехмерных задачах тестируется на известных задачах о свободной конвекции в кубической полости или параллелепипеде при подогреве одной грани [16-18].

2. Постановка разностной задачи

Начально-краевую задачу для системы уравнений Обербека—Буссинеска (1)-(3) будем рассматривать в параллелепипеде П = {[0, х [0,уо] х [0,г0]} с границей Е, определяемой соотношениями: х = 0, х = хо, у = 0, у = у0, г = 0, г = г0. Считая конвективный и диффузионный процессы переноса количества движения аддитивными на малых временных интервалах ¿к < £ < ¿к+1 (т = ¿к+1 — ¿к), сформулируем задачу конвективного движения жидкости как последовательную реализацию двух вышеназванных этапов.

Перепишем уравнения движения (1) в виде

V — v

+ K V = 0,

(4)

V — V

—Vp' + DV + F.

(5)

Аналогичное расщепление на этапы проводится и в уравнении переноса тепла (3).

Применим к уравнению (5) операцию rotor и запишем его в терминах "векторный потенциал (Ф)—ротор скорости (W)":

W — W

DW + rot F.

т

При этом выполняется

ДФ = —W.

(6)

(7)

Здесь K — оператор конвективного переноса K = v • V, D — оператор диффузионного переноса D = z/Д, V = Re- . При вычислении температуры на этапе диффузии полагается V = (RePr)-1. Кроме того, v = rot Ф, W = rot v, W = rot v, и выполнено div Ф = 0. Опустим "крышку" при обозначении искомых функций на этапе диффузии и представим искомые функции в компонентах:

W = (wl,w2,w3)

v = (vi,v2,v3)

dv3 dv2 dvl dv3 dv2 dvl dy dz^ dz dx} dx dy

дфз дф2 d^i дфз дф2 d^i

дУ

dz dz dx dx

дУ

Ф = (ф1,ф2,ф3).

Для функций /ш1,/ш2,/ш3 и ф1,ф2,ф3 имеют место следующие граничные условия:

^дфф i д 2ф2 x = 0, x = 1 : —— = 0, wl = 0, ф2 = 0, w2 = —-——, ф3 = 0, w3 = —

У = 0, У = 1 :

z = 0, z = 1 :

dx

дф2 дУ дфз

= 0, w2 = 0, фф3 = 0, w3 = —

dx2

д 2фз дУ2

dz

= 0, w3 = 0, ф1 = 0, wl = —

д 2ф1 dz2

ф1 = 0, wl = —

ф2 = 0, w2 = —

д 2фз дx2

д 2Ф1 ;

дУ2 д 2ф2 дz2

(8)

3. Реализация этапа конвекции

При построении разностной аппроксимации оператор конвективного переноса К из (4) будем рассматривать в виде суммы одномерных кососимметричных операторов, результат действия которых на искомую функцию f (одну из компонент условной конвективной скорости или условную температуру на этапе конвекции) следующий: Kf = К^ + К^ + К^ [4, 9, 10], где

1 ( дf , 6vulf

Kl f =o v^ +

дx дx

1( 9f , дV2f

K2f = 2 Г дУ + дУ

K f 1 ( дf + дv3f\

Kf =2 ГTz +-дГ)' (9)

Начальные условия для компонент конвективной скорости V = (Ъ1,Ъ2,Ъ3) задаются с предыдущего временного слоя и представляют собой соответствующую компоненту диффузионной скорости, рассчитанную на нужной смещенной сетке, а именно: на сетке (п, т', I') — для первой компоненты, на сетке (п', т, I') — для второй компоненты и на сетке (п',т',1) — для третьей компоненты. Если хп = (п — 1)кх, п = 1, 2,N1 (М1 = N + 1), то хп = хп + кх/2, п' = 1, 2,..., N. Аналогично определяются ут', гу.

"Крышка" для обозначения диффузионной скорости опускается, и используются следующие разностные представления компонент диффузионной скорости через соответствующие компоненты векторного потенциала скорости:

^1п,т'Л'

"Фзп.

т+1,1

—Ф

3п,т,1'

Ф

2 п,т' ,1+1

— ф2п.

т',1

,т,1'

^3п',т',1

ку кг

ф1п,т,1+1 — ф1п' ,т,1 ф3 п+1,т,1' — ф3п,т,1'

кг кх

ф2п+1,т',1 ~ ф2п,т,1 ф1 п',т+1,1 — ф1п' ,т,1

кх

К

Заметим, что при использовании указанных конечно-разностных аппроксимаций тождественно выполняется ^у V = 0 на сетке (п',т',1') (см. также [11]).

Для решения уравнения (4) и вычисления компонент вспомогательной скорости на этапе конвекции применяется двуциклический способ расщепления по направлениям [4]. С использованием двуциклической процедуры в случае некоммутативности операторов (здесь К) обеспечивается второй порядок аппроксимации по т [4].

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

Г+6 — fк

т/2

?к+6 + fк + 1 2+1

0,

Г+6 — f

т/2

к+6

+ К

Iк+6 + f

к+6

0,

f к+6 — f т/2

к+6

+ К

f к+6 + f

к+6

0,

f к+6 — f т/2

к+6

+ К

f к+6 + f

к+6

f к+6 — f т/2

к+ 6

+ К2\

к+ 6

2

0,

fк+1 — f т/2

к+6

+ К*(

V f к+1 + f к+6

2

0. (10)

В этих формулах К^К^ ,К3\ представляют собой конечно-разностную аппроксимацию конвективных операторов (9). Во внутренних точках смещенных сеток рассматривается аппроксимация, использующая центральные разности (см., например, [8, 9]). В граничных точках применяется аппроксимация направленными разностями. Заметим, однако, что при вычислении компонент конвективной скорости в коэффициентах учитываются условия прилипания при попадании на физическую границу области. Используемая при построении схемы аппроксимация конвективных операторов является полной, имеет второй порядок и обеспечивает кососимметричность разностных операторов конвекции [4, 9, 10]. Полученные системы разностных уравнений в (10) решаются с помощью метода прогонки. Хотя для матрицы системы и не выполнены условия диагонального преобладания, метод прогонки устойчив [19].

2

0

2

2

Отметим еще раз, что для расчета температуры Т как на этапе конвекции, так и на этапе диффузии в области П вводится основная разностная сетка (п, т, I). При этом сохраняется методика расчета температуры, основанная на последовательной реализации этапов конвекции и диффузии в уравнении переноса тепла (3).

4. Реализация этапа диффузии

Для реализации этапа диффузии (6), (7) используются смещенные сетки (п',т,1) для расчета первых компонент ротора скорости и1 и векторного потенциала фф1, (п, т', I) — для расчета вторых компонент и2, ф2 и (п,т,1') — для расчета третьих компонент и3, ф3. Заметим, что при переходе на этап диффузии ротор конвективной скорости рассчитывается как раз на нужных вспомогательных сетках и имеют место следующие представления ("волна" опускается):

^3п',т',1 ^3п' ,т' — 1,1 ^2п' ,т,1' ^2п' ,т,1' — 1

и1п',т,1

и2п ,т' ,1 изп ,т,1'

Ну

^1 п,т' ,1' ^1п,т',1' — 1 ^3п',т',1 V3п' — 1,т',1

кх Нх

^2п' ,т,1' V2п' — 1,т,1' ^1п,т' ,1' ^1п,т' — 1,1'

Нх Ну

Отметим также, что равенство divW = 0 тождественно выполняется на сетке (п,т,1). Именно так называемый конвективный вихрь и берется в качестве начального условия для этапа диффузии к = Л¥.

В работе используется модификация метода Гаусса (метод прогонки с параметром) как наиболее экономичная по числу операций [13]. Алгоритм прогонки с параметрами на этапе диффузии представим на примере первых конпонент ротора скорости и векторного потенциала и1, ф1 (индекс 1 для удобства изложения опустим).

Пусть имеет место следующее представление первой компоненты ротора скорости:

ип,т,1 = °ип',т,1 +^^',1,1 + вт Шп' ,М 1,1 + ООЦШп' ,т,1 + вЬ^п'^М- (И)

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

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

В1В2В3 ип',т,1= Рп',т,1 (п' = 1,..., N; т = 2,..., М; I = 2,..., Ь),

Р'п'тЛ = В1В2В3ип',т,1 + Т!п',т,Ь

последовательные этапы реализации которой состоят в расщеплении по направлениям

ООО

и итоговом вычислении и на последнем этапе. Имеет место следующая последовательность (очередность) прогонок по направлениям:

1) прогонка по z:

B3 wn',m,i= Fn'imii, W (n,m,i) = 0 (i = {1; LI}),

= BiB2B3Wn'm,i + Tfn',m,l;

2) прогонка по y:

B2 Wn'mi=Wn',m,i, W (n, i, l) = 0 (i = {1; M1});

3) прогонка по x:

ooo oo ooo . , / Гт\тЧ\

Bi wn',m,i=Wn',m,i, w (i,m,l) = 0(i = {1; N}). (12)

Здесь введены обозначения:

t t d2

вг = e - вг = e + /^л*, л* - dxx2 (i = 1' 2> 3)-

Обозначения x* при i = 1, 2, 3 введены для x, y, z соответственно.

Как и в двумерном случае (см. [12]), прежде всего и только один раз могут быть найдены массивы коэффициентов a, ß при параметрах задачи. Эти массивы находятся в результате решения конечно-разностных задач вида

B2a2m = 0, ai = 1, aMi = 0; B2ßl = 0, ß2 = 0, ß*Ml = 1,

ВзаЗ = 0, a\ = 1, ah = 0; B3ß? = 0, ß3 = 0, ßh = 1.

Здесь, как и в двумерном случае, возможность вышеизложенного нахождения a*, ß* обеспечивается перестановочностью разностных операторов B

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

B2BzBi^n',m,i = Fn'mi (n' = 1,..., N; m = 2,..., M; l = 2, ...,L),

Fn',m,i = B2B3B^n',m,l + ^T W°n',m,i -Здесь введены обозначения

Bi = E - д2Ai, B2 = E - д2A2, B3 = E - д2A3, Bi = E + Д2Ai, B2 = E + A2A2, B3 = E + A2аз,

где

(A^)n',m,i = ^^)n',m,i,

аф)^^ = Лф)^^ - am^p]n',i,i - ßm[Л2Ф]п,Мi,i,

(АзФ)п',m,i = - a [Лзф]n',m,l - ßi [ЛзФ]п',m,Li

и A — итерационный параметр.

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

Итак, реализация схемы расчета ф выглядит как последовательное осуществление следующих этапов:

1) первый этап состоит в прогонке с параметрами для решения системы алгебраических уравнений

В2*фп',т,1 = Еп',т,1, фп',г,ь = 0 (г = {1; М1});

2) второй этап расчета состоит в прогонке с параметрами для решения системы уравнений

В3фп',т,1 = фп',т,1, фп'тг = 0 (г = {1; Ь1});

3) третий этап состоит в обычной прогонке для решения системы уравнений

В1фп',т,1 = фп',т,1, (фх)г,т,1 = 0 (г = {1; N}). (13)

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

При реализации граничных условий для вихря и (в данном случае — это первая компонента ротора скорости), например, при у = 0, используется конечно-разностная

д 2ф

аппроксимация условия и = —— из (8). Учитывается также, что при использовании

ду2

для расчета смещенных сеток требуется задавать граничные условия для компонент векторного потенциала и ротора скорости на смещенных границах. Для первых компонент искомых функций смещаться будут границы х = 0 и х = х0. Реализация гранич-дф

ных условий —— = 0, и = 0, заданных при х = 0 и х = х0 (см. (8)), осуществляется на дх

смещенных границах (т. е. при значениях индекса п' =1, N) с помощью представлений фо,т,1 = ф1,т,1 и и0,т,1 = —и1,т,1 для п' = 1 и аналогичных для п' = N.

В трехмерном случае реализация прогонки с параметрами — процедура более сложная, нежели в двумерном случае (см. [8, 12]). Поскольку принципиальный характер приобретает вопрос очередности прогонок по направлениям, повторим эту очередность для всех компонент искомых функций с указанием сетки расчета. Итак, при расчете и1 и ф1 используется сетка (п', т, I) и устанавливается очередность прогонок по г, у, х соответственно для расчета и и по у, г, х — для расчета ф1. При расчете и2 и ф2 используется сетка (п,т',1) и устанавливается очередность прогонок по х, г, у для расчета и и по г, х, у — для расчета ф2 соответственно. При расчете и3 и ф3 используется сетка (п,т,1') и устанавливается очередность прогонок по у, х, г для расчета

ооо

и и по х, у, г — для расчета ф3.

Общая схема решения задачи состоит в осуществлении следующих этапов:

1) переход на новый временной слой tk+1 начинается с расчета температуры Tk+l на основной сетке (n,m,l) по схеме (10) для этапа конвекции и с помощью, например, метода стабилизирующей поправки [2, 3] на этапе диффузии;

2) при найденном Tk+l рассчитываются компоненты конвективной скорости vk+l по схеме (10) и компоненты ротора конвективной скорости на вспомогательных сетках;

3) вычисляются составляющие компоненты диффузионного вихря (первое слагаемое в представлении (11) ) по схеме (12). На каждом k-м временном слое расчета вводится внутренний итерационный процесс расчета компонент векторного потенциала по схеме (13) с соответствующими граничными условиями. После окончания итераций при s = S считается, что с заданной точностью Еф определены значения компонент векторного потенциала на (k + 1)-м временном слое. Окончательно, с использованием представления (11), будут вычислены и компоненты диффузионного вихря;

4) при переходе на следующий временной слой (к этапу 1) насчитывается искомая диффузионная скорость.

Остается заметить, что итерационный процесс для ф (компонента векторного потенциала) считается сошедшимся, если выполнен критерий сходимости, например, следующего вида [14, 20]:

maX 1фП+т,1 - фПш!1 < ЕФ • max 1ф&п+шА,

п,т,1 ' ' п,т,1 ' '

где в — номер итерации, £ф — заданная точность расчета фя+1.

5. Результаты численного анализа

Предложенный метод расщепления для решения трехмерных задач конвекции проверяется на тесте о свободной конвекции в замкнутой кювете П = {[0, х0] х [0, у0] х [0, г0]} при подогреве одной грани х = х0 и теплоизоляции остальных [16, 17]: Т = 0 при х = 0, дТ

Т =1 при х = х0, —— = 0 при у = 0, у = у0, г = 0, г = г0. При этом будем следить за дп

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

¿0 У0

1 г г дТ

=- ^Итеап(г)(г, )= Ш1ос(у,г)(1у, Шгос(у, г) = — |

Уо го ] ] дх

о о

и х = 0 или х = х0. Аппроксимации второго порядка для производной температуры Т на границе используются при вычислении N^0,-. Тест состоит в достижении стационарного решения.

1. Тест 1 ([16])

Сила тяжести направлена против оси Оу, а область П представляет собой куб, размеры которого определяются значениями х0 = 1, у0 = 1, г0 = 1. При тестировании полагается х = х0, Ог = 105, Рг = 1, И,е = 1. Расчеты ведутся на сетках 21 х 21 х 21, 41 х 41 х 41, 81 х 81 х 81 с шагом по времени (фиктивному времени) т = 0.0001. Условием выхода на стационарный режим является выполнение неравенства [9, 20]

^и^1 - Ми^ | < еки,

где = £ф = 10-5.

Т а б л и ц а 1 Таблица2

Источник Сетка Сетка К

[16] 86 х 65 х 65 4.34 21 х 21 х 21 40.83

[18] 62 х 62 х 62 4.36 41 х 41 х 41 31.13

Настоящая работа 81 х 81 х 81 4.6342 81 х 81 х 81 28.56

Сравнение расчетов с известными результатами представлено в табл. 1.

На сетке 21 х 21 х 21 расчеты, демонстрирующие выход практически на один и тот же стационарный режим, проведены также с разными значениями т. Получены значения интегральной величины Nuoг,, равные 6.21, 6.57, 6.63 при соответствующих значениях временного шага 10-3, 10-4 и 10-5 и = 10-4т.

Устойчивость алгоритма и порядок сходимости вычислительного алгоритма проверяются путем вычислительных экспериментов на последовательности сеток 21 х 21 х 21, 41 х 41 х 41, 81 х 81 х 81. Наблюдения ведутся за интегральной величиной К (интегральной нормой скорости), представляющей собой аналог кинетической энергии. Результаты тестовых расчетов приведены в табл. 2. Экспериментальный порядок сходимости, близкий ко второму г ~ 1.92, и приближенно определяемая относительная погрешность е ~ 3.24 % вычисляются согласно правилу Рунге [21].

Значения интегральных чисел Нуссельта Nuoг,, рассчитанных на последовательности пространственных сеток 21 х 21 х 21, 41 х 41 х 41, 81 х 81 х 81 и характеризующих достигнутый стационарный режим, соответственно таковы: 6.5767, 5.0256, 4.6342.

Проведено дополнительное численное моделирование переходного режима и повторный выход на стационарный режим после ступенчатого нагрева стенки Хо = 1, вида:

- Т = И + Ь, где к = (Т - Т0)/(Ь1 - ¿о), Ь = (ТЛ - Т^/^ - ¿о), если ¿о < I < ¿1, и к = (Т - То)/(г2 - ¿3), Ь = (Тог - Т1з)/{12 - ¿3), если Ь <1 < ¿з;

— Т = Т, если ¿1 < I < ¿2, и Т = То, если I > ¿з.

Расчеты, проведенные на сетке 21 х 21 х 21 с различной амплитудой Т (Т = 1.5, 2), показали повторное достижение стационарного состояния, описанного выше. Так, при Т = 1.5 (То = 1), ¿1 - ¿о = ¿2 - ¿1 = ¿з - ¿2 = 5 получены значения Nuoг, = 6.571, К = 40.47.

Изолинии температуры в плоскости г = го/2 представляют собой типичную картину для поля температуры в двумерном случае (на рис. 1 представлены изотермы в плоскости г = го/2; по горизонтали — ось Ох, по вертикали — ось Оу).

Изотермы в плоскости г = го/2

ТаблицаЗ Таблица4

Источник фшах Источник Фшах

[17] 7.37 18.4 [17] 5.37 13.2

Настоящая работа 7.41 17.8 Настоящая работа 5.17 13.1

2. Тест 2 ([17])

Сила тяжести направлена по оси О г. Область течения представляет собой куб, определяемый размерами х0 = 1, у0 = 1, г0 = 1. При тестировании полагается х = 0, Иа = 5 ■ 105 (число Рэлея), Рг = 0.71, Ие = 1/Рг.

Таблица 3 представляет сравнение результатов с расчетами из работы [17], где отмечено, что большинство расчетов проводились на сетке 15 х 15 х 15, а некоторые — на сетке 51 х 51 х 51. Для сравнения берется не только число Нуссельта, но и максимальное значение величины 1ф21п,т,1 при у = у0/2 (фшах). Расчетная сетка — 41 х 41 х 41.

3. Тест 3 ([17])

Как и в предыдущем случае, сила тяжести направлена по оси Ог, а область течения представляет собой прямоугольный параллелепипед с размерами, определяемыми значениями х0 = 1, у0 = 2, г0 = 1. Расчетная сетка остается прежней: 41 х 41 х 41. При тестировании полагается х = 0, Иа = 1.5 ■ 105, Рг = 1, Ие = 1.

В табл. 4 представлены результаты в сравнении с расчетами, проведенными в [17].

Заключение

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

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

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

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

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

[1] Джозеф Д. Устойчивость движений жидкости. М.: Мир, 1981. 196 с.

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

[3] Самарский А.А. О принципе аддитивности для построения экономичных разностных схем // Докл. АН СССР. 1965. Т. 165, № 6. С. 1253-1256.

[4] Марчук Г.И. Методы расщепления. М.: Наука, 1988. 264 с.

[5] Белоцерковский О.М., Гущин В.А., ЩЕнников В.В. Метод расщепления в применении к решению задач динамики вязкой несжимаемой жидкости // Журн. вычисл. математики и мат. физики. 1975. Т. 15, № 1. С. 197-207.

[6] Белоцерковский О.М., Гущин В.А. Моделирование некоторых течений вязкой жидкости. М.: ВЦ АН СССР, 1982. 66 с.

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

[8] ВоЕводин А.Ф., Юшкова Т.В. Численный метод решения начально-краевых задач для уравнений Навье—Стокса в замкнутых областях на основе метода расщепления // Сиб. журн. вычисл. мат. 1999. Т. 2, № 4. С. 321-332.

[9] ВоЕводин А.Ф., Протопопова Т.В. Метод расчета вязких течений в замкнутых областях // Сиб. журн. индустр. мат. 2001. Т. 4, № 1. С. 29-37.

[10] ВоЕводин А.Ф., Гончарова О.Н. Метод расщепления по физическим процессам для расчета задач конвекции // Мат. моделирование. 2001. Т. 13, № 5. С. 90-96.

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

[12] ВоЕводин А.Ф., Гончарова О.Н. Метод расчета двумерных задач конвекции на основе расщепления по физическим процессам // Вычисл. технологии. 2002. Т. 7, № 1. С. 66-74.

[13] ВоЕводин А.Ф., Шугрин С.М. Численные методы решения одномерных систем. Новосибирск: Наука, 1981. 368 с.

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

[15] Исмаил-заде А.Т., Короткий А.И., Наймарк Б.М., Цепелев И.А. Численное моделирование трехмерных вязких течений под воздействием гравитационных и тепловых эффектов // Журн. вычисл. математики и мат. физики. 2001. Т. 41, № 9. С. 1399-1415.

[16] Бессонов О.А., Брайловская В.А., Никитин С.А., Полежаев В.И. Тест для численных решений трехмерной задачи о естественной конвекции в кубической полости // Мат. моделирование. 1999. Т. 11, № 12. С. 51-58.

[17] Davis G. de Vahl. Natural convection of air in a square cavity: a bench mark numerical solution // Int. J. Numer. Methods Fluids. 1983. Vol. 3. P. 249-264.

[18] МызниковА Б.И., Тарунин Е.Л. Процессы установления стационарных конвективных течений в кубической полости при подогреве снизу // Нестационарные процессы в жидкостях и твердых телах. Свердловск: УНЦ АН СССР, 1983. С. 20-29.

[19] ВоЕводин А.Ф. О корректности метода прогонки для разностных уравнений // Численные методы механики сплошных сред: Сб. науч. тр. / АН СССР. Сиб. отд-ние. ИТПМ. 1972. Т. 3, № 5. С. 17-26.

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

[21] Калиткин Н.Н. Численные методы. М.: Наука, 1978. 512 с.

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