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

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

Том 7, № 1, 2002

МЕТОД РАСЧЕТА ДВУМЕРНЫХ ЗАДАЧ КОНВЕКЦИИ НА ОСНОВЕ РАСЩЕПЛЕНИЯ ПО ФИЗИЧЕСКИМ ПРОЦЕССАМ *

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

Новосибирск, Россия e-mail: voevodin@hydro.nsc.ru О. Н. Гончарова Алтайский государственный университет, Барнаул, Россия Институт гидродинамики им. М.А. Лаврентьева СО РАН

Новосибирск, Россия

On the basis of method of splitting into physical processes the numerical method is suggested for investigation of convective fluid flows. It is characterized by energetic neutrality of velocity field. The splitting into two steps (convection and diffusion) is developed in convection equations written in original variables. Convection step is realized for the velocity. The conversion to new functions (vorticity and stream function) is used on the diffusion step. It eliminates the pressure gradient calculation and guarantees the solenoidality of velocity vector. The method is tested on the well-known problems on free convection in square cavity with one-side heating.

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

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

V + V •VV = -Ур' + т^AV + ?, (1)

Ке

а1у V = о, (2)

Тг + V ■ УТ (3)

4 КеРг ' к 1

* Работа выполнена при поддержке Совета Программы поддержки ведущих научных школ, грант №0015-96162.

© А. Ф. Воеводин, О. Н. Гончарова, 2002.

начальным условиям при Ь = 0:

и граничным условиям на Е:

V = 0, Т = То(х)

Прандтля; ^ = — ^g0T — единичный вектор ускорения силы тяжести). В дальнейшем

дТ

V = 0, Т = Тв(х, Ь) или — = 0.

дп

Здесь система уравнений выписана в безразмерном виде. В статье приняты следующие обозначения: х — вектор пространственной переменной; Ь — время; р' — модифицированное давление (отклонение от гидростатического); Ке, Сг, Рг — числа Рейнольдса, Грасгофа, _Сг

Ке2 .

будем предполагать, что сила тяжести направлена против оси Оу (^0 = — _)).

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

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

Рассмотрим задачу для системы уравнений (1) - (3) в прямоугольной области П. Эти уравнения перепишем в виде

д V

— + XV = —Ур' + DV + f, (4)

дЬ

V = 0, (5)

дТ

— + КТ = DT. (6)

Здесь К = V•У — оператор конвективного переноса; D = рД — оператор диффузионного

1 -1ч

переноса (V = — или V =-).

Р ^ Ке Ке Рг

Рассмотрим уравнение (4). Выделим два различных процесса переноса количества движения: конвективный и диффузионный. Считая эти процессы аддитивными на малых временных интервалах Ьк < Ь < Ьк+1 (т = Ьк+1 — Ьк), сформулируем задачу конвективного движения жидкости, исходя из нижеследующих двух этапов: 1. Этап конвекции

V — V

--- = —К]/ (7)

т

(V — "конвективная скорость").

2. Этап диффузии

V - V

-= —Vp' + DF + f (8)

т

(У — "диффузионная скорость"). Последний этап перепишем в терминах "функции тока — вихря", применив операцию "rotor":

& — &

-= D& + rotf,

т

Д* = —&, (9)

при этом на :

* = 0, dí = О, (10)

где & = rotF; & = rotV; V = (U,«2); U = -7—; u2 = — тт-; divV = 0.

dy dx

При построении разностной аппроксимации оператор конвективного переноса будем рассматривать в виде суммы двух одномерных кососимметричных операторов [3, 4, 6]:

= KiV + K2V,

где

^ =2 (¿1 ^-fl^J, K2V =2 (¿2^ + (11)

В области П вводим основную разностную сетку (xn,ym) для расчета температуры и моделирования этапа диффузии (задача (9), (10)):

xn = (n — 1)hx, n = 1, 2, ..., N, = N, N = NV + 1,

ym = (m — 1)hy, m = 1, 2, ..., M, hy = |0, M = M+1

Уо

МТ

и две вспомогательных аналогично [7]: (жп,уш/), ш' = ш+1/2, п = 1, ..., N, ш' = 1, ..., М — для расчета первой компоненты конвективной скорости У и (жп/, уш), п' = п +1/2, п' = 1, ..., N, ш = 1, ..., М — для расчета второй компоненты конвективной скорости г;2 (V = (^1,^2)).

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

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

Фп,ш+1 - Фпт Фп+1 ,т Т'пт

_ /n,m+1

v1n,m+1/2 = Т j v2n+1/2,m

, ' "2n+1/2,m 7

hy hx

Для решения уравнения (7) и вычисления температуры на этапе конвекции используется двуциклический способ расщепления по направлениям [3]. Представим данную схему для искомой функции й на произвольной сетке (г,^), которая может быть основной (п,ш) либо одной из вспомогательных (п',ш), (п, ш'):

3 г3_ + кН | ---3 | = 0,

т/2

2

~к+1/2 ~к+1/4 / ~к+1/2 ~к+1/4~

гз гз _+ К™

т/2

2

~к+3/4 ~к+1/2 / ~к+3/4 . ~к+1/2~ --+ К2

т/2

2

Ч 3 - + кН —-7Г-4- I = 0,

т/2

2

(12)

К^, К2Н представляют собой конечно-разностную аппроксимацию конвективных операторов (11). Во внутренних точках сеток рассматривается аппроксимация, использующая центральные разности вида

кк к / '*+1 — 1 + / / 1 г+1Ч '*+1,3 — V 1 г-1 ,31 ,3

I ^ 7 I +

2ЛХ

2ЛХ

1

(^1к+1/2,з Шг+1,3 — к—1/2,3 Ш*-1,з)

/2

2 Л

к ( Шг,3+1 — Шг,3-1

(13)

*3

2Л

+

кк ^2к,3+1Шг,3+1 — Щ.з-'з-1

2Л„

1

2Лу

кк ^2г,3+1/2Шг,3+1 — -1/2 Шг,3 -1

(14)

Здесь значения компонент скорости с предыдущего к-го слоя в средних точках рассчитываются в зависимости от искомой функции. Для температуры на этапе конвекции в формулах (12)-(14) искомой функцией Ш является Т, а для расчета коэффициентов используются формулы

V1 к+1 / 2,3 = 2(и4 + и1к+1,з)' ^2-3+1/2 = 2(м2кз + и2-з+1)-

2

При этом учитывается, что диффузионная скорость (и^-, и2к3-) уже заранее рассчитана на основной сетке (п, т). Если искомой функцией является первая компонента конвективной скорости /1, рассчитываемая на сетке (п, т'), то при прогонке по х, например, имеет место следующий пересчет коэффициентов в (13) в средней точке по значениям функции тока —, известным на основной сетке (п,т):

/1га+1/2,т' = 2

1 ( —п+1 —п+1

+

(15)

а при прогонке по у для первой компоненты конвективной скорости /1 полагаем уже в (14) к к

/2га,т/+1/2 = и2га,т+1- (16) Для искомой второй компоненты скорости имеют место соотношения, аналогичные (16) при прогонке по х и (15) при прогонке по у.

0

0

В граничных точках используется аппроксимация направленными разностями:

т^Ъ ~ 1 к ~ т^Ъ ~ 1 к ~

= 1+1/2,3^, К1 = - -1/2,,",3 ,

' ¿X 'Ж

т^Ъ ~ 1 к ~ т^Ъ ~ 1 к -

К2 ^«¿,1 = Т ^2^,1+1/2^^^,2) К2 = - Т"1/2^г,М •

'у Ту

Заметим, что для компонент конвективной скорости на смещенных сетках учитываются условия прилипания при выходе за границу вспомогательной сетки и попадании на границу основной сетки. Аппроксимация конвективных операторов в виде (13), (14) обеспечивает кососимметричность разностных операторов конвекции. Полученные системы разностных уравнений (12) решаются с помощью метода прогонки. Хотя для матрицы системы и не выполнены условия диагонального преобладания, метод прогонки устойчив [4]. Вышеизложенные аппроксимации операторов К1, К2 имеют второй порядок [3].

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

При переходе на этап диффузии определяем вихрь конвективной скорости ) на основной сетке

_ ^2га+1/2,т — ^2га-1/2,т ^1га,т+1/2 — ^1га,т-1/2

(«'га,т 7 , )

и именно он берется в качестве начального условия для этапа диффузии ()к = )). При этом граничные условия для конвективного вихря могут быть вычислены, исходя из условия переноса вихря, рассмотренного на границе, с использованием условий прилипания и вида операторов К1, К2, записанных с учетом того, что на нижнем временном слое = 0. Например, при х = 0:

д) 1

Ж + 2о 5Х = 0

или

кк — + 1 + °1,т ^12,т — = о

т 2 2 Тх = •

Этап диффузии реализуется с использованием прогонки с параметрами. Для начала представим разностную схему второго порядка аппроксимации для решения уравнения (9) в виде [4] (крышку над искомыми функциями опускаем):

Т1Т2)П+1 = ^кт, п = 2, ••., NN, т =2, М,

где

Т = Е — 2Т2 = Е — 2= (в + 2(В + 2^ + т/к*1^

Решение системы будем искать в виде

о к+1 _

^П^ =)п,т + ^п^М+т + ат ^Пд1 + Дг^П+М, (17)

п = 2, •••, N4, т = 2, •••, М•

При этом а, в находятся в результате решения задач:

Т1ап = 0, а1 = 1, = 0, Т вп = 0, в1 = 0, вм = 1

Аналогично предыдущему, но с оператором Т2 вычисляются а, /3. Для определения и имеем задачу

о к+1 к

Т1Т2 ипт = (18)

с нулевыми граничными условиями.

С целью нахождения соответствующей функции тока — построим итерационный процесс для решения уравнения

д- = А(Д- + и)

на основе продольно-поперечной конечно-разностной схемы, известной как метод переменных направлений [6] и имеющей второй порядок аппроксимации:

—S+1/2 - —s Т/2

— S+1 — —s+1/2

т/2

Л Л

Л1^+1/2 + Л2—S + C +arawS+1/2 + вп^+т/2 + 1 + Arn^S M

Л1 —S+1/2 + Л2^+Ч C +a„cS+1/2 + encN+m/2 + атЧЙ1 + вт<M

(19)

Здесь Л — итерационный параметр, а вихрь на границе связан с функцией тока условиями типа Тома [2], например, при n = 1:

<1/2 = - </2- (20)

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

Считается, что для диффузионной скорости выполнены условия прилипания на границе области и, как следствие, на границе полагается — = 0.

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

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

1. Переход на новый временной слой tk+1 начинается с расчета температуры Tk+1 на основной сетке по схеме (12) (w = T) для этапа конвекции и по схеме (18) с нулевой правой частью и соответствующими граничными условиями для этапа диффузии.

2. При найденном T осуществляется расчет конвективной скорости у на вспомогательных сетках и вихря скорости с — на основной сетке.

3. Затем определяется компонента w диффузионного вихря cCk+1 из (18). На каждом k-м временном слое расчета с вводится внутренний итерационный процесс расчета —^ из системы уравнений (19). После окончания итераций при s = S считается, что с заданной точностью е^ определены значения — на (k+1)-временном слое: —k+1 = —S. Следовательно, с использованием (17) будет определен и вихрь W.

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

Остается заметить, что итерационный процесс для — считается сошедшимся, если выполнен критерий сходимости [2, 8]

max n+m -—nm <^max n+m1,

где s — номер итерации; е^ — заданная точность расчета —s+1.

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

В качестве основного теста рассматривается задача о свободной конвекции в квадратной полости 0<х< 1, 0<у< 1 при подогреве сбоку [8,9] (Кв = Рг = 1, 0г = 104).

Тест 1. В начальный момент времени жидкость покоится и нагрета по закону Т(х,у, 0) = х. Граничные условия соответствуют твердым непроницаемым границам с заданной температурой (Т(х,у,Ь) = х). Тест заключается в получении стационарного решения до удовлетворения неравенств [4]

|Шк+1 - | < е,

где число Нуссельта Ки определяет интегральный поток тепла на правой границе х = 1. (Число Нуссельта рассчитывается с использованием аппроксимации второго порядка для первой производной на границе.) При тестировании полагается:

е = 10-12, бф = 10-5, т = то = к2/4, Н = Нх = Ну.

Результаты расчетов: характеристика интенсивности течения ^тах и интегральный тепловой поток на правой стенке х = 1 при Ог = 104 — приводятся в табл. 1, где в последней строке даны предельные значения характеристик из [8]. На рис. 1 представлены типичные линии тока и изотермы при Ог = 104 для первого теста.

Рис. 1. Линии тока (а) и изотермы (б). Тест 1. _Таблица 1

Сетка фтах Ми

9 х 9 7.3558 1.9405

11 х 11 7.0415 1.8886

17 х 17 6.6450 1.8052

21 х 21 6.5460 1.7831

33 х 33 6.4360 1.7607

[7] 6.37 [8] 1.752 [8]

Тест 2. В начальный момент времени жидкость покоится и Т(х, у, 0) = 0 или х. Граничные условия для температуры: Т = 0 при х = 0, Т = 1 при х = 1, дТ/ду = 0 при у = 0, 1.

Результаты расчетов при Ог = 104 приводятся в табл. 2. Для проверки работы алгоритма в условиях переходного режима для второго теста был продолжен расчет до повторного выхода на стационар с граничным условием Т = 2 при х = 1, достигаемым линейным по времени нагревом в течение промежутка времени Ъ = 20т. При этом на сетке 21 х 21 получены значения ^>тах = 6-8024, Ми = 6-0907. На рис. 2 и 3 приведены типичные линии

Таблица 2

Сетка ^тах Ш

9 х 9 6.0645 2.7790

11 х 11 5.7397 2.6657

17 х 17 5.3438 2.4558

21 х 21 5.2474 2.3920

[8] 5.072 2.243

Рис. 2. Линии тока (а) и изотермы (б). Тест 2.

Рис. 3. Линии тока (а) и изотермы (б) после повторного выхода на стационарный режим.

тока и изотермы при Ог = 104 для второго теста. На рис. 3 представлены линии тока и изотермы установившегося течения после повторного выхода на стационарный режим.

Устойчивость алгоритма и порядок сходимости (близкий ко второму) проверяются путем вычислительного эксперимента на последовательностях сеток аналогично [5]. Ограничения на шаг т не обнаружены. Подробные теоретические исследования аналогичного этапа диффузии проведены в [4].

Хорошее согласование с тестами показали также расчеты при других значениях диапазона изменения числа Грасгофа (Ог = 1.4 • 104; 8 • 104; 2 • 105.) При этом подтверждается, например, наличие двух вихрей внутри полости и сосредоточение основного течения вблизи стенок при Ог = 2 • 105 (см. [5, 8]).

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

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

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

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

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

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

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

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

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

[9] Математическое моделирование конвективного тепломассообмена на основе уравнений Навье — Стокса / В. И. Полежаев, А. В. Бунэ, Н. А. Верезуб и др. М.: Наука, 1987. 220 с.

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