CC BY
19
ЭЛЕКТРОМЕХАНИКА И ЭНЕРГЕТИКА
УДК 517.926.22
ОБ ОДНОМ РАЗНОСТНОМ МЕТОДЕ КОМПЬЮТЕРНОГО МОДЕЛИРОВАНИЯ СВОБОДНОЙ КОНВЕКЦИИ В ТЕПЛОЭНЕРГЕТИЧЕСКИХ УСТРОЙСТВАХ
© 2009 г. Н.С. Бузало, А.Н. Никифоров, О.К. Нилова, СА. Семенченко
Южно-Российский государственный South-Russian State
технический университет Technical University
(Новочеркасский политехнический институт) (Novocherkassk Polytechnic Institute)
Предложен новый численный алгоритм решения трехмерной нестационарной краевой задачи, описывающей динамику среды в условиях свободной конвекции в приближении Буссинеска, основанный на конечно-разностном методе маркеров и ячеек. Проведен вычислительный эксперимент для модельной задачи, выполнен анализ полученных результатов, определены приемлемые параметры разностной схемы, сделаны выводы о применимости метода в различных задачах.
Ключевые слова: динамика жидкости и газа, естественная конвекция, численные методы, конечно-разностные методы.
New numerical method of solution of three-dimensional nonstationary boundary value problem of fluid dynamics in free convection in Boissinesq approximation based on finite difference method of markers and cells is proposed. Computing experiment for model problem is taken and analysis of received results is performed. Acceptable parameters of difference scheme are identified. Conclusions on applicability of the method are made.
Keywords: fluid dynamics, free convection, finite difference methods, numeric methods.
Введение
Естественная конвекция в замкнутых и незамкнутых областях характерна для многих технических приложений. Например, свободноконвективные циркуляционные потоки в жидкостях оказывают существенное влияние на характер роста кристаллов. Значительный практический интерес представляет анализ процессов конвекции в полостях теплообменных устройств и топливно-энергетических установок. Также конвективные течения широко изучаются при анализе пожаров, при проектировании зданий, печей, систем аккумулирования и отвода энергии, а также некоторых других сооружений и промышленных устройств [1]. К сожалению, работ, в которых исследуются трехмерные задачи тепло- и массопереноса, причем как теоретических, так и экспериментальных, к настоящему времени опубликовано сравнительно немного, хотя уже разработаны некоторые общие методы построения численных алгоритмов решения трехмерных уравнений сохранения массы и количества движения [2]. Идея работы заключается в применении известного метода маркеров и ячеек [2] к уравнениям свободной конвекции в приближении Буссинеска.
Краевая задача
Пусть в трехмерном пространстве задана прямоугольная область О, заполненная жидкостью. Поверхность области состоит из нижнего основания (при у = 0), верхнего основания SH (при у = Н) и боковой поверхности S = SL и SR и SP и SZ , где
, SR, SP, SZ - левая, правая, передняя и задняя грани области соответственно (рис. 1).
Z
H
Q
X
Рис. 1. Схема расчетной области
Предполагается, что границы области теплоизолированы. На нижнем основании S0 имеется изотермический нагревательный элемент ST температуры T
1 hot ■
Плотность жидкости, за исключением плотности в членах архимедовых сил, которые вызывают свободную конвекцию, считаем постоянной, тх. используем приближение Буссинеска. Тогда для свободноконвек-тивных течений с малыми изменениями плотности уравнения сохранения массы, количества движения и
H
T
0
энергии в области О можно записать в следующем виде:
ди д¥ дW
dU dU-+-
dt dx
- +— +-= 0;
dx dy dz
d(UV) d(UW) dp
+ —-- + —-- = —— + vAU :
dy dz dx
(1)
dV d(UV) dV2 — + —-- +-
dt dx dy
d(vw) dp
+ —-- = —— + vAV ;
dz dy
dW d(UW) + d(VW) + dW2 dt dx dy dz
..-<£. + +VAW + ßg(T - Tq) dz
dT + d(UT) + d(VT) + d(WT) dt dx dy dz
PCp
-AT,
V = V0 , T = T0 при t = 0;
V = 0 на S u S0 u SH
dT
— = 0 на S u S0 u SH . dn
(6)
(7)
(8)
граничными условиями при сложной геометрии течения оперировать значительно легче.
Для построения разностной схемы выбрана полностью разнесенная трехмерная сетка, элементарная единица которой (ячейка) представлена на рис. 2.
(2) (3) / 'WK j, k-/
Vi, j+/ / / k Ui+ /, j, k
(4) hz (5) Ui-/, j, k / / У Pi, j, k Ti, j, k
/
Vi, j -'/2, k /
где и, V, W - компоненты вектора V скорости движения среды; р - коэффициент объемного расширения воздуха; g - ускорение свободного падения; р -аналог давления; р - плотность; Т - абсолютная температура; Т0 - среднее значение абсолютной температуры; V - коэффициент вязкости; X - коэффициент температуропроводности; Ср - теплоемкость воздуха
при постоянном давлении.
Кроме того, необходимы начальные условия для температуры и скорости и соответствующие граничные условия. Скорость на границах области обращается в ноль в силу выполнения условий непротекания и прилипания жидкости к стенкам. Градиент температуры в области отличен от нуля только в одном направлении, поэтому на границах выполняется условие адиабатичности, т.е. система не получает теплоты извне и не отдает ее через границы данной области. Граничных условий для давления на твердой поверхности задавать не надо. Таким образом, начальные и граничные условия для данной задачи следующие:
hx
Рис. 2. Расположение конечно-разностных переменных в ячейке
Различные зависимые переменные определяются в разных точках сетки: давление определяется в центре ячейки, а компоненты скорости - на границах.
В результате дискретизации уравнения (1), получаем выражение
Un
j,k
- Un
Vn+1 - Vn+V j+/2'k и-Уг,к
W.
n+1
•\k + V
- W
n+1
j k - К
h
= 0,
(9)
Порядок аппроксимации (9) равен О(Кх2, Ку2, Иг2).
Дискретизация уравнения (2) проводится с помощью конечно-разностных выражений, центрированных относительно точки сетки (г +1/2, ], k). Это позволяет представить др / дх в виде выражения (Рг+\,]Л - Рг,)/ К, которое в точке (г +1/2, ], k)
имеет второй порядок точности. Аналогично (3) и (4) аппроксимируются с помощью центральных разностей относительно точек (г, ] +1/2, k) и (г, ],k +1/2) соответственно, и др / ду представляется в виде (Рг,] +и - Рг,)/ КУ , а дР / дг - в виде (Рг,],k+1 - Рг,],k )/ ■
Использование разнесенной сетки дает возможность связать значения и, V, W и р в соседних точках. Это также позволяет избежать появления осцилляций
Разностная схема
Метод маркеров и ячеек (МИЯ) является одним из наиболее ранних и получивших широкое распространение методов решения уравнений сохранения массы (1) и количества движения (3), (4) [2]. В методе используется разнесенная сетка и на каждом шаге по времени решается уравнение Пуассона для определения давления. МИЯ применяет основные переменные, благодаря чему его возможно использовать при решении пространственных задач, так как в этом случае
X
h
у W., k
Y
jk
+
h
h
x
+
в решении, в частности для р, которые могут возникнуть, если центральные разности используются для аппроксимации всех производных на неразнесенной сетке. Осциллирующее решение появляется из-за двух несвязанных на различных точках сетки решений для давления, которое возможно, если центральные разности используются на неразнесенной сетке. Осцилляции, как правило, возрастают при увеличении числа Рейнольдса, поскольку диссипативные члены, посредством которых осуществляется связь значений компонент скорости в соседних точках, в этом случае малы
[3].
Таким образом, при дискретизации уравнения (2) в точке (/' + 1/2, ], к) используются следующие конечно-разностные выражения:
ди ^ Ui"+/2, j,k ~Ui+l/2, j,k dt ~ h
+ O(h);
dU2 ^ UL,j,k -Ujjkk u 2,.
dx h„
- + O(hxz):
d(UV )JUV h+i/u+i/ik -UV )i+Vi,j-_ yi,k +O(hi)
dy
К
d(UW) (UW)j+/i,j,k+/i <UW X+1/2, j,k-i/i „„i,
- » - +O(hz )
dz h.
dlU Ui+3/i,j,k-lUi+iji,j,k +Ui-Vi,j,k , 2
dx2
+O(h2):
d+Vl,j-1,k lUi+Vl,j,k +Ui+1/2,j+1,k +0(h2) .
dy2
h,
d 2u Ui+1/2, j,k -1 2Ui+1/2, j,k +Ui+1/2, j,k+1 dz2 ~ Г2 +O(hz ) ;
где
F'+1/2, j,k-U,'+1/2, j,k + ht [VISX FUX FUY FUZ ]; _ 1 2
FUX (Ui+1/2, j,k + 2Ui+1/2, j,kUi+3/2, j,k +
+U,
i+3/2, j,k 4Ulj,k )
V
(
FUY —
(
U 1 +U 1
i+—, j,k i +—, j+1,k V 2 2 У
Л
V 1 +V 1
i+1, j—,k i, j +—,k V 2 2 У
Л
4hy
U 1 + U 1 i+—, j-1,k i+—, j,k V 2 2 У
Л
V 1 + V 1
i+1, j—,k i, j—,k V 2 2 У .
4h„
U , +U 1
FTI7 —
i+—, j,k i+—, j,k+1 V 2 2
W 1+W 1
i, j,k+— i+1, j ,k+— V 2 2 У
f
\
U 1 + U 1
i+—, j,k-1 i +—, j,k V 2 2 У
4hz
Г
Л
W 1 + W 1
i+1, j,k— i, j,k— V 2 2 У .
4hz
U,+w,, -2U + U ,+3/2, j ,k /+1/2, j ,k i-1/2, j ,k
VISX =v
Ui+1/2, j-1,k 2Ui-+1/2, j,k + Ui-+1/2,j+1,k
U,+1/l,, 1-2U + U i+,1/2, j,k-1 i+1/2, j,k i+1/2, j ,k+1
Аналогично в дискретном виде представляются уравнения (3) и (4):
т^И+1 И+1 И+1 \ . /1 1 л
Ч] +1/2,к = Ч,] +1/2,к ]+1,к " Р 1,к) ' (11)
Зр (Рг+1,],к Ри],к I 2
тт *—^-—+^).
ох hx
В приведенных выражениях присутствуют не определенные на рис. 2 члены типа Ui+1 ]к . Их аппроксимация осуществляется следующим образом:
U
i+1, j,k
- ^^Ш, jkk + Ui +3/2, jkk ).
'[(V
, j+1/2,k + V', j+1/2,k ) /2
Решение уравнения (2) осуществляется по следующей явной схеме:
Г гИ+1 П+1 И+1 \ /1П\
иг+1/2,],к = Л+1/2,],к 7 уРг+1,],к " Р],к ^ , (10)
W
n+1
= H"
ht (pÜM+1- Pj ) , (12)
где
G
i, j+1/ 2,k -V,", , 1/. +ht [VISY -FVY -FVZ-FVX ];
l, J+ 7r.,k
Hi, j,k+1/2 - Wi,j,k+1/2 + +ht[VISZ-FWX-FWY-FWZ-_2_(Ti,j,k +1+ Ti,j,k-2T0 )
Аналогично, например, (ЦУ),+у2 ]+у2 к аппроксимируется выражением
(ЦУ)г +12,]+1/2,к =[(Цг+12,],к + Цг+1/2,] +1,к ) /2
^^УХ , FWX ; ^УГ , ; FУz , FWZ ; УШ , У^SZ опреде-
ляются аналогично FUX , FU7 , FUZ , У}^.
Разностное представление уравнения сохранения энергии (5) может быть записано в следующей форме:
rpn + 1 _гр n .1
Ti, jM~Li, j,k+ht
X-
TSX +TSY +TSZ
PCp
TUX TVY TWZ
(13)
где
+
2
h
X
h
2
h
X
h
z
tux —
Ui+1/2, j,k (Ti+1,j,k + Ti, j,k ) ~2h„
TVY —
TWZ -
Ui-1/2, j,k {Ti-1,j,k + Ti, j,k ) ; " 2hx '
Vi, j+1/2,k (Ti, j+1,k + Ti, j,k ) lh„
V', j-1/2,k (Ti, j-1,k + Ti, j,k ) ; " 2hy '
Wi, j,k+1/2 {^i, j,k+1 + Ti, j,k ) 2h
Wi, j,k-1/2 (Ti, j,k-1 + Ti, j,k ) ^
TSX -
TSY —
TSZ —
2hz
Ti-1,j,k - 2T,j,k + Ti+1, j,k
К2
Ti, j-1,k - 2T,j,k + Ti, j+1, k
hy 2
Ti, j,k-1 - 2T,j,k + Ti, j,k+1
h
Pi, j,k-1 - 2Pi, j,k + Pi, j,k+1
(14)
T7n —T7n Pi+1/2,j,k Fi-12,j,k Gi,j+1/2,k Gi,j-1/2,k +
ттП ТТП
Hi, j,k +1/2 - Hi, j,k-1/2
2 1
0,25(| U| + \V\ + \W\) ht— < 1 и ht v / Ar2 < 0.25, (15)
где Ax = max |hx, hy, hz j.
Для решения уравнения (14) необходимо поставить граничные условия для давления p (условия Дирихле) или для производных от давления p (условия Неймана) на всех границах. В последнем случае необходимо выполнение глобального граничного условия, т.е.
ш
Q
222 О p О p О p
-7Г + —7Г + —Т-
Ок2 Oy
Oz2
Op
dxdydz — ff — ds ,
dndn
где интеграл по дО - это интеграл по поверхности расчетной области. Левая часть этого уравнения в дискретном виде вычисляется через правую часть (14).
Для замкнутой области глобальное граничное условие соответствует дискретному условию
Q
i, j,k
— 0.
(16)
В выражения (10) - (12) давление р входит неявно; однако р"+1 определяется до решения (10) - (13)
следующим образом. Обратимся к уравнению неразрывности (1), записанному в разностном виде (9). Подстановка уравнений (10) - (12) позволяет представить (9) в виде разностного уравнения Пуассона для давления, т.е.
Рг-1],к - 2Рг,],к + Рг+1,],к
К2 +
+ Рг,]-1,к - 2Рг,},к + Рг,]+1,к +
Для наложения необходимых граничных условий жидкость считается окруженной одним слоем фиктивных ячеек, в которых соответствующие переменные заданы. Виды использованных здесь граничных условий проиллюстрированы на примере левой границы.
1. Ячейка с твердой стенкой и условием прилипания (условие (7)):
U,
— 0,
V1,j,k — V2, j,k
Whhk —-Wj.
К
Уравнение (14) решается на каждом шаге по времени. После того как р"+1 получено, подстановка этого значения в выражения (10) - (12) позволяет определить , }2,к, W}k+У2.
Поскольку выражения (10) - (13) являются явными формулами, имеется ограничение на максимальный шаг по времени, связанное с устойчивостью решения [2]:
2. Ячейка с адиабатической стенкой (условие (8)):
Т1,],к = Т2,],к .
При решении уравнения Пуассона для давления (14) требуются его значения за пределами области расчета. При записи (14) относительно узла (2, 1, 1) требуются значения р2 01 и V2 -1/21. Значение р2 01 получается из уравнения (3), записанного на стенке.
Численный алгоритм
Решение задачи сводится к следующим этапам.
1. Ввод исходных данных:
- размеры области, нагревательного элемента;
- шаги разбиения;
- плотность жидкости (р), коэффициент объемного расширения жидкости (Р), коэффициент вязкости (у), коэффициент теплопроводности жидкости (X), удельная теплоемкость жидкости (Ср);
- значения функций компонент скорости (и,УЖ), давления (р) и температуры (Т) в начальный момент времени из (6).
Для вычисления температуры, давления и скорости на новом временном шаге по значениям на предыдущем временном шаге осуществляется следующий цикл.
2
h
y
2
h
z
h
h
x
+
2. Вычисление температуры Т"+1 по формуле (13).
3. Решение уравнение Пуассона для давления (14). Для этого используем, например, блочно-
итерационный метод решения эллиптических разностных уравнений. Одна итерация включает в себя последовательное решение трех уравнений, получаемых из (14):
(рг-1,],к - 2рг,],к + рг+1,],к ) + 2Рг,],к + 2Рг,],к _
ттп _ТТП s-~<n
Fi+1/2, j,k~ri-1/2, j,k + Gi,j+1/2,k~Gi,j-l/2,k +
Hn - Hn
ni, j,k+1/2 ni, j,k-1/2
pi, j-1,k + pi, j+1,k pi, j,k-1 + pi, j,k+1 K2 h}
(pi,3-1,k - 2pi, j,k + pi, j+1,k ) | 2Pi, jk + 2Рг, j,k
K2 K2
n n n n
Pi+1/2, j,k Fi-1/2, j,k Gi, j+1/2,k Gi, j-1/2,k +
Hn - Hn
i, j,k+1/2 i, j,k-1/2
pi-1, j,k + pi+1,j,k pi, j,k-1 + pi, j,k+1
K
K
(pi, j,k-1 - 2pi, j,k + pi, j,k+1 ) i 2Pi, j,k + 2Pi, j,k
К2 K2
n n n n
pi+1/2,j,k Fi-1/2,j,k Gi,j+1/2,k~Gi,j-1/2,k +
Hn - Hn
ni, j,k+1/2 ni, j,k-1/2
pi, j-1,k + pi, j+1,k pi, j,k-1 + pi, j,k+1
К
К
Здесь значения давления, входящие в левые части уравнений, неизвестны, а значения, входящие в правые части уравнений, считаем известными и берем их с предыдущего подшага. Записывая эти уравнения для каждой точки области, получаем три системы линей-н^1х алгебраических уравнений с трехдиагональными матрицами. Для их решения применяем метод прогонки. Выполняем описанные итерации до достижения приемлемой невязки уравнения (14).
При использовании граничных условий Неймана
др
— _ 0 матрица системы (14) является вырожденной.
дп
Для того чтобы избавиться от вырожденности матрицы, в нескольких точках заменяем граничные условия Неймана на граничные условия Дирихле р _ 0 . Количество точек замены зависит от метода решения системы линейных алгебраических уравнений (14). При рассмотренном блочно-итерационном методе необходимо задать условия Дирихле в точках, лежащих на тех границах расчетной области, где х = 0, у = 0, г = 0.
4. Проверяем глобальное условие для давления (16). Его выполнение характеризует зону устойчивости построенной разностной схемы.
5. Полученные значения температуры и давления подставляем в выражения (10) - (12) для получения значений компонент скорости иЩ^}- к , V"+l1/2 к,
Wn+1 УЧ,],к+1/2 .
6. Выполняем проверку необходимого условия устойчивости решения (15).
7. Используем вычисленные поля в качестве исходных для вычислительного цикла на следующем временном шаге.
Вычислительный эксперимент и анализ результатов
Проведен расчет полей температуры, давления и скорости движения среды в ограниченной области при свободной конвекции для следующей модельной задачи. Область расчета представляет собой прямоугольную замкнутую поверхность размером 1 х 1 х 0,2 м. Изотермический нагревательный элемент в форме квадрата расположен в центре нижней стенки полости. Площадь квадрата равна 0,04 м2. Стенки полости (включая оставшуюся часть нижней) поддерживаются при постоянной температуре Т0 _ 0. Температура нагревательного элемента Т _ 20 °С.
Для данной области оптимальными являются шаги по пространству порядка 0,01 м и шаг по времени К _ 10-5 с. Уменьшение как пространственных, так и временных шагов ведет к увеличению времени работы программы до 10 раз. Увеличение шагов - к невыполнению условий устойчивости решения. Модельное время расчета составило 600 с. За этот период времени установившийся режим не наступает. Для достижения установившегося режима необходимо нецелесообразное количество времени.
На рис. 3. приведено распределение поля температуры и (и, IV) составляющей вектора скорости движения среды в плоскости у _ 0,5 м (] _ 50), г _ 600 с.
Тестовые расчеты показали, что построенная разностная схема пригодна для решения многих задач о свободноконвективных течениях жидкости в замкнутых областях. Однако в алгоритме имеются жесткие ограничения на размерность шагов по времени и пространству.
h
h
h
X
z
1
h
x
h
z
2
h
h
h
X
h
z
2
2
2
h
z
1
h
h
h
X
X
h
z
2
2
Рис. 3. Распределение поля температуры и векторного поля скорости в плоскости у = 0,5
По всей видимости, это вызвано тем, что, во-первых, используемая разностная схема является явной, во-вторых, отсутствует этап предиктор - корректор. Поэтому для расчетов характеристик течения в отдаленные моменты времени необходимо выполнить большое число итераций, что приводит к накоплению вычислительной ошибки и требует большого времени счета. Таким образом, разработанный алгоритм пригоден для выполнения только локальных расчетов движения среды и только на малых интервалах времени. К сожалению, по этой причине алгоритм не позволяет определять характеристики течений, близких к установившимся.
Однако для локальных по времени расчетов то, что используемая разностная схема является явной и однослойной, дает ей преимущество во времени расчетов по сравнению с известными алгоритмами типа SIMPLE («предиктор - корректор») и неявными схема-
Поступила в редакцию
ми. Благодаря быстродействию составленный численный алгоритм можно использовать как составляющую часть других алгоритмов, например оптимизационных или мультипликационной визуализации данных.
Литература
1. Российская наука на заре нового века: сб. науч.-популярных ст. / под ред. В.П. Скулачева. М., 2001. 496 с.
2. Флетчер К. Вычислительные методы в динамике жидкости: в 2 т. Т. 1 : пер. с англ. М., 1991. 504 с.
3. Вабищевич П.Н., Павлов А.Н., Чурбанов А.Г. Численные методы решения нестационарных уравнений Навье -Стокса в естественных переменных на частично разнесенных сетках // Мат. моделирование. 1997. Т. 9 № 4. С. 85-114.
19 марта 2009 г.
Бузало Наталья Сергеевна - канд. техн. наук, доцент, кафедра «Прикладная математика», Южно-Российский государственный технический университет (Новочеркасский политехнический институт). Тел. (8635) 25-53-09. E-mail: buzalo.n.s@mail.ru
Никифоров Александр Николаевич - канд. техн. наук, профессор кафедра «Прикладная математика», ЮжноРоссийский государственный технический университет (Новочеркасский политехнический институт). Нилова Ольга Константиновна - студентка, кафедра «Прикладная математика», Южно-Российский государственный технический университет (Новочеркасский политехнический институт).
Cеменченко Степан Анатольевич - канд. техн. наук, доцент, кафедра «Общая и прикладная физика», ЮжноРоссийский государственный технический университет (Новочеркасский политехнический институт). Buzalo Natalya Sergeevna - Candidate of Technical Sciences, assistant professor, department «Applied mathematics», South-Russian State Technical University (Novocherkassk Polytechnic Institute). Ph.(8635) 25-53-09. E-mail: buzalo.n. s. @mail.ru
Nikiforov Alexandr Nikolaevich - Candidate of Technical Sciences, professor, department «Applied mathematics», South-Russian State Technical University (Novocherkassk Polytechnic Institute).
Nilova Olga Konstantinovna - student, department «Applied mathematics», South-Russian State Technical University (Novocherkassk Polytechnic Institute).
Semenchenko Stepan Antonovich - Candidate of Technical Sciences, assistant professor, department «General and applied physics», South-Russian State Technical University (Novocherkassk Polytechnic Institute).
