УДК 519.6: 536.24
Метод численного моделирования трехмерного тепло- и массообмена при различных режимах течения.

Для численного решения задач естественной конвекции в двухмерной постановке при значительных числах Релея в работе [2] предлагается новый подход, который базируется непосредственно на уравнениях Навье-Стокса и методе торможения скоростей изменения искомых функций в некоторых узловых точках области, в которых эти скорости превышают допустимые с физической точки зрения значения. Идея метода торможения успешно использовалась при численном моделировании некоторых двумерных задач теплообмена при естественной [2] и вынужденной [3] конвекции несжимаемой и сжимаемой жидкости [4]. При моделировании течений на основе использования метода искусственной сжимаемости точность расчетов существенно снижается за счет того, что принимается постоянной во всей области течения характерная скорость, которая в действительности может существенно меняться в расчетной области.
В настоящей работе излагается метод численного моделирования динамики трехмерного течения и тепломассообмена вязкой жидкости при ламинарном, переходном и турбулентном режимах движения. Для расчета переходного и турбулентного режимов движения жидкости привлекается метод торможения скоростей изменения искомых функций.
Математическая модель и метод расчета. Система уравнений, описывающая трехмерное течение неоднородной жидкости и тепломассоперенос в приближении Буссинеска в переменных - вектор вихря го, вектор потенциала у, температура Т, концентрация примеси С, строится путем исключения из уравнений Навье-Стокса для несжимаемой жидкости функции давления [2], и может быть записана в декартовых координатах для случая, когда можно пренебречь учетом двойного электрического слоя, следующим образом:
дюг дю, дюг дюг ди ди ди (д 2®. д2 ®г д 2о>, dt дх ду dz дх у ду dz	дх2	ду2 dz2 J
+•„ [	де _ - асї,(1)

V g ду g oz) g ду g dz)
дюг	дюr	дюr	дюr	dv dv	dv	\	s2®.	d2®r	d2®r
—- + и—- + V—- + w—- — ю,	ю ,	ю, — = vl	Z2- +	
dt	dx	dy	dz	dx	y dy	dz	I	dx2	dy2	dz2
+ gPr
gx ST
g 8z
g 8x )
(2)

дю,	дю,	дю,	дю,	dw dw	dw	| д2ю,	д2ю,	d2®,
dt dx	dy dz	dx	y dy	dz dx2 dy2 dz2
+ gPrf— -——V gPcf,	(3)
I g dx g ду ) I g dx g dy )

	52Ух t 82Ух
	8 2V х „ 4	= ®х ,
8'	1
	
	(4)

	8Уу 8Уу
	д 2V у
+ я 2 ~®У ’ ду
	
	(5)

	82уг | 82уг 8х2 8у2
	8 2V г + ®. , бу
	
	(6)

	_ (дТ	дТ	дТ	8Т J 1
р J St	дх	ду	dz JI
	2 8 2Т 8 2Т 8 2Т J ^дх2+ дуdz2 ] ’
	(7)

	8С 8С
	Н и	1- V
St дх
	8С 8С Ґ 8 2С — + w— = D\ —- ду	dz J 8х2
	' д2С д2С\
"н	Г н	1 ’
8у 2	dz2 )
	(8)

	8w Sv	ди 8w
® Г =	, ® V =	, ® 7 =
ду	8z	у dz дх
	ду _8у.
8х 8у
	(9)

	и	х	.
ду dz	dz дх
	SV у Sv г w = — 	— ,
8х	8у
	(Ю)


Здесь t - время; и, v, w - составляющие скорости в проекции на оси
декартовых координат х,у, z; &х, &у,	- составляющие вектора вихря; ух, -
составляющие векторного потенциала; Т - температура; С - объемная концентрация компонента; v - коэффициент кинематической вязкости; g - ускорение, создаваемое массовыми силами; [V — температурный коэффициент объемного расширения; рс - диффузионный коэффициент объемного расширения; р - плотность среды; ср - удельная теплоемкость; /. коэффициент теплопроводности; D - коэффициент диффузии.
Граничные условия для уравнений (1) - (6) записываются в предположении, что компоненты вектора скорости на внешних границах области являются заданными. Для области с непроницаемыми стенками значения компонентов вектора скорости в соответствии с условием прилипания равны нулю. Следует отметить, что процессы тепло - и массообмена при естественной конвекции в ограниченном пространстве встречаются во многих технических и биологических объектах. Также они происходят в помещениях зданий.
Для случая, когда область имеет форму параллелепипеда и его стенки непроницаемы, граничные условия уравнений (1) - (3) для составляющих вектора вихря с учетом уравнений (9) записываются следующим образом
= О,
®х
8у
-to.
dz
дх
dv
ОХ
при х = const;
(11)


ди

(12)
= 0 ,
=-
~~
8у
при у = const;

ди
Л OZ
-
0,
при z = const.
(13)

о.>
® у
® у

	Граничные условия на стенке с прилипанием уравнений (4) - (6) для

	составляющих векторного потенциала с учетом выражений (10) следующий вид [6]
	имеют

	= 0, &
	Ч у = 0 ,
	Ч х = 0
	при х = const;
	(14)

	V х = °,
	8ч у п
—- = 0 , 8У
	Ч х = 0
	при у = const ;
	(15)

	V х = 0,
	Ч у = о,
	д-^^- = 0, dz
	при z = const.
	(16)


Для уравнений переноса энергии (7) и переноса массы вещества (8) могут быть заданы граничные условия тепло и массообмена первого, второго и третьего рода.
Для решения системы уравнений (1) - (10) при заданных краевых условиях предлагается явный разностный метод, излагаемый ниже для случая, когда область, в которой протекают процессы трехмерного течения и тапломассообмена, представляет собой прямоугольный параллелепипед с непроницаемыми стенками. Для аппроксимации дифференциальных уравнений разностными вводится пространственно-временная сетка с координатами
xi = ihx , У J = jhy , zm = mhz , tn = nl ,	(17)
где hx, hy , h2 — шаги по координатам x , у, z ; I - шаг по времени; і = 0, =	m=0,l,...,M;n = 0,1,... ;Xl=X,yj=Y,zM=Z.
Численное решение уравнений переноса вихря, энергии и массы проводится на основе трехслойной пересчетной разностной схемы Никитенко, предложенной в [7]. При этом дифференциальному уравнению переноса ставятся в соответствие два разностных уравнения и искомая функция на каждом временном шаге вычисляется в двух приближениях. Разностное уравнение для первого приближения является двухслойным. Оно аппроксимирует неполное уравнение переноса, в котором сохраняются только конвективные члены и временная производная. Для второго приближения используется трехслойное разностное уравнение, строящееся путем аппроксимации всех членов исходного дифференциального уравнения. В соответствии с этой схемой разностные уравнения, аппроксимирующие уравнение (1) переноса вихря в проекции на ось х иуравнение переноса энергии (7) записываются следующим образом:
8,+ и8го, + v8юr + w8„о. -шг8,и -а„8,и -т8м =
1л	л л	ул	Z л	л л	У У	z z
= gPrI —SyT - ' S2T | + gpcI	8yC-'8zc|,	(18)
IS S )	(, g g )
(1 +	)8,- ^..8.<4 + u8x&x + v8уax + w8zax - ax8xu - &y8yu - ®zy.u =
= v(8„- +3„шх +8zzйх)+ gPr(8yT -'SZT\gPcf 5yc -S2C1 ; (19)
V g g ) V g g )
8tT + и8 ХТ + v8 уТ + w8 ZT = 0,	(20)
(1 +	18,Т-ат8,Т"-' + «8ХТ + v8уТ + w8ZT = (з„+ 8„<5X + 8яЙ,>/(cpp) . (21)
В разностных уравнениях (18) - (21) искомые сеточные функции ф",-.„ и ф",(ф = ОУХ, oyy,oyz, и, v, w, Т, О для точки (Xi,yj, zm, tn) записаны для простоты без индексов, т.е. ф = ф",-...., ф = ф”^; (ф = m, Т, С) - весовой параметр разностного уравнения;	8. фй = (ф^ - фйу. m )/I;	8,ф = (ф”7+... - ф”7> ) / і;
8. Ф = (Ф-+1.Л» -Ф”-1,7,т)/(2йх); У.о (Ф”+1,7,„ + Ф!\7.„ -2ф»Л„)/h2.	Разностные
аппроксимации дифференциальных уравнений (2), (3) и (8) записываются аналогично (18), (19) и (20), (21).
Необходимые условия устойчивости решения разностных уравнений находятся при помощи метода условного задания некоторых искомых функций [8]. При П” = 0, когда уравнения (19), (21) являются двухслойными [7], шаг по времени I должен удовлетворять условие
I < min{/r, 7Ш, 1Т, 1С},	(22)
( „	„	\-1	Г / о о оїі-l
где ly =	/hx + v"m Ihy +	Ihz) ,	. [2v(l/h2 +1/h2 +1/k2 If ,
lT = [2(1/h2 +1/h2 +1/h2\	.lD = [2в(1/h2 +1/h2 +1/h2z )]■*.
Если lv > ly = min{/m, lT, lc }, тогда можно найти такие значения параметров П”, при которых шаг по времени I = lv. Согласно условию устойчивости явной трехслойной разностной схемы это достигается, когда [7] 1= (lv / /ф -1) / 2 . Если ly /lv< 1, тогда следует положить = 0.
Уравнения (4) - (6) для составляющих векторного потенциала ух, фу , vz решается на каждом временном слое методом установления с использованием трехслойной явной разностной схемы [7]. На сетке, отличающейся от (17) тем, что вместо реального времени t, вводится дискретная переменная т к = kl^, к = 0, 1, 2, ...,	> 0 , разностная аппроксимация уравнения (4) может быть
записана в виде
(1 + Иф )8,у,кх1]т -Иф8,=8
ххVxijm	Vxijm + ^xxVxijm ^xijm (23)
где £ф, - весовой параметр, О,, - 0. Значения весового параметра £ф, находятся в соответствии с условиями устойчивости уравнения (23) [7]:
£7.7. / /ф0 -1)/2 при /ф> /ф0; <>.. Опри < /ф0,	(24)
причем величина временного шага /^0 = [2(1/hz +1/h; +1/hL_ )]"' отвечает условию устойчивости явной двухслойной разностной схемы, являющейся предельным случаем трехслойной схемы при су, = 0. Процесс установления решения уравнения (23) считается завершенным при удовлетворении условия
-V)//vo *Д,	(25)
І j W
где Д - малое положительное число. В этом случае полагается, что =	■ В качестве начального приближения, отвечающего значению к = О,
принимается	. Численные эксперименты показали, что максимальная
скорость установления решения уравнения (23) достигается при значениях = 2 + 2,5, которым отвечает увеличение временного шага по сравнению с максимальным для явной двухслойной схемы в 5'6 раз.
Значения vх в граничных узловых точках на слое к + 1, в соответствии с условиями (14) - (16), определяются по следующим сеточным уравнениям
= (4^Х>	)/3’	= (VXu,» )/3,	(26)
УХ.» = ЧХ,т - <-,о = УХЖ = 0 ■	(27)
Уравнения (26) аппроксимируют первое из условий в (14) с погрешностью порядка O(hz).
Сеточные функции И дЦ,,1,, на временном слое п +1 определяются с использованием уравнений, которые строятся так же, как и уравнения (23) - (27) при нахождении .
Температура Т и концентрация С в граничных узловых точках при условиях теплообмена первого рода считаются заданными. При граничных условиях второго и третьего рода входящие в эти условия производные дТ / 7; и дС / 7;, где с — нормаль к граничной поверхности, заменяются односторонними разностными отношениями, которые, как и разностные уравнения (26), (27), имеют погрешность порядка O(hz). Составляющие вектора скорости , г",,,1 и іг",,,1 во внутренних узловых точках области определяются по разностным уравнениям, аппроксимирующим выражения (10):
„»+1 _ я	я ...»+г	,л+1 я „л+1	я „л+1	я „л+1	я „л+1	/эех
Uijm ~ °у У zijm	Уі' yijm	Vijm ~ °zV xijm	У I' zijm >	Wijm ~ M yijm	У I' xijm ■	(У)
Цикл вычислений на слое n +1 завершается определением значений составляющих векторного потенциала у х, у у , у 2 в граничных узловых точках по разностным уравнениям, аппроксимирующим условия (11) - (13) с погрешностью порядка (М1, + к1. + к].'). Для узловых точек, расположенных на плоскостях х = 0 и х = X эти уравнения имеют следующий вид:

®Х> = 0;®Х> =—12hx J
,	4г"'1 - v”+1
n+1	2,J,m	3,J,m	п /оп\
mz,i,7> =	~Y,—— при х = О, (29)
£FLV
л,,,»+1	_ «>и+1
„и+1 л . „и+1	WI-2,j,m
,j,m - U Л»у,I,j,rn~
ЮЧ,j,m -
4v”+1 T'7 .
- ‘ - ■
w	1 -2, j.m
2>)
при

x = X.	(ЗО)
Для точек, расположенных на плоскостях у = 0, у = Y и z = 0, z = Z, разностные аппроксимации граничных условий (12) - (13) записываются аналогично выражениям (29), (30).
Процессы тепло- и массообмена при естественной конвекции в ограниченном пространстве встречаются во многих технических и биологических объектах, в светлопрозрачных ограждающих конструкциях. Для практических приложений представляет интерес случай переноса энергии и массы через щелевые прослойки, заполненные жидкостью или газом. Температуры и концентрации примеси Т\, Q вблизи левой стенки щели и Тг,
С2 вблизи правой стенки принимаются обычно постоянными. Для характеристики процессов тепло- и массообмена используются следующие критерии подобия:	Gr = gftX’’(72 - 7))/ v2	число Грасгофа;
Gro = gPсх3(Q - Cj)/v2- диффузионное число Грасгофа; Pr = vpcp / Z - число
Прандтля; Sc = v / D - число Шмидта; Ra = GrPr - число Релея. За характерный размер принимается толщина прослойки X .
Локальные значения теплового потока изменяются по поверхности стенки щели вследствие сложной вихревой структуры течения в прослойке. 13 связи с этим при обобщении экспериментальных или числено найденных данных вводится эквивалентный коэффициент теплопроводности Z.,,. Отношение у eq!> гДе - коэффициент теплопроводности среды, заполняющей прослойку, характеризует влияние конвекции на перенос энергии через щель при различных Ra.
Предотвращение развития неустойчивости численного решения при значениях числа Релея, которым отвечает переходной или турбулентный режим течения, может быть достигнуто путем наложения на временные производные от искомых функций в некоторых узловых точках области, в которых эти скорости превышают допустимые значения, ограничений следующего вида:
8W л dW dW
	= А^	1	
dt dt dt
если
dW
dt
>AW,
W = u, v, w,	rny rnz, T, C.
(31)

Здесь Aw - положительная величина, которая может быть выбрана исходя из требования минимизации числа узловых точек Bw, в которых на данном временномслое реализуетсяусловие (31).
Степень влияния ограничения (31) на результаты решения характеризуется величиной bw определяемой как отношение числа узловых точек Bw, в которых производится ограничение скорости роста функции W на временном слое п, к общему числу узлов разностной сетки - В = I х J х М. Численные эксперименты показали, что Aw слабо зависит от числа узловых точек пространственной сетки, геометрических параметров области течения и числа Ra. При расчете турбулентного переноса относительное число bw = Bw / В узлов, в которых осуществляется коррекция скорости, вначале достаточно быстро возрастает, затем, достигнув максимума, монотонно снижается. По мере уменьшения шагов разностной сетки все более мелкомасштабные вихри оказываются учтенными при численном решении уравнений (1) - (10) и это приводит к снижению величины bw. Для достижения лучшего согласования результатов расчета эквивалентного коэффициент теплопроводности Xeq с опытными данными, целесообразно в режиме развитой ламинарной конвекции (при Ra < Ю5) ограничивать скорости изменения функции и, v, w, Т, С , а при переходном и турбулентномрежиме - функции сосоу, соz, Т, С.
Результаты конкретных расчетов иллюстрируются на примере тепловой конвекции в замкнутой щелевой области в виде параллелепипеда. Две вертикальные стенки области, образующие щель, имеют постоянные температуры Тх на левой и Т2 на правой стенке. Остальные грани теплоизолированы.
На рис.1 приведены данные о зависимости величины ак от числа Рэлея, которые получены в результате численного решения (они представлены точками •) и путем обобщения уравнениями подобия [9] экспериментальных результатов Бояринцева, Муль-Рейера, Девиса, Бекмана, Крауссольда и др. для вертикальных и горизонтальных плоских щелей, кольцевых и сферических слоев, заполненных газом или капельной жидкостью. При определении чисел подобия принимается толщина прослойки, а за определяющую температуру — величина Т = (Тх + Т2)/2 . При малых значениях числа Релея (Ra<Ю3) величина ак « 1, т.е. влияние конвекции практически не проявляется. ПриЮ3 <ак < 106 величина ак = O,lO5Ra0,3 (линия 1) и при 106 <ак < 1О10 величина ак = 0,40Ra0,2 (линия 2). Михеевым для всей области значений Ra > 103 предлагается использовать зависимость ак = 0,18 Ra0,25 (линия 3).

Рис.1- Эквивалентная теплопроводность ак слоя жидкости в зависимости от числа Рэлея
Из рис. 1 следует, что результаты численного моделирования течения и теплообмена достаточно хорошо согласуются с экспериментальными данными при различных режимах движения жидкости.
Выводы. 1. Представлена система уравнений, описывающая трехмерное течение неоднородной жидкости и тепломассоперенос в приближении Буссинеска в переменных - вектор вихря т, вектор потенциала у, температура Т, концентрация примеси С.
2. Для решения системы уравнений при заданных краевых условиях предлагается метод численного моделирования динамики трехмерного течения и тепломассообмена вязкой жидкости при ламинарном, переходном и турбулентном режимах движения. Для расчета переходного и турбулентного режимов движения жидкости привлекается метод торможения скоростей изменения искомых функций. Численное решение уравнений переноса вихря, энергии и массы проводится на основе трехслойной пересчетной разностной схемы Никитенко.
3. Результаты расчетов иллюстрируются на примере тепловой конвекции в замкнутой щелевой области в виде параллелепипеда. Результаты численного моделирования течения и теплообмена достаточно хорошо согласуются с экспериментальными данными при различных режимах движения жидкости.
Литература
1. Пасконов В.М., Полежаев В.И., Чудов Л.А. Численное моделирование процессов тепло-и массопереноса.-Москва: Наука,-1984.-288с.
2. Никитенко Н. И., Нальчик Ю. Н., Сороковая Ю.Н. Метод канонических элементов для моделирования гидродинамики и тепломассообмена в областях произвольной формы. // Инженерно - физический журнал. -2002. -Т. 75, -№ 6. - С. 74-80.
3. Никитенко Н И., Снежкин Ю.Ф., Сороковая Ю.Н., Нальчик Ю. Н. Численный метод моделирования тепло- и массообмена при различных режимах течения в канале с проницаемыми стенками. Инженерно - физический журнал. - 2006. -Т. 79, -№ 3. -С. 91-101.
4. Никитенко Н И., Снежкин Ю.Ф., Сороковая Ю.Н., Нальчик Ю. Н Метод расчета теплообмена при различных режимах течения вязкого сжимаемого газа. Промышленная теплотехника. -2007. -Т. 29, - № 5. -С. 17-23.
5. Калинин Э.К., Лобанов НЕ. Проблемы исследования теплообменных процессов при течении однофазных сред на этапе успешного развития численного моделирования. Тезисы докладов и сообщений VI Минского международного форума по тепло- и массообмену. -2008. - Т.1. -С. 101-103.
6. Роуч П. Вычислительная гидродинамика.-Москва: Мир, -1980,- 616 с.
7. Никитенко НИ. Сопряженные и обратные задачи тепломассопереноса. -Киев: Наукова думка, -1988. -240 с.
8. Никитенко НИ. Теория тепломассопереноса. -Киев: Наукова думка, - 1983.-352 с.
9. Исаченко В.П., Осипова В.А., Сукомел А.С. Теплопередача.-Москва: Энергия, -1975,- 488 с.
Метод чисельного моделювання тривимірного тепло- і масообміну при різних режимах течії.
|М.І.НІкітенко|, Ю.М.Кольчик, Н.М.Сорокова.
Наведено метод чисельного моделювання динаміки тривимірної течії і тепломасообміну в'язкої рідини на базі рівнянь Нав'є-Стокса при ламінарному, перехідному і турбулентному режимах руху. Результати чисельного моделювання достатньо добре узгоджуються з дослідними даними.
Ключові слова: Тримірна течія і тепломасообмін, перехідний і турбулентний режими руху, метод гальмування швидкостей зміни шуканих функцій
Method of numerical modeling of three-dimensional heat and mass transferfor different flow modes.
|N.Nikitenko|,Y.Kolchik,N.Sorokovaya
Presented a method for the numerical simulation of the dynamics of three-dimensional flow and heat and mass transfer of a viscous fluid on the basis of the Navier-Stokes equations for laminar, transitional and turbulent regimes of movement. The results of the numerical simulation is sufficiently good agreement with the experimental data.
Keywords: Three-dimensionalflow andheat and mass transfer, transitionand turbulentmotion, the method ofbrakingrates of change ofthe desiredfunctions
Вентиляція, освітлення та теплогазопостачання. Вип. 18, 2015

Надійшла до редакції 22.04.2015 p.
#

|Н. И.Никитенко * 1|, Ю.Н.Кольчик 2,Н. Н.Сороковая 3
1д.т.н., професор, Институт технической теплофизики НАН Украины
2 К.Т.Н., доцент, Киевский национальный университет строительства и архитектуры, yulia@orblink.kiev.ua
Зк.т.н., С.Н.С. Институт технической теплофизики НАН Украины, n.sorokova@yandex.ua
Излагается метод численного моделированиядинамики трехмерного течения и тепломассообмена вязкой жидкости на базе уравнений Навъе-Стокса при ламинарном, переходном и турбулентном режимах движения. Результаты численного моделирования достаточно хорошо согласуются с опытными данными.
Ключевые слова: Трехмерное течение и тепломассообмен, переходной и турбулентный режимы движения, метод торможения скоростей изменения искомых функций.
Вступление. Математические модели для описания тепло - и массообменных процессов при турбулентном режиме конвекции обычно базируются на идее Рейнольдса об усреднении уравнений Навье - Стокса (RANS). Согласно этой идеи мгновенные значения скорости, давления, плотности и температуры представляются суммами их средних и пульсационных значений. В результате исходные уравнения Навье - Стокса с использованием некоторых дополнительных условий осреднения искомых функций во времени, преобразуются в уравнения относительно осредненных значений искомых функций. В осредненных уравнениях переноса пульсации скорости, температуры и концентрации компонента вызывают появление членов, которые трактуются как турбулентное трение, турбулентная теплопроводность и турбулентная диффузия.
Установление взаимосвязи между характеристиками осредненного и пульсационного переноса является достаточно сложной задачей. Ее решение сопряжено с необходимостью использования большого объема эмпирической информации, что приводит к снижению достоверности и универсальности получаемых результатов расчетов турбулентной конвекции. В последние годы все чаще высказывается заключение, что математическое моделирование процессов турбулентного течения и тепломассообмена на базе осредненных уравнений переноса импульса, массы и энергии является малоперспективным, и более предпочтительными представляются методы математического моделирования на основе исходных дифференциальных уравнений движения, неразрывности, тепло- и массопереноса [1-5].
image1.jpeg