close

Вход

Забыли?

вход по аккаунту

?

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

код для вставкиСкачать
МАТЕМАТИКА И МЕХАНИКА
УДК 532.516
О.Н. Гончарова
Расщепление по физическим процессам для
расчета задач конвекции в двумерных
областях с криволинейной границей
Для исследования конвективных движений
жидкости в двумерных областях с твердой непроницаемой границей, часть которой является
криволинейной вставкой, предлагается численный метод, основанный на идее расщепления по
физическим процессам. Расщепление по физическим процессам в задачах гидродинамики
базируется на методе слабой аппроксимации и
обосновании аддитивности этих процессов при
достаточно малых шагах по времени [1-3].
Расщепление на два этапа (конвективный и
диффузионный переносы) проводится для уравнений конвекции, записанных в физических переменных. Этап конвекции реализуется на смещенных сетках для компонент скорости и состоит в вычислении вспомогательной функции. На
этапе диффузии осуществляется переход к новым искомым функциям: вихрю и функции тока. Такой подход позволяет исключить расчет
градиента давления и обеспечить автоматическое выполнение условия соленоидальности вектора скорости.
Тестирование метода проводится на известных задачах о свободной конвекции в полости
при подогреве сбоку.
Введение. Для моделирования тепловой
гравитационной конвекции в замкнутых областях с твердой непроницаемой границей применяются уравнения Обербека-Буссинеска,
которые в безразмерной форме имеют следующий вид:
1
~vt + ~v ╖ ?~v = ??p0 +
~v + F~ ,
(1)
Re
div ~v = 0,
1
T.
Tt + ~v ╖ ?T =
ReP r
(2)
(3)
Задача состоит в определении скорости жидкости ~v, давления p и температуры T , удовлетворяющим при t = 0 начальным условиям
= 0, T = T0 (x), x ? ,
на границе области | условиям прилипания
для скорости
~v
~v
= 0,
x ? , t ? [0, t? ]
34
и условиям первого, второго или третьего рода
для температуры:
T
?T
?n
= TB (x, t),
?T
?n
+ ?(T ? TB ) = 0;
= TB (x, t),
x ? , t ? [0, t? ].
Здесь p0 | модифицированное давление (отклонение от гидростатического давления); F~ =
?(Gr/Re2 )g~0 T { сила плавучести; g~0 { единичный вектор ускорения силы тяжести; ? { безразмерный коэффициент, характеризующий теплообмен с внешней средой.
Конвективный и диффузионный переносы
выделяются в уравнениях движения (1) и в
уравнении переноса тепла (3). Идея расщепления по физическим процессам, как основа для
организации расчета свободной конвекции, реализуется на каждом временном слое. Метод расщепления по физическим процессам базируется на методе слабой аппроксимации [1], обосновании аддитивности процессов для достаточно
малых шагов по времени [2] и доказательстве
суммарной аппроксимации исходного уравнения
вследствие расщепления. Общая теория расщепления наиболее полно изложена в [3], а особенности расщепления для задач конвекции в прямоугольных областях и параллелепипедах { в [4{
6].
В данной работе численный метод на основе
расщепления по физическим процессам разрабатывается для исследования конвективных движений жидкости в замкнутой двумерной области с криволинейной границей. Метод характеризуется свойством энергетической нейтральности
поля скоростей [4]. Расщепление на два этапа
(конвективный и диффузионный переносы) проводится для уравнений конвекции, записанных
в физических переменных. Этап конвекции реализуется для вектора скорости, условно названной конвективной скоростью. Выделение этапа конвекции в исходных уравнениях позволяет обеспечить корректность расщепления с точки зрения удовлетворения граничным условиям.
Из условий прилипания и гиперболичности системы уравнений следует, что граничные условия на этапе конвекции являются следствием
Метод расщепления для двумерных задач конвекции . . .
уравнений. На этапе диффузии осуществляется переход к новым искомым функциям: вихрю и функции тока. Данный переход позволяет исключить расчет градиента давления, обеспечить автоматическое выполнение условия соленоидальности вектора скорости, так как диффузионная скорость принимается в итоге в качестве искомой. Для функции тока считаются выполненными следствия условий прилипания для физической скорости. В дальнейшем
будем опускать кавычки в терминах конвективная или диффузионная скорость. Отметим, что
принципиальным моментом в организации расчета является введение смещенных сеток [7].
Тестирование метода проводится на известной задаче о свободной конвекции в квадратной
полости при подогреве сбоку [8].
Постановка разностной задачи. Пусть
искомые функции ~v, p0 , T удовлетворяют в области = {(x, y) : 0 < x < 1, 0 < y < f (x)}
системе уравнений (1){(3). Перепишем данные
уравнения в виде
?~v
?t
+ K~v = ??p0 + D~v + F~ ,
?T
?t
(4)
div ~v = 0,
(5)
+ KT = DT.
(6)
Здесь K { оператор конвективного переноса K =
~v ╖ ?; D { оператор диффузионного переноса
D = ?, ? = 1/Re или ? = 1/(ReP r).
В уравнении (4) выделим два различных процесса переноса количества движения: конвективный и диффузионный. Считая эти процессы
аддитивными на малых временных интервалах
tk ? t ? tk+1 (? = tk+1 ? tk ), сформулируем
задачу конвективного движения жидкости, исходя из нижеследующих двух этапов (опустим
стрелочку над векторами). Назовем этапом конвекции движение, описываемое уравнениями
v~ ? v
?
= ?K v~,
(7)
где v~ { конвективная скорость. При построении
разностной аппроксимации оператор конвективного переноса будем рассматривать в виде суммы двух одномерных кососимметричных операторов: K v~ = K1 v~ + K2 v~ [3].
Этап диффузии определяется уравнением
v^ ? v~
?
= ??p0 + Dv^ + F,
35
которое в результате применения операции rotor
преобразуется к виду:
?
^ ? ?~
= D?^ + (rot F, k).
(8)
?
Здесь под v^ понимается диффузионная скорость процесса и используются обозначения ?^ =
^ ,
(rot v^, k), ?~ = (rot v~, k), v^ = (^u1 , u^2 ), u^1 = ? ?/?y
^
^
u
^2 = ?? ?/?x. При этом функция тока ? удовлетворяет в области уравнению Пуассона:
?^ = ??^ ,
(9)
которое решается в дальнейшем методом установления [9]. На границе области имеют место
условия
? ?^
?^ = 0,
= 0,
(10)
?n
представляющие собой следствия условий прилипания для физической скорости.
Расщепление на конвективный и диффузионный этапы в уравнении (3) приводит к последовательному решению следующих разностных
уравнений:
T~ ? T
= ?K T~,
(11)
?
T^ ? T~
?
= DT^.
(12)
Для решения задачи в области с криволинейной границей осуществим, как и в [10], отображение области течения на квадрат = {0 ?
? ? 1, 0 ? ? ? 1}, применив преобразование
? = x, ? = y/f (x). Перепишем уравнения (8),
(9), (12) в переменных ?, ? в следующем общем
виде:
i
Bh
t = U? + V? + R,
(13)
где = ?, =
= T . При этом
U
f
?
(крышка опускается) или
= B11 ? + B12 ? , V = B12 ? + B22 ? ,
а коэффициенты вычисляются следующим образом:
B11 = f, B12 = ??f 0 ,
2
B22 = (1 + ? 2 f 0 )/f, B = 1/Re,
б
в
R = (Gr/Re2 ) ?T /?? ? ? (f 0 /f )?T /?? ,
если = ?. Если же = ?, то B = ? и R = ??
(? | итерационный параметр). Если = T , то
B = 1/(ReP r), а R = 0.
Конвективные операторы K1 , K2 в переменных ? , ? имеют вид:
1 h ?v ? (v1 v) i
+
,
K1 v = v 1
2 ??
??
МАТЕМАТИКА И МЕХАНИКА
K2 v
=
= K 2 v +
f0
v v
2f 1
+
=
1h
?v
(v ? ?f 0 v1 ) +
2f 2
??
В области вводится основная разностная сетка (?n , ?m ) для расчета температуры и реализации этапа диффузии (9), (12): ?n = (n ? 1)h? ,
n = 1, 2, ..., N ; ?m = (m ? 1)h? , m = 1, 2, ..., M ;
+ 1; h? = 1 , M = M + 1. В
h? = N1 , N = N
M
расчетной области вводятся также две вспомогательных разностных сетки: (?n , ?m ), (m0 =
) { для расm + 1/2, n = 1, ..., N, m0 = 1, ..., M
чета первой компоненты конвективной скорости
, m =
v~1 , и (?n , ?m ), n0 = n + 1/2, n0 = 1, ..., N
1, ..., M ) { для расчета второй компоненты конвективной скорости v~2 . Здесь v~1 , v~2 { компоненты конвективной скорости v~.
~ij
K1h w
2h w~ij
K
v2 n+1/2,m
??m
?n,m+1 ? ?nm
,
fn
h?
│?
=?
n+1,m
? ?nm
h?
f 0 n ?n,m+1 ? ?nm ┤
.
fn
h?
Для решения уравнений (7), (11) используется
двуцикличеcкий способ расщепления по направлениям на основе элементарных схем Кранка{
Николсона [3]. Представим данную схему для
искомой функции w~ на произвольной сетке
(i, j ), которая может быть основной (n, m) для
температуры либо одной из вспомогательных
(n0 , m), (n, m0 ) для компонент конвективной
скорости:
│w
~ijk+1/4 ? wijk ┤ h │w~ijk+1/4 + wijk ┤
+ K1
=0,
? /2
2
│w
~ijk+1/2 ? w~ijk+1/4┤ h │w~ijk+1/2 + w~ijk+1/4┤
+ K2
+
? /2
2
f0
v w
~ k+1/4 = 0,
2f 1 ij
│w
~ijk+3/4 ? w~ijk+1/2┤ h │w~ijk+3/4 + w~ijk+1/2┤
+ K2
+
? /2
2
+
36
2h ?
│
v1 ki+1/2,j w
~i+1,j ?
=
1 │ k
v2 i,j +1/2 w
~i,j +1 ?
fi 2hy
fi+1 ? fi?1 ?j v1 kij
?
2h ?
+ ?j +1 v1 ki,j +1
w
~i,j +1+
2
fi+1 ? fi?1 ?j v1 kij
2h ?
┤
+ ?j?1 v1 ki,j?1
w
~i,j?1 .
2
+
?
1
?v2 ki,j?1/2 w
~i,j?1 ?
Реализация этапа конвекции. Началь-
1
=
┤
?v1 ki?1/2,j w
~i?1,j ,
0
ные условия для конвективной скорости v~ задаются с предыдущего временного слоя и представляют собой диффузионную скорость, рассчитанную на вспомогательных сетках. При этом
(14)
Здесь K1h , K 2h представляют собой конечноразностную аппроксимацию конвективных операторов. Во внутренних точках сеток (n, m),
(n0 , m) или (n, m0 ) рассматривается аппроксимация, использующая центральные разности вида
0
=
= 0,
│w
~ijk+1 ? w~ijk+3/4┤ h │w~ijk+1 + w~ijk+3/4┤
+ K1
=0.
? /2
2
i f0
?
+ (v2 ? ?f 0 v1 )v + v1 v.
??
2f
v1 n,m+1/2
f0
v w
~ k+3/4
2f 1 ij
Значения компонент скорости с предыдущего
k -го слоя рассчитываются по-разному для каждой расчетной сетки в зависимости от искомой функции. Заметим, что диффузионная скорость (u1 , u2 ) уже заранее рассчитана на основной сетке (n, m). В граничных точках используется аппроксимация направленными разностями
аналогично тому, как это было сделано в [5]. Заметим также, что при расчете компонент конвективной скорости на смещенных сетках учитываются в коэффициентах условия прилипания
при выходе за границу вспомогательной сетки и
при попадании на границу основной сетки. Используемая аппроксимация конвективных операторов обеспечивает кососимметричность разностных операторов конвекции. Полученные системы разностных уравнений в (14) решаются
с помощью метода прогонки. Хотя для матрицы системы и не выполнены условия диагонального преобладания, метод прогонки устойчив [11]. Вышеизложенные аппроксимации операторов K1 , K2 имеют второй порядок [3]. Используемая схема обладает свойством энергетической нейтральности поля скоростей, т.е. сохранением нормы скорости в L2 (для доказательства см.: [4]).
Метод расщепления для двумерных задач конвекции . . .
Реализация этапа диффузии. При переходе на этап диффузии насчитывается вихрь
конвективной скорости ?~ на основной сетке
w
~n,m
v2 n+1/2,m ? v2 n?1/2,m
?
h?
=
??m
╖
?
fn+1 ? fn?1 1
╖
f h?
fn
v2 n,m+1/2 ?v2 n,m?1/2
?
h?
?n,1
1
v1 n,m+1/2 ? v1 n,m?1/2
,
fn
h?
?
i
Bh k
U ? + V k+1/2 ? +
f
=
+Rk+1/2 ,
k+1 ? k+1/2
?
=
i
B h k+1
k
U
? ?U ? .
f
(15)
Здесь используются представления
U k+1
= B11 k? +1 + B12 k?+1/2 ,
Uk
= B11 k? + B12 k? ,
V k+1/2
= B12 k? + B22 k?+1/2
и аппроксимация входящих в (15) производных
следующего вида:
б ?U k+1 в
n,m
??
=
hб
1
h?
(B )11
в
n+1,m
=
k+1 n+1,m ? k+1 n,m
k+1
б
)11 в nm
? (B
n,m
h?
k+1
?
n?1,m
h?
h
б
в
1
+
(B )12 n+1,m ╖
2h?
╖
╖
= 0,
Vn,2
и именно он берется в качестве начального условия для этапа диффузии ?^ k = ?~ .
Этап диффузии реализуется с использованием метода стабилизирующей поправки. Представим разностную схему второго порядка аппроксимации для решения уравнения (13) в следующем виде [1, 10]:
k+1/2 ? k
Для нахождения ? используется вышеприведенная схема для организации итерационного
процесса.
Воспользовавшись условиями (10) и уравнением (9), получим аналог условий Тома для вихря (см. также: [10]). В частности, на границе
? = 0 (при m = 1) эти условия имеют следующий вид:
?
i
+
k+1/2 n+1,m+1 ? k+1/2 n+1,m?1
?
2h?
б
)12 в
? (B
╖
n?1,m
k+1/2 n?1,m+1 ? k+1/2 n?1,m?1 i
.
2h ?
37
+
?n,1
=?
1 1
h? fn
Vn,2 ,
= ??2 (f 0 )n (?? )n,2 +
1 + ?2 2 (f 0 )2 n
fn
(?? )n,2 ,
где для первых производных функции тока
используются аппроксимации второго порядка.
Заметим, что условия превращаются в классические условия Тома в случае прямолинейной границы f (x) = 1.
Результаты численного анализа. Общая
схема решения задачи состоит в осуществлении
следующих этапов.
1. Переход на новый временной слой tk+1
начинается с расчета температуры T k+1
на основной сетке по схеме (14) для этапа
конвекции и по схеме (15) с соответствующими температурными граничными условиями на этапе диффузии.
2. При найденном T k+1 осуществляется
расчет конвективной скорости v~k+1 на
вспомогательных сетках и вихря скорости
?
~ k+1 | на основной.
3. На каждом k-м временном слое расчета ?^
по схеме (15) вводится внутренний итерационный процесс расчета ?^ . После окончания итераций при s = S считается, что
с заданной точностью ?? определены значения ?^ на (k + 1)-м временном слое:
?^k+1 = ?^S .
4. При переходе на следующий временной
слой (к этапу 1) искомая скорость полагается равной v^ (она известна на основной
сетке и может быть насчитана на вспомогательных).
Итерационный процесс для ?^ считается сошедшимся, если выполнен критерий сходимости [5,
8]
s+1
s
s+1
max
|?n,m
? ?n,m
| < ?? ╖ max |?n,m
|.
n,m
n,m
МАТЕМАТИКА И МЕХАНИКА
Рис. 1: Линии тока. Изотермы (справа)
Рис. 2: Линии тока. Изотермы (справа)
В качестве основного теста рассматривается
задача о свободной конвекции в полости 0 ?
x ? 1 ; 0 ? y ? f (x) при подогреве сбоку. При
этом Re = P r = 1, Gr = 104 , а вектор силы тяжести направлен против оси Oy. В начальный
момент времени жидкость покоится и нагрета по
закону T (x, y, 0) = x. Граничные условия соответствуют твердым непроницаемым границам с
заданной температурой T (x, y, t) = x . Тест заключается в получении стационарного решения
до удовлетворения неравенства [4, 5]
|N uk+1 ? N uk | < ?,
где число Нуссельта N u определяет интегральный поток тепла на правой (горячей) границе
x = 1 (см. [4, 5]). При тестировании полагается
? = 10?8 , ??
?
= ?0 = h2 /4,
= 10?5 ,
h = h?
Результаты тестирования алгоритма, а именно:
характеристика интенсивности течения ?max и
интегральный тепловой поток на правой стенке
x = 1 в случае прямолинейной границы f (x) = 1
приводятся в таблице, в последней строке которой представлены предельные значения характеристик из [8].
Сетка
9╫9
11 ╫ 11
17 ╫ 17
21 ╫ 21
33 ╫ 33
41 ╫ 41
[8]
?max
7.3558
7.0415
6.6450
6.5460
6.4360
6.4115
6.37
Nu
1.9405
1.8886
1.8057
1.7831
1.7607
1.7539
1.752
Устойчивость алгоритма и экспериментальный порядок сходимости проверяются путем вы-
= h? .
38
Метод расщепления для двумерных задач конвекции . . .
числительного эксперимента на последовательностях сеток 11 ╫ 11, 21 ╫ 21, 41 ╫ 41, при этом ведутся наблюдения за величиной ri , характеризующей для каждой сетки интенсивность движения max
|?n,m |. Экспериментальный порядок
n,m
сходимости r и приближенно определенная относительная погрешность вычисляются, согласно правилу Рунге (см., например: [12]). Таким
способом определяется r ? 1.9 и относительная
погрешность ▓ ? 2.8%.
Для криволинейных границ, задаваемых как
f (x) = 1 + 0.1 sin ?x и f (x) = ?0.1x + 1.1, результаты расчетов приведены на рисунках 1 и 2,
соответственно. На рисунках представлены семейства линий тока и изотерм.
Таким образом, проведено построение и тестирование численного алгоритма для исследования задач конвекции в двумерных областях с
криволинейной границей. Отметим также, что
реализация этапа диффузии может быть осуществлена на основе метода прогонки с параметрами [5].
Библиографический список
1. Яненко, Н.Н. Метод дробных шагов решения многомерных задач математической физики
/ Н.Н. Яненко. { Новосибирск, 1967.
2. Самарский, А.А. О принципе аддитивности для построения экономичных разностных
схем / А.А. Самарский // ДАН СССР. 1965. {
Т. 165, №6.
3. Марчук, Г.И. Методы расщепления / Г.И.
Марчук. { М., 1988.
4. Воеводин, А.Ф. Метод расчета вязких течений в замкнутых областях / А.Ф. Воеводин,
Т.В. Протопопова // Сиб. журнал индустриальной математики. { 2001. { Т. 4, №1.
5. Воеводин, А.Ф. Метод расчета двумерных
задач конвекции на основе расщепления по физическим процессам / А.Ф. Воеводин, О.Н. Гончарова // Вычислительные технологии. { 2002.
{ Т. 7, №1.
6. Гончарова, О.Н. Метод расщепления по
физическим процессам для расчета трехмерных
задач конвекции / О.Н. Гончарова // Известия
АлтГУ. 2007. { Вып. 1 (53).
7. Белолипецкий, В.М. Математическое мо-
39
делирование течений стратифицированной жидкости / В.М. Белолипецкий, В.Ю. Костюк,
Ю.И. Шокин. { Новосибирск, 1991. 176 с.
8. Тарунин, Е.Л. Вычислительный эксперимент в задачах свободной конвекции / Е.Л. Тарунин. { Иркутск, 1990.
9. Самарский, А.А. Теория разностных схем
/ А.А. Самарский. { М., 1977.
10. Овчарова, А.С. Метод решения термоконвективной задачи в многослойной среде с
криволинейными границами раздела / А.С. Овчарова // Динамика сплошной среды : cб. науч.
тр. АН СССР / Сиб. отд-ние. Ин-т гидродинамики. { 1994. { Вып. 106.
11. Воеводин, А.Ф. О корректности метода
прогонки для разностных уравнений / А.Ф. Воеводин // Численные методы механики сплошных сред: Сб. науч. тр. АН СССР / Сиб. отдние. Ин-т теор. и прикл. механики. { 1972. {
Т. 3, №5.
12. Калиткин, Н.Н. Численные методы /
Н.Н. Калиткин. { М., 1978.
Документ
Категория
Без категории
Просмотров
25
Размер файла
139 Кб
Теги
физическая, областям, конвекция, расщепление, двумерные, процесса, границе, расчет, задачи, криволинейных
1/--страниц
Пожаловаться на содержимое документа