Единое уравнение для расчета коэффициента вязкости диоксида углерода Текст научной статьи по специальности «Физика»

УДК 546.817+532.13+548.57

А. Б. Каплун, А. Б. Мешалкин ЕДИНОЕ УРАВНЕНИЕ ДЛЯ РАСЧЕТА КОЭФФИЦИЕНТА ВЯЗКОСТИ ДИОКСИДА УГЛЕРОДА

Ключевые слова: вязкость, газ, жидкость, диоксид углерода, уравнение состояния, внутренняя энергия.

На примере диоксида углерода анализируется проблема аналитического описания зависимости коэффициента вязкости жидкостей и газов от параметров состояния. Отмечается, что уравнением, связывающем вязкость с внутренней энергией термодинамической системы, хорошо описывающем коэффициент вязкости в широком диапазоне параметров состояния, можно воспользоваться только при наличии термического и калорического уравнений состояния. Показано, что, совершая различные приближения в разных областях термодинамической поверхности, из этого уравнения в качестве частных случаев можно получить ряд известных уравнений для вязкости жидкости и газа.

Keywords: viscosity, gas, liquid,carbon dioxide, state equation, internal energy.

The problem of analytical description of liquid and gas viscosity vs. the state parameters was analysed by the example of the carbon dioxide. It was noted that the equation, which unites the viscosity and the internal energy of a thermodynamic system, perfectly describes the viscosity in a wide area of state parameters, but it can be used only if the state equation is available. It was shown that for any approach in different areas of thermodynamic surface, a number offamous equations for liquid and gas viscosity can be obtained as a special case of this equation.

Известно, что в строгой кинетической теории [1] коэффициент динамической вязкости может быть представлен в виде суммы | = | ш + |ф, где первый член обусловлен

переносом импульса при "столкновениях" частиц, а второй дает вклад в перенос импульса за счет межмолекулярного взаимодействия; при малых плотностях стремится к вязкости разреженного газа

Однако современное состояние кинетической теории (см., например, [2]) не позволяет надеяться на получение сравнительно простых уравнений для описания коэффициента вязкости с малым количеством подгоночных коэффициентов в достаточно широкой области параметров состояния. В связи с этим для описания вязкости флюидов широко используются модельные теории и эмпирические уравнения.

Ранее нами было показано [3], что

избыточная вязкость в широком интервале параметров состояния описывается уравнениями:

х 0

(

Д-q = Ax exp

л

Хо - x

(1)

Д| = г|(Т, Р) - ; |(Т, Р) = + Д| . (1а) Здесь P, V, 7- давление, удельный объем и температура системы; ц(Р,Т) - вязкость системы в заданном состоянии; Д| = ц(Р, Т) - цеп (Т, р) -избыточная вязкость; |0 (Т) - вязкость разреженного газа; |еп (Т,р)=|0 (Т)Т(Ьр) - вязкость системы твердых сфер по Энскогу:

у) = — (1 + 0.8уд + 0.761у2 д2),

д

д = 1 + 0.625у + 0.287у2 + 0.1103у3 + (2) + 0.0386у 4 + 0.0127у5, где у = Ьр, р=1Ж - плотность вещества; х = ' плотность энергии взаимодействия;

Ди = Ug (T) - U(P,T) - энергия взаимодействия; U g (T) - внутренняя энергия идеального газа;

U (P, T) - внутренняя энергия системы; x 0 =

H О

V

Vo = lim (V (T, P))

T ^0, P ^0

гипотетический объем

переохлажденной до абсолютного нуля жидкости, Н0 = ид (0) - и(0,0) - энтальпия испарения при

абсолютном нуле температуры и давлении, равном нулю, и(0,0) - внутренняя энергия тела при абсолютном нуле температуры и нулевом давлении; А, а, Ь - индивидуальные константы вещества. Как следует из (1), избыточная вязкость Д| является однозначной функцией плотности энергии взаимодействия Д| = f (х=Ди/У) (см., напр., Д| для диоксида углерода, рис. 1).

Описание экспериментальных данных, близкое к оптимальному, достигается при Ь=У0. Для хорошо изученных веществ, таких как инертные газы, диоксид углерода и др., физические параметры У0 и Н,° обычно приводятся в таблицах термодинамических свойств веществ. В этом случае уравнения для вязкости (1) и (1а) содержат всего лишь 2 эмпирические константы - А и а, так же, как и известное уравнение Аррениуса. В том случае, если данные о У0 и Н° отсутствуют, количество эмпирических констант, которые надо найти из данных по вязкости, в уравнениях (1) и (1а) возрастает до четырех - А, а, У0 и х0. Уравнения (1) и (1а) описывают зависимость вязкости жидкости и газа широкого круга веществ в пределах погрешности эксперимента (таблиц) во всем изученном диапазоне параметров состояния [3].

В работе [4] было проведено сопоставление уравнений (1) и (1а) с рядом других известных уравнений для описания вязкости жидкостей и газов (Эйринга, Френкеля, Андраде и др.).

а

250 200 150 100 50

ю =

и х = -

0

100

200

300

Рис. 1 - Избыточная вязкость А^ CO2 по [3] как функция плотности энергии взаимодействия AU/V в области температур 220...1000 К и давлений 0,1...100 МПа

В работе [3] расчет коэффициента вязкости ряда веществ проведен нами с использованием табличных данных по плотности, энтальпии и внутренней энергии при заданных температуре и давлении, что, конечно же, является существенным недостатком, поскольку затрудняет применение предложенных уравнений в практических целях.

Так как в уравнения для расчета вязкости (1) и (1а) входят температура и плотность (удельный объем) вещества, а в экспериментах обычно в качестве независимых параметров состояния используются температура и давление, то для вычисления коэффициента вязкости в широком диапазоне параметров надо иметь единое уравнение состояния газа, жидкости и флюида.

Энергию взаимодействия AU (T ,V ) можно найти с помощью известного соотношения термодинамики, если известен фактор сжимаемости PV

Z = -

RT

AU (T ,V )= Ug (T) - U (T ,V ) = R- T2 КЦ 1 ^, (3)

H

газовая

,0Jp p

постоянная, h

'¡9

Я - универсальная молекулярная масса.

Ранее в работе [5] с помощью предложенных в [6] термического и калорического уравнений состояния реальных газов из (1), (1а) получено уравнение для расчета коэффициента вязкости реальных газов при температурах до 1000 К, давлениях до 100 МПа и плотностях до 1,5 рС.

Рассчитанные значения вязкости в указанных диапазонах параметров состояния согласуются с экспериментальными данными по вязкости диоксида углерода и других газов в пределах погрешности экспериментальных и табличных данных.

В настоящей работе для расчета вязкости в широком диапазоне плотностей от разреженного газа до плотности в тройной точке привлечено разработанное нами единое малопараметрическое термическое уравнение состояния, приведенное в [7]. Это уравнение в приведенных переменных

для фактора сжимаемости р с 1с

Z= PV/RT имеет вид (рС, ТС,, РС - параметры критической точки):

Ж = 1 + а1 (ет -1 - т)(ю^ - ю)2 (4ю - ю^ )ю -

- а2ю(е-т -1)- а3ю(е-3'т -1)- а4ют + (е3т -1 - 3т)ю ^ - ю)4 (бю -ю1г )ю + (е6 т - 6 • т)(сог - ю)3 (3ю - ю г )ю2 +

+ a5

+ Эб

a7 ю + —7— + -

1 - Zc ю

(4)

a8 ю

a0 ю

(l - Zc ю)2 (l - Zc ю)3

(1 - Жс ю)4

Уравнение (4) содержит 10 эмпирических подгоночных коэффициентов. Уравнение (4) описывает фактор сжимаемости флюида "нормальных" веществ до 100 МПа и 1000К при плотностях от нуля до плотности в тройной точке в основном в пределах погрешностей опытных (табличных) данных. Поскольку в расчетные формулы для вязкости (1) и (1а) входит плотность, то из уравнения (4) вначале численным методом рассчитывается плотность в зависимости от температуры и давления, а затем по уравнению (3) рассчитывается энергия взаимодействия ДУ (Т V):

ДУ = Я • Тс • {-а! (ет - 1)ю г - ю)3

ю -a2юе X -

- 3 ■ a3юe-3 X -

3 ■ a5 (е3X-1) -ю)5 ю+ (5)

+ а4ю + 6• а6е6т(юг-ю)4ю2 } Заметим, что непосредственно в расчетную формулу (3) для энергии взаимодействия ДП при расчете вязкости входят только 6 из 10 эмпирических коэффициентов, входящих в уравнение (4) для фактора сжимаемости 2.

Затем, подставив х = ДУ/V по (3), (5) в уравнение для коэффициента избыточной вязкости (1а), получим зависимость вязкости от плотности и температуры, которая в приведенных переменных будет иметь вид:

( ^ \

-л(ю, х) = ^AUpc ю exp

X 0

Xо - AUPcю

(6)

Ис"

В (6), как уже упоминалось, в общем случае входят четыре эмпирических подгоночных коэффициента - А, а, х0, а также V0 для расчета вязкости системы твердых сфер (2) по Энскогу, которые должны быть найдены из экспериментальных данных по вязкости. В тех же

и 0 I/

случаях, когда известны Н0 и V0, число эмпирических коэффициентов, которые надо найти из данных по вязкости, уменьшается до двух - А и а. Это позволяет считать, что при расчете этих коэффициентов можно будет ограничиться достаточно малыми объемами исходных экспериментальных данных по вязкости.

В настоящей работе рассчитан коэффициент вязкости диоксида углерода в интервале параметров состояния: по температурам - до 1000 К, по

0

4

a10 ю

давлениям - до 100 МПа, по плотности - от 0 до плотности в тройной точке. В качестве исходных данных по вязкости для расчета коэффициентов А и а были использованы табличные данные из справочника В.В. Алтунина [8] и обзора [9]. В расчетах использованы значения Н0 и У0, приведенные в [8]. В таблице 1 приведены значения коэффициентов уравнений и критические параметры диоксида углерода.

Таблица 1 - Константы уравнения (1), критические параметры диоксида углерода

а1 3.67103-10"3

а2 1.76986

аз 5.42282-10-2

а4 2.64972

а5 9.87467-10-5

ае 4.43953-10-5

а7 4.82120-10-1

аз 8.25279-10-3

ад 2.38006-10-3

аю -2.04736-10-4

R, кДж/(кг-К) 0.188924

Zc 0.274586

ratr 2.520381

Тс, К 304.1282

Pc, МПа 7.3773

рс;, кг/м3 467.6

100-

10

200

400

600

800 T, K

1000

Рис. 2 - Сравнение табличных значений вязкости (символы) [8, 9] с расчетом по (3) - (6) (линии) для С02 при плотностях от 0 до плотности в тройной точке р^

На рис. 2 приведены рассчитанные по (6) значения коэффициента вязкости в зависимости от температуры и давления на изобарах. Как показал анализ, расхождения между рассчитанными по (6) и табличными значениями вязкости С02 по [8, 9] в основном не превышают погрешностей экспериментальных данных по вязкости разных авторов - порядка (2-5)%, кроме критической области, где коэффициент вязкости зависит, в частности, от способа измерений (рис. 3).

<

о о

500

1000

р, kg/m

Рис. 3 - Отклонение рассчитанных по уравнению (6) значений вязкости диоксида углерода от табулированных данных работы [9]

Покажем, что ряд известных уравнений для расчета коэффициента вязкости можно рассматривать как частные случаи уравнения (6), см. также [4]. Так, для нахождения Ди можно воспользоваться каким-либо хорошо

зарекомендовавшим себя эмпирическим уравнением состояния, например, уравнением состояния Ван дер Ваальса

P + -ar|(v - b) = RT ;

или

Р=ж. - а

V - Ь V2 '

Подставив (7а) в (3), получим дР

и = J

дТ

- P

dv = -- + Ug (Т) v у

(7)

(7а)

(8)

Результат, содержащийся в (8), хорошо известен, (см., например, [1]), тогда

ди = а = ар

V

Подставив это выражение в (1), получим

( \

(9)

Д-q = Ax exp

x 0

а-

v хо - x ,

= Aap2 exp

i a p 2

1--p

x 0

(10)

Обозначив

С = Аа и р = —, будем иметь х0

ДЦм = Ср2 exp

а

1 -Р2р2

^vw (Т, p) = ^en + Ср2 exp

а

22

1 -Р2р

(11)

(12)

Таким образом, в некотором интервале параметров состояния, в котором выполняется уравнение (9), избыточная вязкость Д^, в соответствии с (11), есть однозначная функция плотности. Этот экспериментальный факт хорошо известен и часто используется для описания вязкости газов малой и умеренной плотности.

0

v

а

и

Отличие состоит в том, что в качестве избыточной вязкости в литературе принимают

Дц = ц(Т, P) - ц0 (Т). Обычно зависимость избыточной вязкости от плотности описывают с помощью полиномов по степеням плотности.

Как видно из уравнений (11) и (12), коэффициент вязкости стремится к бесконечности

при плотностях р ^ —. Как известно, вязкость

стремиться к бесконечности при полном отсутствии трансляционного движения молекул, то есть при абсолютном нуле. Таким образом, параметр р в уравнениях (11) и (12) есть ни что иное, как удельный объем У0 при абсолютном нуле.

1

р 0 = — может 0 V

Численное значение

быть

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

Тогда уравнения (11) и (12) могут быть представлены в виде:

(., Л

ДЦуи = Ср2 ехр

(Т, р) = Леп + Ср2 ехР

б

22

1 - V 02р

1 - V02р2

(13)

(13а)

где р заменено на У0, а а заменено на б. Здесь, так же как и в уравнении (1а), в выражении для Цеп полагается, что Ь = V0.

Отметим также, что параметр а в уравнениях (11) и (12), вычисленный при описании опытных данных по вязкости может не совпадать с параметром а в "точных" уравнениях (1) и (1а), поскольку при выводе уравнений (11) и (12) использовано весьма грубое уравнение состояния Ван дер Ваальса.

Таким образом, в той области состояний, в которой избыточная вязкость может считаться однозначной функцией плотности, для описания опытных (или надежных табличных) данных в соответствии с уравнениями (7) и (9) необходимо иметь 2 подгоночных константы, а при отсутствии достоверных сведений о V0 - 3 подгоночных константы.

Нами проведена проверка применимости уравнений (13) и (13 а) для описания экспериментальных данных по вязкости различных веществ, и определена область параметров состояния, в которой избыточная вязкость по уравнению (13) в пределах погрешности экспериментальных данных остается однозначной функцией плотности.

На рисунке 4 приведена зависимость Дц = f (р) для диоксида углерода (исходные данные взяты из [8]). Как видно из рисунка 4, избыточная вязкость Дц действительно остается (в пределах некоторого интервала неопределенностей) однозначной функцией плотности вплоть до ртах ~ 0.8 г/см3, то есть примерно до 1,7ю. При более

высоких плотностях наблюдается существенное расслоение данных о Дц по изобарам, возрастающее по мере приближения к температуре плавления. Коэффициенты уравнения (13) для диоксида углерода, полученные методом наименьших квадратов, имеет следующие значения: С=19,704 мкПа-сек, б=1,087, У0 = 0,5848 см3/г; плотность р - в г/см3, значение У о для расчетов было взято из [8].

300-

200-

100

0

0,0 0,5 1,0 1,5 2,0 2,5

Рис. 4 - Избыточная вязкость диоксида углерода, как функция плотности

Проведенные расчеты показали, что уравнения (13), (13а) хорошо описывают экспериментальные данные по вязкости умерено плотного газа в достаточно широкой области параметров состояния.. Как показали наши расчеты, максимальные отклонения между рассчитанными по уравнению (13а) значениями вязкости и табличными значениями достигают 8% для диоксида углерода (в диапазоне температур от 220 до 1000 К и давлений от 1 до 1000 бар).

Если осуществить те или иные приближения в уравнениях (1), (1а) для соответствующих областей состояния вещества, можно получить и другие хорошо известные уравнения для расчета вязкости.

Рассмотрим жидкость вблизи температуры плавления в малом интервале изменения давления. В этом случае вязкость газа (при той же температуре) много меньше вязкости жидкости, то есть, как следует из уравнения (1а)

^ = ДЦ, (14)

И, поскольку зависимость вязкости жидкости от давления слабая и интервал давлений мал, можно считать, что вязкость жидкости зависит только от температуры. Как показали расчеты, в достаточно широком интервале температур выше тройной точки функция х = ДУ/v = х(Т, р) также слабо зависит от давления и в относительно широком интервале температур хорошо описывается зависимостью:

х = Х0-уТ, (15)

или

х = хг -у(Т - Тг), 15а где Хг - значение функции х в тройной точке, Тг -температура в тройной точке.

и

б

Тогда, с учетом (1), (14) и (15а) будем иметь:

цL = Aq = Ax exp

axn

V x0 - x у

= akr -у(( - Ttr )]exp

ax0

It

(16)

Пренебрегая температурной зависимостью в предэкспоненциальном множителе и объединяя постоянные коэффициенты, получим уравнение типа уравнения Аррениуса:

Е

ц = В • еЯТ, (17)

где В = А[хг - у(Т - Тг)], а энергия активации Н0Я

вязкого течения Е = —11

a .

v о У

Чтобы (17) выполнялось, необходимо, чтобы к = УТ—Тг) << 1. Для диоксида углерода

хг

х0 =1020 MPa, хг =405.2 MPa, и у = 2.839 К-1 в диапазоне температур от Ъг=216.58 К до Т=260 К и для давлений до 5 MPa. Считая достаточным условие к=0.1^0.2 получим, что интервал температур, в котором (16) переходит в (17), равен для диоксида углерода 15^30 К.

Как показали расчеты, ДУ в достаточно широкой окрестности тройной точки мало отличается от Н00 и слабо зависит от температуры и давления. В этом случае можно положить, что ДУг « Н0, и

'tr ~ ' '0

(

\

x 0 - x

AH 0/v 0

H 0

H 00/v 0 - AU tri v

= 1 + -

(18)

v0 H0 (v - v0) v - v0

v - v0

Подставляя полученные выражения в уравнение (1) с учетом уравнения (14), будем иметь:

-L = A—exp v

V v - v0 у

(19)

Пренебрегая слабым изменением плотности (объема) в достаточно малом диапазоне температур и давлений вблизи тройной точки, можно записать:

qL = C exp

B

V v - v 0 у

(20)

Н°

где С = А——; В = аv0, то есть получим известное V

уравнение Дулиттла. Один из вариантов этого уравнения, полученного в модели свободного объема, приведен в монографии [10]. Это уравнение по (10) имеет вид:

ЯТ 3 й

-L =

v - v 0 2 ^е

exp

f Nv' Л

V v - v 0 у

(21)

Здесь й = Ь/2л - постоянная Планка, N - число

Авогадро, к - постоянная Больцмана, 6 -температура Дебая, V'- объем частицы,

принимающей участие в вязком течении, V и v0 -мольный объем вещества при температуре T и абсолютном нуле. Формула (17) не имеет эмпирических констант и удовлетворительно выполняется, в частности, для жидких металлов: натрия, калия, олова и ртути [10] При этом v'- объем структурной единицы течения - принят равным объему иона соответствующего металла.

Показано, что избыточная вязкость жидкости, газа и флюида является однозначной функцией плотности энергии взаимодействия. Получена расчетная формула для описания коэффициента вязкости диоксида углерода. Для расчета плотности и внутренней энергии вещества, которые входят в расчетную формулу для вязкости, необходимо располагать термическим и калорическим уравнениями состояния либо табулированными значениями плотности и внутренней энергии в зависимости от температуры и давления (см, например, [3, 5, 7, 11-13]) В настоящей работе мы использовали предложенное нами малопараметрическое единое уравнение состояния жидкости и газа [7], содержащее 10 эмпирических коэффициентов, часть из которых имеет физический смысл. В уравнение для вязкости входят от двух до четырех индивидуальных эмпирических коэффициентов, которые

рассчитываются по экспериментальным (надежным табличным) данным по вязкости. Проведен расчет коэффициента вязкости в широком интервале параметров состояния для диоксида углерода. Отклонение рассчитанных значений вязкости от табличных в основном не превышают расхождений между экспериментальными данными разных авторов.

Показано, что ряд наиболее известных уравнений для расчета вязкости жидкости и газа в различных областях термодинамической поверхности являются частными случаями предложенных уравнений. Проведенный анализ позволяет лучше установить область применимости этих уравнений.

Литература

1. J. O. Hirschfelder, C. F. Curtiss, and R. B. Bird. Molecular Theory of Gases and Liquids. Wiley, New York, 1954. p.162-165.

2. H.N.V. Temperley, J.S. Rowlinson, and G.S. Rushbrooke, Eds. Physics of simple liquids. North-Holland, Amsterdam, 1968. p. 36-43.

3. А.Б. Каплун. ТВТ. 27, №5, 884-888. (1989).

4. A.B. Kaplun and A.B. Meshalkin, High Temp. High Press., 31, 253-257. (1999).

5. A.B. Kaplun, A.B. Meshalkin. J. Mol. Liquid., 120, 103106. (2005).

6. А.Б. Каплун, А.Б. Мешалкин. ДАН, 48, 490-493. (2003).

7. А.Б. Каплун, А.Б. Мешалкин. ЖФХ, 87, №8, 12941299. (2013).

8. В.В. Алтунин. Теплофизические свойства двуокиси углерода. Из-во стандартов. М., 1975.

9. A. Fenghour, W.A. Wakeham, V. Vesovic. J. Phys. Chem. Ref. Data. 17, N1, 31-65. (1998).

0

v0v

v

v

0

0

a

10. А.Н. Соловьёв, А.Б. Каплун. Вибрационный метод измерения вязкости жидкостей. Новосибирск, Наука. 1970. с. 135-140.

11. А.Б. Каплун, А.Б. Мешалкин. Вестник Казанского технологического университета, № 8, 88-91 (2010).

12. Дьяконов Г.С., Динмухаметова Р.А., Клинов А.В., Дьяконов С.Г. Вестник Казанского технологического университета, 16, № 4, С. 33-37. (2013).

13. Дьяконов Г.С., Казанцев С.А., Дьяконов С.Г. Вестник Казанского технологического университета, 16, № 4, С. 42-44. (2013).

© А. Б. Каплун, доктор технических наук, главный научный сотрудник лаборатории термодинамики веществ и материалов Институт теплофизики СО РАН, kaplun@itp.nsc.ru; А. Б. Мешалкин, доктор физико-математических наук, ведущий научный сотрудник лаборатории термодинамики веществ и материалов Институт теплофизики СО РАН, kaplun@itp.nsc.ru.

© A. B. Kaplun, Dr., chief scientist of substances and matter thermodynamics laboratory of Institute of Thermophysics SB RAS, kaplun@itp.nsc.ru; A. B. Meshalkin, Dr., leading scientist of substances and matter thermodynamics laboratory of Institute of Thermophysics SB RAS, kaplun@itp.nsc.ru.