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

УДК 532.516

О.Н. Гончарова

Метод расщепления по физическим процессам для расчета трехмерных задач конвекции

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

1. Введение. Основные уравнения.

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

Уг + V -Чр'+^т Др+у; (1)

Не

сііу V = 0;

Ті + V -УТ =

1

-АТ.

(2)

(3)

ЯвРт

Будем исходить из заданного при і = 0 начального состояния V = 0, Т = Т0(ж) (ж Є О). Граничные условия для V на £ соответствуют условиям прилипания: V = 0, (ж Є £Є [0Для температуры может быть рассмотрено одно из условий первого, второго или третьего рода:

дТ

т = Тв(х,г), — = Тв(х,г);

дп

дТ

— + а(Т - Тв) = 0; ж є£,і Є [0Л]. дп

Здесь р' - модифицированное давление (отклонение от гидростатического давления); р =

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

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

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

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

векторному потенциалу скорости. Поскольку вычисляемыми функциями являются как компоненты скорости, так и новые искомые функции, то принципиальным моментом в организации расчета является введение смещенных сеток (см., например, [10]). Смещенные сетки вида (п',т,1), (п,т',1), (п,т,1г) вводятся для расчета соответствующих компонент векторного потенциала и ротора скорости. Сетки (и, т1,1'), (п1, т, V) и (и', т', I) будут использоваться для расчета компонент конвективной скорости, а (и, т, I) - для расчета температуры. Смещенные сетки позволяют проследить за автоматическим выполнением условия несжимаемости поля скоростей и организовать восстановление скорости через векторный потенциал. Для реализации этапа диффузии применяется метод стабилизирующей поправки [1, 3]. Для новых искомых функций на твердых непроницаемых границах используются условия, выражающие равенство нулю касательных составляющих векторного потенциала и производной по нормали его нормальной составляющей [11-13]. Тестирование метода расщепления по физическим процессам в трехмерных задачах проводится на известных задачах о свободной конвекции в кубической полости или параллелепипеде при подогреве одной грани [13-15].

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

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

вективного движения жидкости как последовательную реализацию двух вышеназванных этапов. Перепишем уравнения движения (1) в виде

V-V

KV = 0:

V-V

= -Vp' + dV + f.

(4)

(5)

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

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

W — W

= DW + rot F.

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

ДФ = -W.

(6)

(7)

Здесь К - оператор конвективного переноса К = V ■ V, В оператор диффузионного переноса В = РД, V = Не-1. При вычислении температуры па этапе диффузии полагается V = (КеРг)-1.

Кроме того, V = пйФ, V = пйУ, V = гсЛУ, и выполнено сИуФ = 0. Опустим ’’крышку” при обозначении искомых функций на этапе диффузии и представим искомые функции в компонентах:

V = (т1,т2, т3) =

' дщ дщ дщ дщ дщ дщ\

ду дг1 дг дЖ дж ду

V = (щ,щ,щ) =

(дф дгф2 д'фх дф дгф2 дф\ \

\ ду дг , дг дж , дж ду ) ,

Ф = (ф1,ф2,фз)-

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

дф

ж = 0,ж = 1 : —— = 0,т = 0, ф = 0, дж

w2 = —■

д2ф2

дж2

,фз=0,ыз = -

д2фз_ дж2 ’

дф2

у = 0,у=1:-^. =о,т2= о, ф3 = 0, ду

дф , п д2 ф\

т ~-----ду2 , ф1 _ ^,т ~--ду2 ’

дф

г = 0,г = 1 : —— = 0,т = 0, ф = 0, дг

дф , п д2ф2 /оЛ

т = — -о^-,ф2=0,т = —^т- (8)

3. Реализация этапа конвекции. При построении разностной аппроксимации оператор

конвективного переноса К из (4) будем рассматривать в виде суммы одномерных кососимметричных операторов, результат действия которых па искомую функцию / (одну из компонент условной конвективной скорости или температуру) следующий: К/ = К/ + К2/ + К/ [3, 7, 16], где

і/ 2 д^ дх

К/ = -^ 2

К/=-2

1 ( д/ , дщ/

ду

/

дх

ду

дщ/

дх

(9)

фз

ф2

Ну

¡+1 - ф2и,т',1

Нг

ф1 п,т,1+1 ф1 и' ,т,1

Нх

щ3 и',т,1

ф3п+1,т,1' ф3и,т,1' _

Нх ’

ф2 п+1,т' ,1 ф2 и,т,1

Нх

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

Начальные условия для компонент конвективной скорости V = {щ,щ,Ъ3) задаются с предыдущего временного слоя и представляют собой соответствующую компоненту диффузионной скорости, рассчитанную на нужной смещенной сетке, а именно: на сетке {и, т', I') - для первой компоненты, па сетке (п', т, ¡') - для второй компоненты и на сетке (п', т', I) - для третьей компоненты. ’’Крышка” для обозначения диффузионной скорости опускается и используются следующие разностные представления компонент диффузионной скорости через соответствующие компоненты векторного потенциала скорости:

/*+* - /к т/2

/к+| - /к+1 т/

/к+| - /к+| Т/

/к+5 - /к+|

Т/

/ к+| - / к+! Т/

/к+1 ___ /к+1

Т/

К*

К*

к+її

/к

= 0:

К

и /к+~в +/к+г

К

и /к+її +/к+?

= 0:

= 0;

= 0:

= 0:

К

к

-)=0. (10)

ф1 п',т+1,1 ф1 и' ,т,1

Ну

Заметим, что при использовании указанных конечно-разностных аппроксимаций тождественно выполняется сИу V = 0 на сетке

(п', т', V).

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

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

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

4. Реализация этапа диффузии. Для реализации этапа диффузии (6)-(7) используются смещенные сетки (п', т, I) для расчета первых компонент ротора скорости т и векторного потенциала -ф, (п, т',1) - для расчета вторых компонент т, ф и (п, т, ¡') - для расчета третьих компонент т, фз-

Заметим, что при переходе на этап диффузии ротор конвективной скорости рассчитывается как раз на нужных вспомогательных сетках, и имеют место следующие представления (’’волна” опускается):

«3 и' ’ ,т' ,1 «3 и' ;т'

,1 Ну

V и ' ,т,1' ' «2 и' ,т,1 —1

«1 и, Нг ,т' ,1 «1 и,т',

,1 — «3 и' ' ,т' ,1 Нг ! «3 и'—1 ,т',1

«2 и' Нх ’,ш,1' «2 и' —1

' «1 и ,т' ,1' Нх «1 и,т' — 1,1'

Ну

-ц

Равенство сііуТУ = 0 тождественно выполняется на сетке (и, ш,1).

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

этапа диффузии Жк = Ж. Продемонстрируем метод стабилизирующей поправки [1,3] на примере первой компоненты ротора скорости и векторного потенциала (индекс 1 для удобства изложения опустим). Представим разностную схему второго порядка аппроксимации для решения уравнения (6) в следующем виде:

,,,к+1/3 _...к

-------------= Ліадк+1/3 + А2тк+

,к+2/3 _ ,,,к+1/3

,,к+1 _ тк+2/3

= л2(.^2/3

— А3 (■

тк+1 -тк

(П)

т = ———— из (о) используется аппроксимация д'Г

следующего вида:

2 і 1 і

ти', 1,1 — т^'Фи',2,1 — 77 Фп',3,1-

Ч Ч

Эф

(12)

Здесь введены обозначения Л* (г = 1^2, 3) для операторов второй разностной производной по пространственной переменной, умноженных на -/к -д2 ^ -д2 ^ - д,

V (Л1 ~ ~ "д^1 3 ~ удгУ'а правая

часть представляет собой конечно-разностную аппроксимацию компоненты (го!;^1)!.

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

Реализация граничных условий = 0, т = О,

дх

заданных при х = 0 и х = хо, осуществляется на смещенных границах с помощью представлений ф0,т,1 = фх,т,1 И то,т,1 = -т1 ,т,1 ДЛЯ ’’левой” смещенной границы и аналогичных для ”правой”.

Соответствующая компонента векторного

ф

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

на основе схемы (11). В операторах Л* (см. схему (11)) нужно использовать вместо V коэффициент Л, являющийся итерационным параметром. Заметим, что итерационный процесс вводится на каждом временном слое к.

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

у, г, х г, х, у

для третьих. Поскольку аппроксимации граничных условий вида (12) связывают значения компонент ротора скорости и векторного потенциала, то проводится дополнительная проверка сходимости итераций (см., например, [8, 18]).

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

1. Переход на новый временной слой Ьк+1 начинается с расчета температуры Тк+1 па основной сетке (и,т,1) по схеме (10) для этапа конвекции и по схеме (11) с нулевой правой частью и соответствующими температурными граничными условиями на этапе диффузии.

Тк

чет компонент конвективной скорости Vк+1 по схеме (10) и компонент ротора конвективной скорости на вспомогательных сетках.

к

3. Вычисляется диффузионный вихрь (ротор диффузионной скорости) по схеме (11). На

к

дится внутренний итерационный процесс расчета компонент векторного потенциала по схеме, аналогичной (11), с соответствующими граничными условиями. После окончания итераций при в = S считается, что с заданной точностью е^ определены значения компонент векторного по-к

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

Остается заметить, что итерационный проф

ала) считается сошедшимся, если выполнен критерий сходимости, например, следующего вида [11,18]:

max; - Ф-

nml \ <єФ 'тах; 11\>

где в - номер итерации, е^ - заданная точность расчета ф5+1. Помимо этого, осуществляется проверка малости величины

є = max \^n',і , і(Фк+1)

n',l

, 1, l

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

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

хх

дТ

14]: Т = 0 при х = 0; Т = 1 при х = х0, -г— = 0

дп

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

Nunv —

1

У 0*0

zo

lNumean( z)dz,

где

Уо

NUm„A z>=JNul°c{ y,*dy,

. . дТ,

Нщ°с{ у,г) = дх *=*

и х = 0 или х = х0. Аппроксимации второго порядка для первой производной температуры на границе используются при вычислении Мщ°с.

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

Сила тяжести направлена против оси Оу, а область Л представляет собой куб, определяемый значениями щ = 1 ,у = 1,г = 1 • При тестировании полагается х = х0, От = 105, Рт = 1, Яв = 1. Расчеты ведутся па сетке 41 х 41 х 41, т = 0.0001, а условием выхода на стационарное решение является выполнение неравенства

\Мпк+г — N4^ \ <еИи,

где £ми = 1 о-, £.ф = 10-4, Є = 10-4. Сравнение расчетов с известными результатами представлено в следующей таблице.

Таблица 1

цитир. сетка Nuov

[13] 86 х 65 х 65 4.34

[15] 62 х 62 х 62 4.36

наст, работа 41 х 41 х 41 4.38

Получено также, что изолинии температуры в плоскости я = zo/2 представляют собой типичную картину для поля температуры в двумерном случае.

Тест 2 ([14]).

Сила тяжести направлена по оси Oz. Область течения представляет собой куб, определяемый размерами щ = 1, уо = 1, zo = 1. При тестировании полагается x =

0, Ra = 5 • 105 (чисто Рэлея), Pr = 0.71, Re = 1/Pr.

Таблица 2 представляет сравнение результатов с расчетами, проведенными в [14], где отмечено, что большинство расчетов проводились на сетке 15x15x15, а некоторые - на сетке 51 x 51 x 51. Для сравнения берется не только чисто Нуссельта, но и максимальное значение величины ф2 при у = у0/2 (ф

max)-

Таблица 2

цитир. Nuov Фтах

[14] 7.37 18.4

наст, работа 7.30 18.78

Тест 3 ([14]).

Как и в предыдущем случае, сила тяже-

Ог

чения представляет собой прямоугольный

х,

уо = 2, го = 1. Расчетная сетка остается хх

х = 0, Яа = 1.5 • 105, Рт = 1, Яв = 1.

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

n

Таблица 3

цитир. Nuov Фтах

[14] 5.37 13.2

наст, работа 5.01 13.12

Отметим, что этап диффузии представляет возможность расщепления по направлению, ко-

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

Литература

1. Яненко П.П. Метод дробных шагов решения многомерных задач математической физики.

- Новосибирск, 1967.

2. Самарский A.A. О принципе аддитивности для построения экономичных разностных схем // ДАН СССР. - 1965. - Т. 165. - ,V'6.

3. Марчук Г.11. Методы расщепления. - М., 1988.

4. Белоцерковский О.М. Метод расщепления в

применении к решению задач динамики вязкой несжимаемой жидкости / О.М. Белоцерковский, В.А. Гущин, В.В. Щенников // Журнал вычислит, математики и матем. физики. - 1975. - Т. 15. .Y'l.

5. Белоцерковский О.М. Моделирование некоторых течений вязкой жидкости / О.М. Белоцерковский, В.А. Гущин. - М., 1982.

6. Белоцерковский О.М. Численное моделирование в механике сплошных сред. — М., 1984.

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

8. Воеводин А.Ф. Метод расчёта вязких течений в замкнутых областях / А.Ф. Воеводин, Т.В. Протопопова // Сиб. журн. индустриальной математики. - 2001. - Т. 4. - .V" 1.

9. Воеводин А.Ф. Метод расчета двумерных задач конвекции на основе расщепления по физическим процессам / А.Ф. Воеводин, О.Н. Гончарова // Вычислительные технологии. - 2002. - Т. 7. - т.

10. Белолипецкий В.М. Математическое моделирование течений стратифицированной жидкости / В.М. Белолипецкий, В.Ю. Ко-стюк, К).II. Шокин. - Новосибирск, 1991.

11. Роуч П. Вычислительная гидродинамика. — М., 1975.

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

- Т. 41. - ,V'9.

13. Бессонов O.A. Тест для численных решений трехмерной задачи о естественной конвекции в кубической полости / O.A. Бессонов, В.А. Брайловская, С.А. Никитин, В.И. Полежаев // Математическое моделирование. — 1999. - Т. 11. .Y'l2.

14. 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.

15. Мызникова Б.И. Процессы установления стационарных конвективных течений в кубической полости при подогреве снизу / Б.И. Мызникова, Е.Л. Тарунин // Нестационарные процессы в жидкостях и твердых телах. - Свердловск, 1983.

16. Самарский A.A. Разностные схемы с операторными множителями / A.A. Самарский, П.П. Вабищевич, П.П. Матус. - Минск, 1998.

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

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