CC BY
34
Вестник Кузбасского государственного технического университета. 2015. № 6 125
ПРИКЛАДНАЯ МАТЕМАТИКА
УДК 536.24: 541.2
ТЕПЛО-МАССООБМЕН В НЕИДЕАЛЬНОМ АКСИАЛЬНОМ ПАРОГАЗОВОМ ПОТОКЕ
Анисимов Михаил Прокопьевич1 ,
докт. физ-мат. наук, проф., зав. лаб.. Email: anisimovmp@mail.ru
Тукубаев Нурлан Толгатович1, ведущий программист. Email: nurlntt@mail.ru Беликов Дмитрий Викторович2, канд. техн. наук, начальник отдела. Email: dimitrv.belikov@yandex.ru
Овчинникова Тамара Эдуардовна3,4, канд. физ.-мат. наук, ст. н. с., доцент. Email: teonew@iwep.nsc.ru
Пинаев Виктор Алексеевич5, канд. физ.-мат. наук, доцент. Email: pinaev@kemcitv.ru
1 Конструкторско-технологический институт научного приборостроения СО РАН. 630058, Новосибирск, ул. Русская, 41
2 Научно-исследовательский институт - Республиканский исследовательский научноконсультационный центр экспертизы, г. Москва, Антонова-Овсеенко ул., 13, стр. 1.
3Институт водных и экологических проблем СО РАН. Россия, 656038, г. Барнаул, ул. Молодежная, 1.
4 Новосибирский государственный университет, г. Новосибирск, ул. Пирогова, д. 2
5 Кемеровский институт (филиала) Российского экономического университета им. Г.В. Плеханова, 650992, Россия, г. Кемерово, пр. Кузнецкий, 39.
Аннотация. В работе приводится описание численной модели расчета тепломассообмена в аксиальном реальном парогазовом потоке. Модель предназначена для интервала давлений от низких (несколько Торр) до нескольких сотен бар и температурного интервала от эмпирически легко достижимых, порядка -100 °С, до критических условий, достигающих нескольких сотен градусов Цельсия. Модель предназначена для описания условий зародышеобразования в поточной диффузионной камере высокого давления.
Ключевые слова: нуклеация, тепло-массообмен, поточная диффузионная камера
Перспективным подходом для исследования нуклеации в паро-газовой среде является применение ламинарного осесимметричного парогазового потока, реализованного в методе, названнном «поточная диффузионная камера» [1]. Этот метод позволяет определить зависимость скорости нуклеации от различных параметров процесса.
В данной публикации рассматривается задача тепломассообмена в поточной диффузионной камере высокого давления. Это вызвано необходимостью измерений скоростей зародышеобразования при повышенных давлениях (до 100-200 бар) в системах, которые, например, находятся влизи критической линии для бинарной паро-газовой системы. Математическая модель тепломассообмена строится на базе стационарной системы уравнений Навье - Стокса, к которой добавляется уравнение состояния неидеального газа-носителя. Пар является малой добавкой к давлению газа-носителя.
Решение этой задачи для низких давлений (до
5 атмосфер) в приближении идеального газа было представлено нами ранее [2].
Для экспериментов при повышенных давлениях, например, в окрестности критических условий возникает необходимость создания поточной диффузионной камеры для давлений порядка 100 бар и более, что неизбежно приводит к необходимости использовать реальное уравнение состояния парогазового потока и реальные коэффициенты переноса.
В работе приводится оригинальное описание разработанной и реализованной численной модели расчета процессов тепломассообмена в аксиальном неидеальном парогазовом потоке, для интервала давлений от долей до нескольких сотен бар и температурного интервала от эмпирически достижимых, порядка -100 °С и ниже, до критических условий, порядка нескольких сотен градусов Цельсия. Модель предназначена для описания условий зародышеобразования в поточной диффузионной камере высокого давления.
126 М. П. Анисимов, Н. Т. Тукубаев, Д. В. Беликов, Т. Э. Овчинникова, В. А. Пинаев
Математическая модель
Стационарное движение осесимметричного ламинарного потока парогазовой смеси описывается системой уравнений Навье - Стокса в приближении пограничного слоя.
Эти уравнения в системе цилиндрических координат имеют следующий вид [3]:
2-(г-р-и) = 0; (1)
дх
р-и
ди
дР }__д_
дх г дг
Г • TJ
ди
^ дТ дР
PwCp —= м —+ 77
дх дх
/ Я Л ди 2 1 д + ,дТ~ ГА
г дг дг
дС
р-и— = ~ —
\_д_ дх г дг P = pR-TZ(pj\,
,8С г-р-а---
дг
(2)
(3)
(4)
(5)
Здесь х - координата по оси цилиндрической трубы, 0 < х < L ;
г - координата по радиусу, 0 < г < R ; R-радиус трубы
Т (х, г) - абсолютная температура;
С (х,г) - концентрация пара исследуемого вещества;
Р (х,г) - давление ;
Ср - изобарная удельная теплоемкость газа;
р{Т,Р) - плотность газа;
г/{т, Р), Я(Г,Р), d(T,P) - коэффициенты вязкости, теплопроводности газа и диффузии рабочего вещества в газе-носителе, соответственно;
R
g
- газовая постоянная, Ru- уни-
версальная газовая постоянная, М - молярная масса газа;
Z (/Э, Т) - фактор
сжимаемости газа.
Искомыми функциями являются и(х,г) - скорость потока; Р(х,г) - давление; р{х,г) -
плотность; Т(х,г) - температура; С(х,г) - концентрация паров вещества.
Уравнения (1) - (5) имеют следующий физический смысл:
(1) - уравнение неразрывности в стационар-
ном случае при условии
= 0;
(2) - уравнение движения (закон сохранения импульса);
(3) - закон сохранения энергии;
(4) - уравнение переноса примеси;
(5) - уравнение состояния газа.
Для однозначного решения системы уравнений (1)-(5) принимаются следующие граничные условия:
на оси симметрии:
|^-(х,о)=о, ^—■(*,())=0, ^(х,0)=0 ;
дг дг дг
на стенке:
u(x,R) = 0, T(x,R)=Tw(x),
C(x,R)= Psal(Tw(x))/(Kh-TJx))
на входе:
Г(0,г)=Г0(г), и(0,г) = и0(г),
Р{0.г)=Р0. С(0.г)= РАЪЖКь-Тш)
Здесь
Tw(x) - заданная температура стенки трубы, То(г) - начальный профиль температуры потока,
Ро -давление при входе в трубу; ио(г) - начальный профиль скорости потока при входе в трубу;
Кь -константа Больцмана..
Предполагается, что распределение скорости на входе, ш(г\ имеет параболический профиль, согласованный с массовым расходом газа Q :
А
Q = 2тг^гр udr; w0(r) =
2 Q
и R2
'.-4'
* J
Давление насыщенных паров исследуемого вещества РчШ{т')определяется по формуле Антуана :
W-'-tTc’
где А, В, С - эмпирические коэффициенты для данного вещества.
Для фактора сжимаемости неидеального газа используется следующая эмпирическая аппроксимация [4]:
z{co,r) = 1 + а1 (ет -1 - r\cotr - со)2 (4со- cotr )со --а2со(е~т -1)-аъ со(е~3т - l)- а4 сот +
+ а6 (евт -6r\cotr - cof (3со- cotr )со2 +
2
+ а
со
7l-Zcco + a\l-Zcco)2
со
+ а о
со
- + «10-
со
o-z^r “'VW
Здесь
СО — р / рс — приведенная плотность газа:
Т — Тс /Т - приведенная обратная температура;
-Лгу
Ю/г =■
Рс
; ри., рс - значение плотности
газа в тройной и критической точках;
Тс - критические значения температуры газа; Zс коэффициент сжимаемости газа в критической точке;
Вестник Кузбасского государственного технического университета. 2015. № 6 127
Рис. 1.Изотермы р(р) для С02.: 1 : Тк=Тс+50; 2 : Тк=Тс; 3 : Тк=Тс-50
а{,а2 , ...,<з10 - эмпирические коэффициенты для конкретного газа.
Для вычисления плотности р неидеального газа необходимо для данных температуры Тк и давления Рк решить нелинейное относительно р уравнение
F{p) = Pi-RsTkpZ{jk,p) = 0. (6)
В данной работе решение находится методом бисекции [5]. Для того чтобы применить стандартную процедуру поиска корня, необходимо локализовать искомое значение на таком отрезке
[Л>Л]> что значения г(л) и F(p2) имеют разные знаки. Графики изотерм
Р(р) = RgTkpZ(p,Tk) при трех значениях Тк приведены на рис. 1.
Из графиков видно, что при Тк > Тс уравнение (6) давление является однозначной функцией плотности. В этом случае задаем некоторый достаточно малый отрезок Ар и, последовательно перебирая отрезки с концами рх = / Ар, р2= рх+ Ар, / = 0,1,..., находим первый же отрезок [рх р2 ], для которого Р(рг) • Р{р2 ) < 0 .
Далее применяется стандартная процедура бисекции.
При Тк < Тс решение уравнения (6) не единственно. В этом случае нужно искать наименьший из корней. Сначала находим значение р0, такое,
что на отрезке [0, р0 ] Р(р) > 0, а затем уходит в
отрицательную область значений (процедура аналогична предыдущему поиску). Затем на этом отрезке находим точку локального максимума р{.
Физически плотность газа при Тк < Тс определена лишь для давлений в интервале [0,Р(рх)], а искомое значение локализовано на отрезке
[0, рх ], на котором решение единственно и может быть найдено применением стандартной процедуры.
Динамическая вязкость неидеального газа определялась по методике, предложенной в [6]. Для заданного давления и температуры описанным выше способом определяется плотность р . По полученному значению плотности и температуре Т вычисляется плотность энергии взаимодействия:
у(р,Т) = Аи(ы,т)-р, где AU = U,S(T)-U(T,P) - разность внутренней энергия идеального газа и внутренней энергии вещества в заданном состоянии:
AU = R-TC■ {-ар - - of а-
- а2сое~т - 3 • аъсоё~Ъ т + аАсо --Ъ-а5(е3т -1 \cotr - со)5со +
+ 6 • a6e6 r(cotr - cof со2};
Далее рассчитывается вязкость в разреженном газе [.мкПа-с], как функция только температуры:
арМТ
7„(Т) = -
ехр
IЬ,
/=0
In
т
Здесь М - молекулярная масса вещества, а, bi
и 8 - эмпирические коэффициенты для каждого вещества.
Затем вычисляются вязкость плотного газа по Энскогу и избыточная вязкость.
Вязкость плотного газа по Энскогу:
= 1 + 0.8 ,yg(y)+0.7614y2 g2(y). g(y)
gO)=l+0.625>^0.287y2+0.1103y3+0.0386/+
128 М. П. Анисимов, Н. Т. Тукубаев, Д. В. Беликов, Т. Э. Овчинникова, В. А. Пинаев
+0.0127у5;
Избыточная вязкость
пЛТ’Р) = А'УЫ Р
а-
Уъ
Уо-У.
Здесь yQ=Hy^, A,a,H0,V0 ~ индивидуальные эмпирические константы газа; Н0 - энер-
гия сублимации, VQ - гипотетический удельный
объём переохлажденной до абсолютного нуля жидкости.
Вязкость {Па-с) равна сумме вязкости по Эн-скогу и избыточной вязкости:
п{т,р)=[ча,Ат,р)+пЛт,р)]-10-*;
Для вычисления коэффициента диффузии в бинарной системе вещество - газ применена эмпирическая корреляция Фуллера, Шеттлера и Гиддингса [7]:
d(T р)= оооизг'^Км, +м2)/(м,м2)]'
где Р - общее давление, выраженное в барах; М/, М2 - молярные массы исследуемого вещества и соответствующего газа-носителя (г/моль); £Н(2>),- молекулярные диффузионные
объемы компонент.
Теплопроводность газа Л(Т,Р) по эмпирической формуле [8]:
Л(Г,Р) =
= т'15 0к0р + Л,) + Р{к2Р + к3) + кл
где к0,к],...,к4 - эмпирические константы,
приводимые в справочной литературе для каждого газа.
Численный метод
Для численного решения задачи в расчетной области строится сетка узлов в осевом направлении х и радиальном г. По оси х сетка равномерная с постоянным шагом Ах (L - длина трубки):
Х; = iАх, Ах = L/(N — 1), / = 0,1,..., А — 1.
По оси г сетка строится с таким расчетом, чтобы на каждом шаге было постоянным приращение площади сечения:
= Г; У = -1;
М -1
В уравнения (2)-(4) входит дифференциальный оператор вида ^ \г(р - Б (2) в качестве дг дг
(р выступает вязкость Т)(х,г), а в качестве у/ -скорость и(х, г) .В (3) ^^-теплопроводность Л(х,г) и температура Г(х,г), а в (4) (р , у/ -плотность р(х,г) и концентрация С(х, г) .
Этот дифференциальный оператор аппроксимируется следующим образом:
{r<p)i
+ 1/2
П+1-П
¥i+1
М,+1/2 ! Mj-
1/2
L П + \~П П-Г;_ 1
Vi
>/{rj
7+1/2 _r/-l/2
V/+l + Vi
где Vj~у 2 = —^----- ~ значение переменной в
середине интервала сетки.
Для нахождения скорости и и давления Р из (1) - (2) используется метод прогонки с параметром [9].
Для этого разностная аппроксимация уравнения (2) на (/+ 1)-м шаге по х определяется следующим образом:
гК {r)^f^ + rP'r'=D\uM), (8)
Ах
п а а
где Dr = —г л—. ' дг дг
Профиль скорости на (/+1) -м шаге по х представим в виде
им =и‘+zM+Рр'ум, (9)
i+1 /+1 ,
где Z и у - параметры (все эти переменные являются векторами из М элементов).
Подставляя (9) в (8) и группируя члены, содержащие и не содержащие Р , приходим к системе уравнений:
^(/>^ = аФ'Ч|)+А-И; (Ю)
L\ JC
rb0(r)^— = Dr (/+1)- Г. (11)
Ах
Граничные условия для параметров z‘+1 и / +1
у те же, что и для скорости и :
z'+l(Л) = У+|(Л) = О, ^(0) = ^21.(0) = 0.
дг дг
Для определения Р‘х+Х воспользуемся формулой (9) и массовым расходом, который не зависит от длины трубки, х .
В результате получим соотношение
R
J rp‘z,+]dr
if1 =-21---------, рм=р‘+руАх. (12)
| rplyl+ldr о
После того, как в (/+1) -м сечении найден профиль скорости, для этого сечения вычисляются
Вестник Кузбасского государственного технического университета. 2015. № 6 129
Temperature
Рис.2. Поле температуры
профили температуры и плотности.
Поскольку коэффициенты вязкости, теплопроводности и диффузии зависят от температуры и давления, то для раскрытия нелинейностей используется итерационный процесс, на каждом шагу которого строится аппроксимация уравнений по координате г, а для решения полученной трехдиагональной системы линейных алгебраических уравнений (СЛАУ) используется метод прогонки [10].
Итак, для каждого сечения трубы
/ = 1,..., N — 1; выполняется цикл из трех этапов:
1) итерационный процесс, состоящий из 5 шагов (а)-(е):
a) В узлах (i, j) j = 0,1,..., М — 1; текущего сечения находим значения параметровz\ и уУ
путем решения СЛАУ полученных из (10) и (11), соответственно;
b) Находим давление в этих точках по формуле (12);
c) Находим профиль скорости потока
и* у = 0,1,...,М-1 по формуле (9);
d) Определяем в этих точках температуру из уравнения (4);
e) Вычисляем новое значение плотности газа по уравнению состояния (5).
Если значение плотности отличается от значения, полученного на предыдущей итерации, меньше чем на заданную величину или количество итераций превосходит предельное значение, итерационный процесс заканчивается, иначе процесс возвращается к шагу а).
2) На основе уравнения (4) выполняется рас-
чет /-го профиля концентрации С{ паров вещества.
3) В каждом узле /-го профиля вычисляется пересыщение пара по формуле:
S(x,r)=C(x,r) / С°(х,г),
где С°(х,г)-концентрация насыщенного пара исследуемого вещества при температуре Т(х,г). Таким образом,
S(x,r)=C(x,r) Kh Т(х,г) / (Рш Т{х,г)).
После расчета всех сечений получаем поля
температур 7V , С/ концентраций паров вещества и пересыщений Sj.
На рис. 2 и 3 приведены рассчитанные по этим алгоритмам распределения температуры и пересыщения паров глицерина в конденсирующем устройстве для системы глицерин-двуокись углерода (распределение концентраций пара топологически не отличается от распределения температуры на рис. 2 и здесь не приводится). Расчет выполнен для трубки с радиусом 0.6 см и длиной 500 см, при температуре стенки 40°С, температуре насыщения 110°С, объемном расходе газа 6.0 см3/с и давлении 20 атм. Константы уравнения состояния (5) и критические параметры для диоксида углерода приведены ниже:
«i=3.67103-10'3; 02=1.76986; «3=5.42282-10'2; «4=2.64972;
«5=9. 87467-10-5; «6=4.43953-10 5; «7=4.8212-10-1; «8=8.25279-10'3; «9=2.38006-10'3; «ю=-2.0473-10'4;
Д*=0.188922; Zc=0.274586; Гс=304.1282;
Рс.=0.4676 г /см3; Ри=1.17853 г /см3;
Коэффициенты для расчета теплопроводности
130 М. П. Анисимов, Н. Т. Тукубаев, Д. В. Беликов, Т. Э. Овчинникова, В. А. Пинаев
Suprsalurafen
30—I
Radius г. сл|
Tube length х, cm
Рис.З. Поле пересыщений
диоксида углерода:
ко=АА 1995-10'9; к\= 5.7498-10'7; к2= 9.5477-10'7; къ= 1.551 * 10"4; к$= 4.0693-10'3
Коэффициенты формулы Антуана для вычисления давления насыщенных паров глицерина в мм рт. ст.:
А=\ 1.27243, 5=4480.0, С=0.
Выводы.
В данной публикации приведен подробный алгоритм и результаты численного решения задачи Навье-Стокса для потока паров глицерина в среде углекислого газа в приближении реальных газов. Получены распределения температур и
концентраций пара, а также пересыщений паров глицерина по потоку. Приведенный алгоритм может быть применен для описания условий заро-дышеобразования в поточной диффузионной камере высокого давления, аналогичной ранее описанной в [1] камере низкого давления.
Статья подготовлена по материалам научно-исследовательской работы, выполненной ФГБНУ НИИ РИНКЦЭ по заданию № 2015/Н7 Министерства образования и науки РФ на выполнение работ в рамках государственного задания в сфере научной деятельности (проект № 3246).
СПИСОК ЛИТЕРАТУРЫ 1 2 3 4 5 6 7 8 9 10
1. Anisimov М.Р. and Cherevko A.G. Gas-flow diffusion chamber for vapor nucleation studies. Relations between nucleation rate, critical nucleus size and entropy of transition from a metastable into a stable state // J.Aerosol Sci., 1985. Vol.16, No. 2. P. 97-107.
2. Анисимов М.П., Овчинникова Т.Э., Пинаев В.А.. Тепло-массообмен в аксиальном парогазовом потоке в приближении идеального газа // Вести. КузГТУ, 2013. № 3. С.88-91].
3. Лойцянский Л.Г. Механика жидкости и газа. М.: Наука, 1987.
4. Каплун А.Б., Мешалкин А.Б. Малопараметрическое уравнение состояния для расчета термодинамических свойств веществ в жидком и газообразном состоянии // Журнал физической химии, 2013. Т. 87, №8. С. 1294-1300.
5. Каплун А.Б., Мешалкин А.Б. Единое уравнение для расчета коэффициента вязкости диоксида углерода // Вестник Казанского технологического университета, 2014. Т. 17, № 21. С. 108-113.
6. Волков Е. А. Численные методы. Учеб, пособие для вузов. 2-е изд., испр. М.: Наука, 1987.
7. Poling В.Е., Prausnitz J.M., O’Connel J.P. The properties of gases and liquids. Fifth edition. McGraw-HILL, 2001.
8. Варгафтик Н.Б., Филиппов Л.П., Тарзиманов А., Тоцкий E.P.Справочник по теплопроводности жидкостей и газов. М.:Энергоатомиздат, 1990.
9. Воеводин А.Ф., Гончарова О.Н. Метод расщепления по физическим процессам для расчета задач конвекции // Мат. моделирование, 2001. Т. 13, № 5. С. 90-96.
10. Самарский А.А., Гулин А.В. Численные методы. М.: Наука, 1989.
Поступило в редакцию 19.09.2015
Вестник Кузбасского государственного технического университета. 2015. № 6 131
UDC 536.24: 541.2
HEAT-MASS TRANSFER WITHIN AN AXIAL VAPOR-GAS FLOW IN THE REAL GAS APPROXIMATION
Anisimov Mihail P.1,
Dr. Sc. (Phys.-Math.), professor, head of Department, Email: anisimovmp@mail.ru
Tukubaev Nurlan T.1, lead programmer, Email: nurlntt@mail.ru Belikov Dmitriy V.2, C. Sc. in Engineering, Email: dimitrv.belikov@yandex.ru
OvchinnikovaTamara E.3,4, C. Sc. (Phys.-Math.), Associate Professor, Email: teonew@iwep.nsc.ru
Pinaev Viktor A.5,
C. Sc. (Phys.-Math.), Associate Professor, Email: pinaev@kemcitv.ru
1 Design-technological Institute of scientific instrument making RAS. Russian str., 41, Novosibirsk, 630058, Russia
2 Research Institute-Republican research scientific-consulting Centre for expertise. Antonov-Ovseenko str., 13, pp. 1. Moscow, Russia.
3HoInstitute for water and environmental problems SB RAS. Molodezhnaya str., Barnaul, 1656038, Russia.
4 Novosibirsk State University. Pirogov str, 2. Novosibirsk. Russia.
5G.V. Plekhanov Russian Economic University, Kemerovo Institute (branch), Kuznetskii pr, 39, Kemerovo, 650992, Russia
Abstract. The digital model of the heat-mass transfer problem for the real vapor-gas axial flow is presented. The model covers the pressure interval from low (some Torr) to several hundred bars and the temperature gap from 100 Centigrade below the ice melting point to the critical conditions around some hundred of the Celsius degrees. The model is aimed for the nucleation condition description within the High Pressure Flow Diffusion Chamber.
Keywords: nucleation, heat-mass transfer, flow diffusion chamber
REFERENCES
1. Anisimov M.P. and Cherevko A.G. Gas-flow diffusion chamber for vapor nucleation studies. Relations between nucleation rate, critical nucleus size and entropy of transition from a metastable into a stable state // J.Aerosol Sci., 1985. Vol.16, No. 2. P. 97-107.
2. Anisimov M.P., Ovchinnikova T.Je., Pinaev V.A.. Teplo-massoobmen v aksial'nom parogazovom po-toke v priblizhenii ideal'nogo gaza // Vestn. KuzGTU, 2013. № 3. S.88-91].
3. Lojcjanskij L.G. Mehanika zhidkosti i gaza. M.: Nauka, 1987.
4. Kaplun A.B., Meshalkin A.B. Maloparametricheskoe uravnenie sostojanija dlja rascheta termodina-micheskih svojstv veshhestv v zhidkom i gazoobraznom sostojanii // Zhurnal fizicheskoj himii, 2013. T. 87, № 8. S. 1294-1300.
5. Kaplun A.B., Meshalkin A.B. Edinoe uravnenie dlja rascheta kojefficienta vjazkosti dioksida ug-leroda // Vestnik Kazanskogo tehnologicheskogo universiteta, 2014. T. 17, № 21. S. 108-113.
6. Volkov E. A. Chislennye metody. Ucheb. posobie dlja vuzov. 2-e izd., ispr. M.: Nauka, 1987.
7. Poling B.E., Prausnitz J.M., O’Connel J.P. The properties of gases and liquids. Fifth edition. McGraw-HILL, 2001.
8. Vargaftik N.B., Filippov L.P., Tarzimanov A., Tockij E.R.Spravochnik po teploprovodnosti zhidkostej i gazov. M.:Jenergoatomizdat, 1990.
9. Voevodin A.F., Goncharova O.N. Metod rasshheplenija po fizicheskim processam dlja rascheta zadach konvekcii // Mat. modelirovanie, 2001. T. 13, № 5. S. 90-96.
10. Samarskij A.A., Gulin A.V. Chislennye metody. M.: Nauka, 1989.
Received 19 September 2015
