CC BY
74
УДК 519.6+532.5
А.Ф. Воеводин, О.Н. Гончарова, Т.В. Протопопова
Численные методы исследования конвективных течений: реализация метода расщепления по физическим процессам*
A.F. Voevodin, O.N. Goncharova, T.V. Protopopova
Numerical Methods of Investigation of the Convective Flows: Realization of the Splitting Method into Physical Processes
Для расчета конвективных течений предлагается численный метод, в котором реализована идея расщепления по физическим процессам. Предлагаемая схема гарантирует энергетическую нейтральность поля скоростей, обладает свойством устойчивости в линейном приближении. Расщепление на конвективный и диффузионный переносы проводится в уравнениях конвекции Обербека-Буссинеска. В работе описаны особенности реализации этапов, приведены тестовые задачи и результаты тестирования численных алгоритмов.
Ключевые слова: конвективные течения, конвекция, диффузия, методы расщепления, расщепление по физическим процессам, численный алгоритм.
It is proposed to use a numerical method based on realization of splitting into physical processes for numerical investigations of the convective flows. The scheme guarantees an energetic neutrality of the velocity field and has a property of stability in linear approximation. Splitting into convective and diffusive transfer is realized in the equations of convection. The paper describes the features of realization of the stages, test problems and results of testing of the numerical algorithms.
Key words: convective flows, convection, diffusion, splitting methods, splitting into physical processes, numerical algorithm.
Введение. Методы численного решения уравнений движения вязкой несжимаемой жидкости можно условно разделить на две группы. Первую составляют методы в естественных переменных «скорость-давление». Основное затруднение при таком подходе заключается в определении граничного условия для давления. Ко второй группе относятся методы, связанные с введением вихря скорости и функции тока (в двумерных областях) или ротора скорости и векторного потенциала (в трехмерных). Достоинство такого подхода в том, что нет необходимости заботиться о соленоидальности вектора скорости (условие выполняется автоматически). Это особенно важно в случае задач конвекции, когда нет притока или оттока жидкости в область течения. Однако при таком подходе возникают трудности, связанные с заданием и численной реализацией граничного условия для вихря (ротора скорости) на стенке с прилипанием. Кроме того, в трехмерном случае (в отличие от двумерного) в уравнениях переноса компонент ротора скорости появляются слагаемые, описывающие усиление
*Работа выполнена в рамках проекта № 7.3975.2011 Алтайского государственного университета (поддержан Министерством образования и науки РФ).
вихря за счет растяжения вихревых нитей. Наличие этих членов вносит дополнительные трудности (по сравнению с двумерным случаем). Для расчета конвективных движений жидкости в двумерных (прямоугольники) и трехмерных (прямоугольные параллелепипеды) областях авторами предлагается численный метод, в котором реализована идея расщепления по физическим процессам. Предлагаемая схема расщепления является физически оправданной (гарантируется энергетическая нейтральность поля скоростей) и обладает свойством устойчивости в линейном приближении.
1. Расщепление по физическим переменным. Этап конвекции. Расщепление на конвективный и диффузионный переносы в двумерном случае может быть реализовано для уравнений конвекции Обербека-Буссинеска [1, 2], записанных как в переменных «вихрь - функция тока», так и в физических переменных. В трехмерном случае расщепление проводится в физических переменных. Конвективный и диффузионый переносы естественным образом выделяются как в уравнениях движения, так и в уравнении переноса тепла. Идея метода расщепления по физиче-
ским процессам базируется на методе слабой аппроксимации и аддитивности этих процессов для достаточно малых шагов по времени [3, 4]. Общая теория расщепления наиболее полно изложена в
[5] и развивается в [6-8].
1.1. Особенности реализации этапа конвекции в двумерном случае (переменные «вихрь — функция тока»). Выделение этапа конвекции в уравнении переноса вихря позволяет получить на этапе диффузии для вихря линейную систему Стокса. Реализуется рассматриваемый этап на основе схем Кранка-Николсона, обладающих свойством энергетической нейтральности. Постановка граничных условий не требуется в силу гиперболичности уравнения и условий равенства нулю компонент скорости на границах области, занятой жидкостью.
1.2. Особенности реализации этапа конвекции (физические переменные). Выделение этапа конвекции позволяет исключить расчет градиента давления, обеспечивает корректность расщепления для выполнения граничных условий. Из условий прилипания и гиперболичности системы уравнений следует, что граничные условия на этапе конвекции являются следствием уравнений. В этом случае этап конвекции реализуется для компонент условной скорости (вспомогательной функции) также на основе элементарных схем Кранка-Николсона. Принципиальным моментом в организации расчета является введение смещенных сеток (см.: [9]). Сетки (п,т',1'), (и',т,1') и () используются для расчета компонент вспомогательной конвективной скорости в случае трехмерных задач, а (п, т, I) - для расчета температуры.
2. Этап диффузии. Если расщепление по физическим процессам проведено для системы уравнений конвекции, записанной в физических переменных, то на этапе диффузии осуществляется переход к переменным «вихрь - функция тока» для двумерных задач и «ротор скорости - векторный потенциал скорости» для трехмерных задач. При этом в двумерном случае получается система Стокса, а в трехмерном случае система векторных уравнений распадается на три независимых системы Стокса для каждой пары: «компонента ротора скорости - соответствующая компонента векторного потенциала». Основная трудность построения численных методов для решения таких задач заключается в том, что соответствующие разностные уравнения оказываются перевязанными условиями на неподвижных твердых стенках. Для преодоления этой трудности разработан оригинальный двухполевой безытерационный способ расчета (прогонка с параметрами) с точным выполнением (в разностном смысле) следствий условий прилипания.
2.1. Одномерная линейная задача на отрезке 0 < х < 1 для иллюстрации прогонки с параметрами. Пусть
дш д2ш г
дЬ дх2 + ’
д2ф
я 2 + Ш = 0 дх2
(1)
(2)
где /(х) - заданная функция.
Уравнения (1), (2) можно рассматривать как одномерный аналог системы уравнений Стокса. Для полной постановки задачи к уравнениям нужно добавить начальные и граничные условия. Рассмотрим их в виде:
ф(х, 0) = фо(х), ш(х, 0) = ш0(х), (3)
ф(0,Ь) = ф(1,Ь) = 0, (4)
дФ (о.о=дф (і.<)=о. (5)
Такой выбор граничных условий моделирует условия прилипания и непротекания в двумерном случае.
На отрезке [0,1] введем равномерную сетку с шагом Н =1/1 : Ои = [х% = іН; і = 0,..., І}, а по переменной Ь - сетку ОТ = [Ьп = пт; п = 0,1,... }. На сетке О и х ОТ аппроксимируем уравнение (1) следующей неявной разностной схемой:
(Е - тК)шП+ = шП + т/П+1/2, (і = 1,.., І - 1),
(6)
где Л = д2/дх2 + 0(Н2) - оператор второй раз-
хП+1/2
ностной производной по направлению х, / =
(/П + /П+1)/2.
Систему (6) в стандартном виде можно записать как систему разностных уравнений второго порядка
ЬшП+1 = -аш—1 + Ь^1 - аи^1 = и,
I = < + т/П+1/2, (г = 1,...,1 - 1). (7)
При этом выполняются следующие условия на коэффициенты разностных уравнений:
а = —77 > 0, Ь =1 + 2а. Н2
(8)
Уравнения (7) не замкнуты, и их решение зависит от двух параметров. Обозначим эти параметры ш1, ш2, тогда решение этой системы можно представить в виде
шП^1 = Ш0П+1 + а4ш1 + вгш2, (г = 0,. .. ,1). (9)
Если выбрать в качестве параметров значения Ш0П+1 и шП+1, то неизвестные ш0П+1, а*, в* в силу линейности системы (7) определяются путем решения следующих задач:
Ьш0п+1 = /і , ш0п+1 = ш0п+1 = 0; (10)
Ьаі =0, ао = 1, аі = 0;
Щ = 0, во = 0, ві = 1.
(11)
(12)
Системы уравнений (10)—(12) при условиях (8) эффективно решаются методом прогонки. Следует учесть, что одномерные векторы-столбцы а, в не зависят от индекса п и определяются вне цикла по времени.
Для аппроксимации уравнения (2) используем следующую разностную схему:
ЛФП+1 + шП+1 =0, (г = 1,...,1 - 1). (13)
Система уравнений замыкается с помощью граничных условий
ФП+ = ФП+ = 0. (14)
Подставим представление (9) в уравнение (13):
АфП+1 + Ш0П+1 + щт1 + ві.^2 = 0.
(15)
Для исключения неизвестных параметров в (15) используем условия (5). Аппроксимируем эти условия, например, со вторым порядком точности. Следствиями такой аппроксимации являются условия Тома для вихря скорости на твердой стенке [2]:
ФП+1 - Ф,
п+1
+ Н-
ФП+1 - ФП+1
.п+1
.п+1
+ Н-
2
0(Н2),
0(Ь2). (16)
Если в представлении (9) в качестве параметров выбраны значения шП+1, шП+1, то с помощью (16) они легко исключаются. Для фП+1 получаем замкнутую систему уравнений
Лфп +1 = -ш0п+1 + 2а.
п +1
фп+1
1 Н2
первого и последнего столбцов и для ее обращения требуются специальные алгоритмы, например, метод прогонки с параметрами.
Расчет одного шага по времени завершается вычислением неизвестных шгП+1 по формулам (9) с использованием условий (16).
Отметим, что в данном алгоритме уравнения для вихря и функции тока решаются отдельно (двухполевой метод), оба граничных условия, по сути, используются при решении уравнения для функции тока, при этом точным образом (в разностном смысле) реализуются граничные условия для вихря на твердых стенках. Алгоритм, подробно изложенный для одномерной модельной задачи, нашел свое обобщение для двумерных и трехмерных задач. Важным моментом явился тот факт, что уравнение для функции вихря (или каждой компоненты ротора скорости) можно аппроксимировать факторизованными разностными схемами с перестановочными операторами.
2.2. Особенности этапа диффузии для двумерных задач в переменных «вихрь — функция тока». Для двумерных задач алгоритм был обобщен на случай использования схем повышенного порядка аппроксимации: применялись схемы четвертого порядка аппроксимации внутри области и аппроксимация граничных условий с третьим порядком по пространственной переменной (условия Вудса). В двумерном случае методом энергетических неравенств доказана устойчивость разностной краевой задачи [10].
Особенности реализации этапа диффузии для трехмерных задач. Смещенные сетки (п',т,1), (п,т',1), (п,т,1') вводятся для расчета соответствующих компонент векторного потенциала и ротора скорости на этапе диффузии. Так называемый конвективный вихрь берется в качестве начального условия для этапа диффузии. При написании граничных условий для новых искомых функций следует исходить из их выполнения с учетом смещения некоторых границ. Для реализации этапа диффузии строится алгоритм прогонки с параметрами [7, 8], который может быть назван обобщением на трехмерный случай алгоритма, изложенного в разделе 2.1.
Алгоритм прогонки с параметрами для первых компонент ротора скорости и векторного потенциала т1, ф1 (индекс «1» для удобства изложения опустим) предполагает, что имеет место следующее соотношение:
„ . , = „,000 „п' ,т,1 „п ,т,1
+ а2т„п ',0,1 + вт „п' ,М,1 +
0
0
2
і
ФП+1 = ФП +1 = 0. (17)
Решив полученную систему уравнений, найдем ФП+1. Матрица системы (17) отличается от трехдиагональной матрицы наличием ненулевых
+а1 тп ^,т,0 + в1 тп' ,т,Ь. (18)
Первое слагаемое отвечает за удовлетворение однородным граничным условиям для компоненты ротора скорости на двух гранях па-
раллелепипеда, а граничные условия на четырех оставшихся гранях называются параметрами. Коэффициенты при параметрах вычисляются только один раз в результате решения конечноразностных задач (возможность изложенного нахождения аі, ві обеспечивается перестановочностью соответствующих разностных операторов). Пользуясь представлением (18), связью компонент ротора скорости и векторного потенциала, перестановочностью разностных операторов, получаем конечно-разностную схему для вычисления искомых функций. Конечно-разностная схема второго порядка для вычисления слагаемого „°°<3т і предусматривает последовательное осуществление расщепления по направлениям. Имеет место следующая очередность прогонок по направлениям для первых компонент ротора скорости и векторного потенциала: прогонка по г, прогонка по у, прогонка по х. Отметим, что установлена принципиальность очередности расчетов по направлениям, решена проблема постановки граничных условий для вспомогательных функций на внутренних этапах конечно-разностной схемы. Алгоритм расчета реализован таким образом, что граничные условия для вспомогательных функций на внутренних шагах схем являются однородными. Оригинальность алгоритма расчета состоит также в том, что компоненту ф векторного потенциала не вычислить без первого слагаемого „000 разложения (18), а компоненту „ ротора скорости не определить без компоненты векторного потенциала ф.
Общая схема решения задачи состоит в осуществлении следующих этапов.
1. Переход на новый временной слой Ьк+1 начинается с расчета температуры на основной сетке для этапа конвекции и с помощью, например, метода стабилизирующей поправки [5] на этапе диффузии.
2. При найденном распределении температуры осуществляется расчет компонент конвективной скорости и компонент ротора конвективной скорости на вспомогательных сетках.
3. Вычисляются составляющие компоненты диффузионного вихря (первое слагаемое в представлении (18)) с использованием конечноразностной схемы второго порядка. На каждом к-м временном слое расчета вводится внутренний итерационный процесс расчета компонент векторного потенциала также с помощью соответствующей конечно-разностной схемы второго порядка и с соответствующими граничными условиями. После окончания итераций при в = 6” считается, что с заданной точностью є определены значения компонент векторного потенциала на (к+1)-м временном слое. Окончательно, с использованием представления (18), будут вычислены и компоненты диффузионного вихря.
4. При переходе на следующий временной слой (к этапу 1) насчитывается искомая диффузионная скорость. Итерационный процесс для компоненты векторного потенциала считается сошедшимся, если выполнен критерий сходимости [2, 8, 11].
Метод расщепления по физическим процессам обобщается на области с произвольными, достаточно гладкими границами, а также на случай зависимости коэффициентов переноса от температуры.
Тестирование метода проводится на известных задачах о свободной конвекции в кубической полости и параллелепипеде при подогреве одной грани [12].
Пример 1. В качестве теста рассмотрена задача о конвективном течении жидкости в области, представляющей собой квадрат. Считается, что температура задана таким образом, что соответствующее точное решение для функции ф(х,у,Ь) имеет вид:
ф(х, у, ¿) = ехр-4^ 4 х2(1 - х)2у2(1 - у)2.
Такой способ тестирования называют методом пробной функции. Он позволяет исследовать устойчивость и сходимость численного метода на точном решении. Результаты расчетов приводятся в таблице 1. Здесь е^ = \\фе - фп||/||фе||, Фе и ф - точное и численное решения, соответственно, а под ||ф|| понимается тахг^фг^. Значение еш определяется аналогично.
Таблица 1
Сетка єш єФ
16 х 16 0.0042750 0.0008108
32 х 32 0.0005775 0.0001198
64 х 64 0.0000674 0.0000191
Результаты демонстрируют наличие 3-го порядка сходимости по пространству.
Пример 2. Рассматривается задача о конвективном течении в полости квадратного сечения (0 < х < 1,0 < у < 1) при по-
догреве сбоку. Граничные условия соответствуют твердым стенкам с заданной температурой (Т(х,у,Ь) = х, (х,у) € Г). Начальные условия таковы: ш(х,у, 0) = 0, ф(х,у, 0) = 0, Т(х,у, 0) = х. Число Прандтля полагалось равным 1, а число Грасгофа — 104 и 105. Тест заключается в получении стационарного решения до удовлетворения неравенств ¡Ми^1 - Ми^] < ет , (т = 1, 2), где N41, Ми2 - тепловые потоки (числа Нуссельта) соответственно на правой и нижней стенках полости [13]; т - шаг схемы по времени.
Для вычисления чисел Нуссельта использовались аппроксимация производной по пространству трехточечной асимметричной формулой, а также формула трапеций численного интегрирования.
В численном эксперименте полагалось є = 10~4, т = Н2/4, Н - шаг сетки по пространству (сетка квадратная). Здесь же приведены предельные при Н ^ 0 значения характеристик, полученные в работе [13].
Таблица 2
Сетка Фтах N41 N42
О 1 О
16 х 16 6.48128 1.82880 1.26850
24 х 24 6.42051 1.78145 1.24972
32 х 32 6.39910 1.76459 1.22980
[13] Н ^ 0 6.37 1.752 1.175
От = 105
16 х 16 15.17258 3.79884 2.02810
24 х 24 13.60772 3.69271 2.31230
32 х 32 13.58138 3.56760 2.26870
[13] Н ^ 0 13.1 3.43 2.18
Пусть вектор силы тяжести направлен против оси Оу, область П представляет собой куб (х0 = 1, уо = 1, г0 = 1). Принимаются следующие значения безразмерных параметров задачи От = 105 (число Грасгофа), Рт = 1 (число Прандтля), Кв = 1 (число Рейнольдса). Тепловые потоки вычисляются на стенке X = хо. Расчеты проводятся на трех сетках 21 х 21 х 21, 41 х 41 х 41 и 81 х 81 х 81 с шагом по времени (по фиктивному времени для стационарной задачи) т = 0.0001. Условием выхода на стационарный режим является выполнение неравенства [7, 14]
\Мик0+ - Мпкоу | <єКи,
где єКи = 10-5.
Сравнение расчетов с известными результатами представлено в таблице 3.
Таблица 3
Цитир. Сетка N иОУ
[12] 86 х 65 х 65 4.34
[11] 62 х 62 х 62 4.36
Наст. работа 81 х 81 х 81 4.6342
Результаты показывают, что ошибка по каждой из характеристик на сетке 32 х 32 не превышает 5%.
Пример 3. Метод расщепления для решения трехмерных задач конвекции проверяется на тесте о свободной конвекции в замкнутой кювете П = {[0, х0] х [0, у0] х [0, г0]} при подогреве одной грани х = хо и теплоизоляции остальных [11,
дТ
12]: Т = 0 при х = 0, Т =1 при х = х0, —— = 0
дп
при у = 0, у = у0, г = 0, г = 20. При этом будем следить за безразмерными тепловыми потоками на изотермической холодной или нагретой стенке, определяемыми следующими числами Нуссельта:
Ыип
1
у020
Ыи„
г(х)йх,
Уо
ь{г)
N41
ОС (У, z)dy,
. . дТ,
Ми1ос(у,г) = — \х=х
и х = 0 или х = хо. Аппроксимации второго порядка для производной температуры Т на границе используются при вычислении Ыщос. Ставится задача достижения стационарного решения.
На сетке 21 х 21 х 21 расчеты, демонстрирующие выход практически на один и тот же стационарный режим, проведены также с разными значениями т. Получены значения интегральной величины Ыиоу, равные 6.21, 6.57, 6.63 при следующих значениях временного шага 10~3, 10~4 и 10~5 соответственно и е^и = 10~4 • т.
Таблица 4
Сетка К
21 х 21 х 21 40.83
41 х 41 х 41 31.13
81 х 81 х 81 28.56
Устойчивость алгоритма и порядок сходимости вычислительного алгоритма проверяются путем вычислительных экспериментов на последовательности сеток 21 х 21 х 21, 41 х 41 х 41, 81 х 81 х 81. Наблюдения ведутся за интегральной величиной К (интегральной нормой скорости), представляющей собой аналог кинетической энергии. Результаты тестовых расчетов приведены в таблице 4. Экспериментальный порядок сходимости, близкий ко второму г « 1.92, и приближенно определяемая относительная погрешность 3.24% вычисляются согласно правилу Рунге [15].
Библиографический список
1. Джозеф Д. Устойчивость движений жидко- 3. Яненко Н.Н. Метод дробных шагов реше-
сти: пер. с англ. — М., 1981. ния многомерных задач математической физики.
2. Роуч П. Вычислительная гидродинамика. — Новосибирск 1967.
— М., 1975. 4. Самарский А.А. О принципе аддитивности
для построения экономичных разностных схем // ДАН СССР. — 1965. — Т. 165, № 6.
5. Марчук Г.И. Методы расщепления. — М., 1988.
6. Ковеня В.М., Слюняев А.Ю. Схемы оптимального расщепления для решения уравнений Навье-Стокса // Проблемы и достижения прикладной математики и механики / РАН. Сиб. отд-ние.Ин-т теор. и прикл. механики. — 2010.
7. Воеводин А.Ф., Протопопова Т.В. Метод расчета вязких течений в замкнутых областях // Сибирский журнал индустриальной математики.
— 2001. — Т. 4, №1.
8. Воеводин А.Ф., Гончарова О.Н. Метод расщепления по физическим процессам для расчета задач конвекции // Вычислительные технологии.
— 2009. — Т. 14, №1.
9. Белолипецкий В.М., Костюк В.Ю., Шо-кин Ю.И. Математическое моделирование течений стратифицированной жидкости. — Новосибирск, 1991.
10. Воеводин А.Ф., Юшкова Т.В. Численный
метод решения начально-краевых задач для уравнений Навье-Стокса в замкнутых областях на основе метода расщепления // Сибирский журнал вычислительной математики. - 1999. - Т. 2, №4.
11. Мызникова Б.И., Тарунин Е.Л. Процессы установления стационарных конвективных течений в кубической полости при подогреве снизу // Нестационарные процессы в жидкостях и твердых телах. — Свердловск, 1983.
12. Бессонов О.А., Брайловская В.А., Никитин С.А., Полежаев В.И. Тест для численных решений трехмерной задачи о естественной конвекции в кубической полости // Математическое моделирование. — 1999. — Т. 11, №12.
13. Тарунин Е.Л. Вычислительный эксперимент в задачах свободной конвекции. — Иркутск, 1990.
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. Калиткин Н.Н. Численные методы. — М., 1978.
