Приближенная теория термического уравнения состояния Текст научной статьи по специальности «Математика»

УДК 541.12.011

Г. С. Дьяконов, Р. А. Динмухаметова, А. В. Клинов, С. Г. Дьяконов

ПРИБЛИЖЕННАЯ ТЕОРИЯ ТЕРМИЧЕСКОГО УРАВНЕНИЯ СОСТОЯНИЯ

Ключевые слова: термическое уравнение состояния, принципы замыкания аналитической термодинамики.

Предложено общее решение задачи о термическом уравнении состояния в терминах давление, энтропия с использованием принципа обратимости в виде соотношения Максвелла.

Keywords: thermal equation of state, the principles of closure of the analytical thermodynamics.

It is proposed a general solution for the problem of the thermal equation of state in terms ofpressure and entropy with the principle of reversibility such as Maxwell relationship.

Ранее [1] нами было предложено уравнение состояния, в котором давление выражалось через разность текущей и идеальной энтропии:

Zv(n*,T*) = 1 + n*B2 (T*)

+ an

f —AS (n* ,T *) A a — 1

(1)

где Zv =

P

pkj

- фактор сжимаемости; B2 - второй

вириальный коэффициент; p = N - числовая плотность;

N - число частиц; kB - константа Больцмана; S(n',T') — S0(n',T')

AS(n" ,T *) =

Nka

разница

между

реальной энтропией и идеальной энтропией; для Леннард-Джонсовых (ЛД) систем параметр а = 2,5 . Здесь и далее используются следующее безразмерные величины: кТ

T =

n = pa

a - эффективный диаметр

молекулы, е - глубина потенциальной ямы.

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

dZv

dT *

+ Zv — 1 = —n*

dAS

dn*

(2)

для ее замыкания.

Таким образом, имеем два уравнения (1), (2) с двумя неизвестными функциями - давление и энтропия, которые могут быть разрешены.

Подстановка (1) в (2) приводит к одному уравнению в частных производных первого порядка:

¡1 \ \ (п*,Т*)

-(IV(,Т )-1) + „ П '

Zv(n*,T*) — 1 + (a — B2 (T*))

dZv (n* ,T*)

= T*-i-'- + Zv(n*,T ).

лт* V /

(3)

дТ

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

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

Ч -=п, = 1 + П0В2 . (4)

Необходимо отметить, что используемое уравнение (1) содержит параметр а, как множитель подобия, значение которого было принято равным 2,5 в соответствии со скейлинговым законом [2].

Однако если следовать результатам, которые изложены в предыдущей статье [1], нетрудно увидеть, что постоянство параметра а , который является степенью при комплексе

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

Для удобства в выражении (3) введем новую искомую функцию, которая связана с соотношением

IV -1

u = T*

- +a — B,

n

После преобразований получим нелинейное уравнение в частных производных первого порядка

du

1 du

dT* pu dn* с граничным условием

u\ *

In

Решение уравнения (5) имеет вид:

d (B T *) dT*

* = aT *

Zv(n*,T*) = 1 + T_ B2(t)t,

T*

n*- n*0 = i 1

dx

(5)

(6)

(7)

1

t-^x) - a)x + —- B2(t)t

a a

Сравнение результатов расчета давления по (7) с данными численного эксперимента представлено на рис. 1 и в таблице. Расхождение F в таблице вычисляется по формуле

(п*,Т *)- (п*,Т *)

F=100-

(8)

ZvMD (п*,Т *)

где ZvMD (п* , Т *) - значение сжимаемости, полученное с

помощью вычислительного эксперимента методами молекулярной динамики и Монте-Карло, аппроксимированные Бенедиктом-Веббом-Рубином (БВР) [3], Zv (п*,Т *) - найденное по уравнению (7).

Рис. 1 - Фактор сжимаемости для Леннард-Джонсовых флюидов. Линии - расчет по (7), геометрические фигуры - данные численного эксперимента [3]

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

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

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

В частности при решении гидродинамических задач часто используется коэффициент гидравлического сопротивления

С,:

P

ро2/2

= Cf,

где P - давление, р - плотность, и - скорость движения среды.

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

Cf = ^е). (9)

Вид этой функции одинаков для подобных явлений. Поскольку число Рейнольдса это комплекс подобия, определяемый отношением сил инерции к силам вязкости и определяющейся через масштабную скорость, масштабный размер и кинематическую вязкость, нетрудно увидеть аналогию вида eAS = Re. Если при этом принять, что характерная скорость явления равна средней скорости молекулярного движения, то есть #7

v =-

m

характерная длина это

характерный размер объема, в котором происходит явление, то есть x = ^V, то эти значения определяют величину идеальной

энтропии = -ln(n*) + Cv ln(T*) + const.

В то же время текущая энтропия S определяет текущее значение характерной скорости и длины. Для подобных явлений Cf = Л Re*, где А и k константы. В том случае, когда зафиксировали

показатель

ДБ

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

3

2

0

0

0.2

0.4

0.6

0.8

n

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

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

dU

dV T

dU T dS

dV T dV T

тт ди

где U - внутренняя энергия системы, ^^

по объему

T

при

производная

постоянной

внутренней энергии температуре.

Вводя его (9) также в виде степенной функции, получаем окончательный вид общего решения.

Zv = 1 + n* B2 +

1 */ +—n (e

P У

,-AS P

" 1

du

dV dU_ dV

"dU / '( dU dS

dV rl [dV T dV T У

(10)

+ T

dS dV

= P,

= T dP

= dT

- P,

откуда

dU dP

dV T dT V

dU T dS - P

dV T dV T

- P

= 1 -

d In (P)

вид

d In (T)'

После преобразований уравнение (10) принимает

dIn (P,) "

Zv = 1 + n*B2 +

P n*( ASP- 1)

P

1 -

d In (T *)

(11)

где добавка в виде разницы берется по модулю.

Полученная система уравнений оказывается нелинейной и достаточно сложной. Чтобы избежать численного решения, используем приближение, в котором поправку, введенную выше, будем вычислять на основе исходного невозмущенного решения (7) и тогда получим учет первого приближения к найденному ранее решению. То есть Р0 и Д50 - вычисляется с помощью (7) и (2) соответственно.

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

термодинамики).

Для нахождения решения фазовое пространство разбивается на три области по температуре: 0,75 < T*< 1,0 , 1,0 < T*< 1,33 , 1,33 < T*< 5,0 .

Первая - область низких температур 0,75 < T*< 1,0 , где для нахождения параметров P и у уравнения (11) используется асимптотика S - функции при температуре T =0,8. Решение с S - функцией было рассмотрено в ранее выполненной работе (Моделирование замыкания аналитической термодинамики). Для этого параметры общего решения необходимо изменять до тех пор, пока искомое решение не будет иметь пересечение с асимптотическим решением Z при Т =0,8. Находим, что при плотностях равных 0,845 и 0,895 P = 0,4256, у = -0,0797. Средняя ошибка общего решения (10) относительно численного экспериментом [3] составляет 6,2%.

Следующая - область средних температур, которая является наиболее сложной, и для определения параметров используются обе асимптотики. Сначала оба параметра необходимо найти на критической изотерме, путем минимизации функции вида

sum = £| F (n* ,T *j )|, (12)

i, j

где F определяется по (8) и в качестве ZvMD выступает решение по модели квадратичной связи функционалов, тогда P и у оказываются равными 0,449 и -0,0922 соответственно. Фиксируется полученное значение у и на изотерме T =1,0 путем минимизации функции (12), где в качестве ZvMD выступает решение с S - функцией, находим значение P = 0,4381. Линейно аппроксимируя P в рассматриваемой области, получим функцию вида p(T*) = 0,0335 T* + 0,4046 . Затем на изотерме Т =1,34 вычисляется степень у минимизацией выражения (12), где интервал по плотности берется от 0,25 до 0,6. Искомый параметр оказывается равным у = -0,09435 . Тогда средняя ошибка общего решения (10) относительно численного эксперимента [3] в рассматриваемой области составляет 6,14%.

В области температур выше критической для определения параметров общего решения (11) используется асимптотика квадратичной связи

функционалов на температуре T =1,33 и на интервале по плотности 0.25 < n* < 0.6 . Причем значение P остается равным 0,4

1

(P = — = 0,4), а параметр у находится путем

а

минимизации выражения (12), откуда

вычисляется у = -0,062498. Максимальная ошибка в

рассматриваемой области составляет 17%, средняя -4,3%. Результаты сравнения решения, найденного по (10) с численным экспериментом [3] по всей области фазовой диаграммы можно увидеть в таблице.

Проверка общего решения (11) на основе того же массива опытных данных молекулярной динамики и Монте-Карло в аппроксимации БВР [3] показала, что среднее расхождение не превышает 6%.

Аналогичным образом по зонам найдено решение уравнения, рассмотренного в ранее выполненной работе (Принципы замыкания термодинамики Леннард-Джонсовских флюидов), в которой избыточная энтропия аппроксимируется степенной функцией вида:

АП-к

АБ = ААП—, (13)

Т-т ' 4 '

где А, к и т параметры, то есть уравнение (1) принимает вид

( -АБ (п* ,Т*) А

(п*,Т*) = 1 + п*В2 (Т*)-

) + аП

-1

(14)

Для нахождения параметров фазовое пространство разбивается на четыре области по плотности и температуре: первая - Ж <п* < 0,9, 0,7 < Т* < 1,2, вторая - NL < п < 0,6, 1,2 < Т < 1,31, третья -0,0 < п < 0,6, 1,31 < Т < 5,0 и четвертая - 0,6 < п < 1,0, 1,2 < Т < 5,0, где МЬ - плотность на линии кипения.

В первой области используются обе асимптотики для нахождения параметров: ниже температуры 1,2 -асимптотика с 8 -функцией, а при температуре 1,2-асимптотика с квадратичной связью функционалов. После минимизации функции (12) значения констант в уравнении (14) оказываются равными A=-4,26775, m=0,432026, q=1,68986.

Для нахождения параметров во второй области, как и в последующих двух, используется квадратичная асимптотика, и минимизируется выражение (12). Тогда для второй области константы уравнения (14) равны -A=-3,23293, m=0,593341, q=1,0147, для третьей - A=-3,443, гп=0,47452, q=1,13906 и для четвертой - A=-4,24731, m=0,215733, q=1,68029. Результаты вычислений представлены в таблице.

Выводы

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

Таблица 1 - Результаты расчетов по предложенным уравнениям состояния. Расхождение Р (%) вычислено по (8). Двойными линиями обведены области фазового пространства, где расчеты выполнялись с фиксированными параметрами, значения которых приводятся выше в тексте

* п Т г(3) Р (7) Р (11) Б (20)

0,83 0,75 0,187 -190,5 11,287 -0,90

0,85 0,75 0,674 -28,8 -6,285 -0,01

0,87 0,75 1,228 -0,12 -6,116 0,63

0,89 0,75 1,856 12,0 -5,007 1,12

0,9 0,75 2,198 15,7 -4,560 1,31

0,76 0,9 0,126 -232,3 47,135 -20,95

0,78 0,9 0,415 -50,2 9,970 -4,99

0,8 0,9 0,748 -14,0 5,630 -1,07

0,82 0,9 1,128 1,67 4,962 1,00

0,84 0,9 1,558 10,48 5,071 2,39

0,86 0,9 2,043 16,1 5,278 3,42

0,87 0,9 2,307 18,2 5,341 3,84

0,56 1,2 0,102 -163,0 8,913 -11,80

0,59 1,2 0,171 -107,9 13,04 2,99

0,68 1,2 0,628 -22,3 10,70 -8,31

0,76 1,2 1,450 5,7 1,89 -3,387

0,84 1,2 2,790 18,4 1,499 0,60

0,88 1,2 3,706 22,1 2,183 1,97

0,05 1,33 0,850 0,324 0,055 0,21

0,15 1,33 0,593 0,863 0,263 0,82

0,25 1,33 0,410 1,886 3,931 3,17

0,35 1,33 0,304 -1,238 11,319 5,46

0,45 1,33 0,269 -21,81 5,216 -1,29

0,55 1,33 0,371 -34,7 13,424 0,27

0,65 1,33 0,775 -13,4 11,128 -7,71

0,75 1,33 1,684 7,788 2,324 -4,67

0,85 1,33 3,310 19,45 5,397 -0,91

0,95 1,33 5,962 24,9 6,068 1,21

1 1,33 7,772 25,7 4,227 1,27

0,05 4 1,017 0,134 -0,03 0,36

0,1 4 1,043 0,269 -0,14 1,13

0,15 4 1,078 0,318 -0,39 1,97

0,2 4 1,123 0,274 -0,76 2,69

0,25 4 1,180 0,183 -1,49 3,18

0,3 4 1,251 0,125 -2,21 3,39

0,35 4 1,340 0,196 -2,57 3,33

0,4 4 1,451 0,483 -2,55 3,01

0,45 4 1,588 1,053 -2,22 2,51

0,5 4 1,756 1,928 -1,70 1,88

0,55 4 1,962 3,085 -1,19 1,15

0,6 4 2,212 4,454 -0,87 0,35

0,65 4 2,513 5,929 -0,95 -0,54

0,7 4 2,872 7,387 -1,59 -1,53

0,75 4 3,297 8,705 -2,98 -2,66

0,8 4 3,798 9,771 -5,40 -3,97

0,85 4 4,384 10,485 -9,94 -5,47

0,95 4 5,863 10,479 -11,34 -9,16

1 4 6,782 9,521 -8,67 -11,40

Работа выполнена при поддержке РФФИ грант №12-08-00465-а.

Литература 2. И. З. Фишер Статистическая теория

жидкостей, Наука, Москва, 1961. 280с; 1. А. В. Клинов, С.А. Казанцев, Г.С. Дьяконов, С.Г. Дьяконов, 3 J. Johnson , J. Zoll^^ К. Gubbins, Mol. Phys.:: 78,

Вестн. Казан. технол. ун-та, 13, 1, 22-27 (2010); 3 591-618 (1"3).

© Г. С. Дьяконов - д-р хим. наук, проф., член-корреспондент АН РТ, вице-президент АН РТ, ректор КНИТУ, [email protected]; Р. А. Динмухаметова - асп. каф. процессов и аппаратов химической технологии КНИТУ, [email protected]; А. В. Клинов - д-р техн. наук, проф., зав. каф. процессов и аппаратов химической технологии КНИТУ, [email protected]; С. Г. Дьяконов - д-р техн. наук, проф., акад. АН РТ, советник ректората КНИТУ.