Моделирование замыкания аналитической термодинамики Текст научной статьи по специальности «Математика»

УДК 541.12.011

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

А. В. Клинов, С. Г. Дьяконов

МОДЕЛИРОВАНИЕ ЗАМЫКАНИЯ АНАЛИТИЧЕСКОЙ ТЕРМОДИНАМИКИ

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

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

Keywords: thermodynamics, the equation of state, potential of the Lennard-Jones.

We propose five methods for closing the equations of thermodynamics and obtaining the equation of state. Calculations are performed for different points of the phase diagram in comparison with the results of numerical experiments for molecules with Lennard-Jones potential.

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

дЦ_

дV

= T

дР_

~дГ

- Р.

(1)

Кроме этого, при математическом описании процессов химической технологии используются соотношения статистической физики, которые

являются определением энергии системы и и ее давления Р:

U = - NkbT + 2nNp\p(r)g(r)r2dr .

P -=pKT-^p f p>....

2

Ф'

і dr

-g(r)r3dr.

(2)

(З)

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

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

p(r) = 4є

f Ґ \ 12 а

где ст — эффективный диаметр молекулы, е — глубина потенциальной ямы. Выбор такого потенциала

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

полученных из численных экспериментов [1] Монте-Карло (МК) и молекулярной динамики (МД), что позволяет проводить проверку результатов.

Если перейти к приведенным переменным согласно соотношениям Р' = Рс3/в, Т = кьТ / в,Р=рс3, и' = и/в и использовать термодинамические функции:

коэффициент сжимаемости

Z = ■

P

pkbT

внутреннюю энергию, приходящуюся на одну молекулу - Е = N, безразмерные функционалы сил отталкивания и сил притяжения

Фі (p ,T') = j[xj g(x,p',T')x2dx, Ф2 (p,T') = f f-1 g(x,p , T* )x2dx,

(4)

то выражения (1), (2) записываются в виде

E' = ф T' + 8p (Фі -Ф2 ),

(5)

Z = 1 + 1bnT- (2Ф1 -Ф2),

а выражение (З) можно представить следующим образом

' дZ 1 . fдE'

= -p •V \-r- /T"

^ I = -p I . (б)

Т чдТ" )р ' удр

Если в уравнение (6) подставить выражения (5), то получится квазилинейное неоднородное дифференциальное уравнение первого порядка в частных производных для используемых функционалов Ф-|, Ф2:

4• T

= 2 • T

дФ1

дг'

+p

дФ 2

дТ

+ p

дФ1 дР JT

дФ„

- 3<Ф =

(7)

КдР УТ

Впервые понятие краевой задачи в термодинамике было использовано в работе [2].

2б

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

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

Первая модель замыкания [3] базируется на конкретном виде функционалов для потенциала Леннарда-Джонса (4) и позволяет предположить, что приближенно замыкающим уравнением соотношение

2

Ф1 =Ф 2

будет

(В)

В этом случае уравнение (7) принимает вид [2•T• (4Ф2 - Df®2

-[p (2 -DHt

дT

= ФФ22 -Ф2,

(9)

которое идентично дифференциалах [4] dT*

системе уравнений в полных

dp

.Ф„

ФФ22 - Ф2

2 • T* (4Ф2 -1) p (2Ф2 -1)

данных уравнений являются

а интегралами функции

Т = —Ф

1 т*0.25

p

*4 (

p

1

Чф 2

- Ф

(10)

две

(11)

Для получения конкретного решения дифференциального уравнения первого порядка в частных производных от двух переменных требуется граничное условие первого типа, т. е. знание искомой функции на какой-либо линии. Согласно работе [З] для уравнения состояния простого вещества имеется характерная линия в виде прямой, называемая линией Zeno, для точек которой Z( T, p) = 1 , и согласно (З)

ф2 (ІЗ)

Ф1 =■

2

В частности, для вещества типа Леннард-Джонса уравнение этой прямой и ее параметры будут следующими

Т-т + Р = 1, Г; = 1,1, Рв' = 3,418 . (13)

Тв рв

Используя соотношение (8) и граничное условие (12) для линии Zeno, находим, что на этой линии Фх1 = 0,25, Фх2 = 0,5 и, согласно (8), (9), можно получить аналитическое неявное выражение для расчета искомой функции Ф2 в следующем виде

Q- 3Ф2 ) Р

0,25

0,25 T

л

= 1. (14)

^(1 - зф2 )2 -2 Т ‘ ‘

После нахождения Ф2, по формуле (12) можно найти Ф1, а по формулам (5) 7 и Е. Вместо значений

Фх1=0,25, Фх2 = 0,5 для линии Zeno можно взять экспериментальные данные MD, которые близки к этим значениям, при этом результаты расчетов совпадают в пределах погрешностей данных MD.

Вторая модель замыкания [6] основана на предположении о том, что уравнения (5), (6) имеют фундаментальный характер и их структура не зависит от вида потенциала взаимодействия. Это означает, что они справедливы и для крайних случаев: 1 - имеются только силы притяжения (например, система гравитационно взаимодействующих тел в космосе), 2 - имеются только силы отталкивания [3] (например, система одноименно заряженных частиц в плазме). При таком предположении для функционалов Ф1 и Ф2 должны быть справедливы независимые уравнения,

позволяет замкнуть исходные соотношения

что

дФ1

2 T

дT

дФ 2

~дГ

дФ

4 •T I —г I +p I—11 - ЗФ1 = 0, (ІЗ)

+ p

дp

дФ 2 дp*

- Ф 2 = 0 .

(1б)

Для упрощения решения этих уравнений сделаем переход к другим функционалам Р-|, Р2.

16^p F1 = ^“ ф1,

16^p

F2 =^~ Ф 2

(І7)

T 1 2 T

тогда уравнения (19), @0) примут более простой вид

Р_

4

2

дP

дF2

дp*

+ T

3Fl

дT *

дF2

дТ

= 0,

= 0.

(18)

Для этих уравнений решение получается в виде характеристических функций

F1 = F

1

V T*/

F2 = F2

fpp21 v T* J

(19)

которые имеют постоянное значение вдоль своих характеристических линий, проходящих через произвольную точку р , Т , для которой ищется решение. Уравнения этих линий

pp p2

!—*г = Const, !—*г = Const. T T

(З0)

В качестве граничного условия используем также линию Zeno со значениями на этой линии: Фх1 = 0,25, Фх2 = 0,5 или данные эксперимента. При использовании данной

методики сначала рассчитываются p.,, T, , p2, T2 - точки пересечения характеристических линий с линией Zeno путем решения четырех уравнений

(ЗІ)

*4 *4 ^* *

p________= p1 T1 + p1 = 1

p2

pB

p2

(ЗІ)

р ________

T * T2* 'TB

Далее, зная значения Фх1(1), Фх2(2) в точках пересечения с линией Zeno, находим согласно

З7

(17) значения характеристических функций Рх1, Рх2 в этих точках, которые совпадают со значениями

характеристических функций Р1, Р2 в точке р , Т .

После этого по (17) находятся значения функционалов

Ф1, Ф2 в точке р , Т , что приводит к уравнениям

-1 = Фх 1(1), Ф2 = ІРф х2 (2), (22) Т1 Р Т2 Р

и далее по формулам (5) можно найти 7 и Е*.

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

д (х ) = £( х - X (р,Т)) (23)

В этом случае функционалы сил будут выражаться через одну функцию Х(р, Т),

( „ Л10 ( „ Л4

Ф1 (р,т ) =

1

, Ф2 (pT) =

1

, (З4)

что позволяет замкнуть исходные соотношения и получить одно уравнение для Х(р, Т)

[в-Т*(6 -5)}ГдХ

it

[ 2p(2X6 -5 )}[fp I =X (Ф -х6),

IX

|p

которое идентично системе уравнений

dI dp

(ЗЗ)

dX

8 • —* (X6 - 5) 2p (2X6 - 5) X (Ф - X6)'

(Зб)

а интегралами данных уравнений являются две функции

—*9/8 X15

;г-г, Т2 =p*9XФ0 (X6 - Ф) 07)

X6 - Ф

На линии Zeno искомая функция X = 2 уравнения

9 15 -T-*8 v^15

t *8 26 = - X

1/6

что дает два

Ф - X6

f

1 —

pB

(ЗВ)

|

B|

t5 = p9 XФ0 (Ф - X6),

совместное решение которых позволяет найти функцию Х(р,Т) и далее Ф-|, Ф2 , 7 согласно (24), (5).

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

пропорциональности С.

pc (r) = -C 4є\а Если ввести понятие функционала Сьюзерленда

(З9)

“ f 11(

Фс (p,T*) = J|-' g(x,p ,T )x2dx, (З0)

то уравнение (З) примет вид

dB d-

.p C(

Ф dB

E* = ФT* - T*p ^ - 8*p^c, 2 d— c

(З1)

Z = 1 + pB2 - 8^-*- Фc, и 2 T Ф c

После их подстановки в (б) получим

1' • ФФ

или в полных дифференциалах dT = dp = dФс

p Ф Iі-1 ф c V ф

(ЗЗ)

интегралами данных уравнений являются функции

з з

ф у-Ф —у

C , Т2 = p.

p

Ф

T *(

(З4)

(ЗЗ)

На линии Zeno функционал описывается уравнением

T'B2 (Г) 3 8nCy

и из (34) получим два уравнения

f tB2 (t )з V-3 8яСу

Фс (p,—*)'

3_

у-3

ф

T *у

Ф tу

1 -л 1

—в I

pB

ф

Iу

p

(Зб)

совместное решение которых позволяет найти функцию Фс(р,Т) и далее 7 согласно (31), при

этом, как оказалось, 7 не зависит от

коэффициента пропорциональности С.

Наилучшее совпадение с данными МЭ

получилось при у = 28.

Пятая модель замыкания базируется на возможности представления решения уравнения (7) в следующем виде

Ф1 (р*, Т•) = ф10 (р*, Т•) + Р Ф0 (р*, Т* ) ,

Ф 2 (р, Т *) = Ф 20 (р, Т *) + д Фо (р, Т *),

где Ф10 и Ф20 являются решениями уравнений (15), (16) во второй модели и выражаются соотношениями (22). Функция Ф0 (р*,Т*)

является добавочной поправкой, для которой с помощью (7) можно получить

дифференциальное уравнение

(З7)

2 •—* (2p - q)

|ф0

И

-p(p - q)

іф0

dp*

■. (ЗВ)

- (3p-q)Ф0 = 0

16np'

Если сделать замену F0 =-*—Ф 0, то получим

9

ЗВ

+ 2T • (2p - q Ш

/t V /p

и решение в виде характеристической функции

С ,*2p(2p-q)Л

Fo = Fo

p

r*(p-q)

= 0 (39)

(40)

которая имеет постоянное значение вдоль

характеристической линии, уравнение линии имеет вид

р*2(2 p-q)

• = Const. (41)

p

T "(p-q)

В качестве граничного условия используем также линию Zeno со значениями на этой линии: Фх-| = 0,25, Фх2 = 0,5 или данные эксперимента. При

использовании данной методики необходимо, как и во второй методике, рассчитать p1,T1, р2,Т2 - точки

пересечения характеристических линий с линией Zeno путем решения четырех уравнений (21), а также точку

пересечения характеристической линии p0, Т0 согласно уравнениям

p

•2(2 p-q)

Po

*2(2 p-q)

r*(p-q)

r*(p-q)

T0 + p0 = 1 ’t; * •

(41)

(42)

'0 ■ в Pb

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

*2 p(2 p-q) т *

Fo =

p

Ф o = -

(43)

T*(p-q) ’ _0 16яр В этом случае из (37) получим для линии Zeno

Ф10 (р*,т*) = Ф Х1 + p Ф0 (p*,t*) ,

Ф20 (р*,т*) = ф х 2 + q ф0 (p*,t*)

и, используя свойства характеристических функций F1 F2, F0, получим выражения для функционалов Ф10, Ф20

Ф0 в произвольной точке р , T

(44)

Ф1

(

p

Ф20 (P*,T* ) = TP

V ’ P

p P1

p Ф x 1 -T x 1 16^

q P2

pp ф x 2 -t2 x2 16^

■*(p-q)

I

2(2p-q) Л

t "(p-q) 2

(45)

Ф0 (p*,t * )=:

p0

2 p(2 p-q)

■*(p-q)

16яр

Далее, согласно (37) и (5), можно найти 7 и Е . В данной методике имеется два подгоночных параметра р и д , если минимизировать погрешности расчетов по этой методике и МБ для 7 на линии равновесия двух фаз то получим р=3,909206, д=2,271164.

Результаты расчетов по данным методикам даны в таблице и на рис.1. В табл. 1 для разных точек

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

фактора сжимаемости 7 и результаты расчетов: по изложенным выше пяти методикам.

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

Рис. 1 - Значения фактора сжимаемости 7 для границы двухфазной области.

Экспериментальные значения - 7уМй [1], результаты расчетов: по первой методике - 1, по второй - 2, по третьей - 3, по четвертой - 4 и по пятой методике - 5

Таблица І

Результаты расчетов по

P* Т * Z[1] s\ % S2, % S3, % S4, % f/? vO S%

1 2 3 4 5 б 7 8

0,04 1,08 0,82 3,21 6,52 -2,65 -0,25 5,55

0,04 1,45 0,9 0,8 2,33 -2,4 -0,15 1,73

0,04 1,81 0,94 0,13 0,91 -1,75 -0,08 0,52

0,04 2,18 0,96 -0,06 0,33 -1,16 -0,04 0,09

0,04 2,54 0,98 -0,09 0,09 -0,69 -0,02 -0,05

0,04 2,91 0,99 -0,06 0,01 -0,31 -0,01 -0,06

0,04 3,27 1,0 0,0 0,0 0,01 0,0 0,0

0,04 3,64 1,01 0,07 0,04 0,27 0,0 0,09

0,04 4,0 1,01 0,14 0,09 0,49 0,01 0,19

0,2 1,45 0,57 1,3 2,17 -18,8 -4,22 1,86

0,2 1,81 0,76 -0,51 -0,46 -8,93 -1,14 -0,77

0,2 2,18 0,88 -0,56 -0,66 -4,19 -0,3 -0,86

0,2 2,54 0,96 -0,25 -0,31 -1,33 -0,05 -0,39

0,2 2,91 1,02 0,14 0,18 0,6 0,01 0,22

0,2 3,27 1,06 0,52 0,7 2,0 0,0 0,85

0,2 3,64 1,1 0,87 1,2 3,07 -0,03 1,45

Окончание табл. 1

1 2 3 4 5 б 7 8

0,2 4,0 1,12 1,2 1,68 3,91 -0,08 2,02

0,36 1,45 0,43 -5,89 -21,2 -39,0 -13,4 -14,1

0,36 1,81 0,74 -1,72 -5,51 -10,2 -1,72 -3,99

0,36 2,18 0,94 -0,36 -1,01 -1,81 -0,13 -0,79

0,36 2,54 1,08 0,5 1,25 2,17 0,02 1,04

0,36 2,91 1,18 1,17 2,72 4,55 -0,13 2,35

0,36 3,27 1,26 1,73 3,81 6,14 -0,33 3,41

0,36 3,64 1,32 2,21 4,68 7,3 -0,53 4,29

0,36 4,0 1,37 2,63 5,4 8,18 -0,69 5,06

0,53 1,45 0,53 -0,96 -18,1 -17,3 -3,33 -7,69

0,53 1,81 0,99 -0,02 -0,17 -0,18 0,01 -0,08

0,53 2,18 1,28 0,63 4,15 4,37 -0,6 2,16

0,53 2,54 1,48 1,3 6,19 6,67 -1,24 3,58

0,53 2,91 1,62 1,97 7,48 8,16 -1,69 4,7

0,53 3,27 1,72 2,58 8,43 9,25 -1,99 5,66

0,53 3,64 1,8 3,14 9,19 10,09 -2,18 6,5

0,53 4,0 1,86 3,63 9,82 10,76 -2,31 7,25

0,69 1,08 0,29 3,73 -43,6 -22,5 -14,5 -3,74

0,69 1,45 1,32 -0,21 3,93 2,19 -0,28 0,3

0,69 1,81 1,88 0,06 7,55 4,57 -2,01 0,89

0,69 2,18 2,22 0,72 9,02 5,9 -3,14 1,67

0,69 2,54 2,43 1,47 9,99 6,97 -3,71 2,58

0,69 2,91 2,57 2,22 10,75 7,9 -3,95 3,52

0,69 3,27 2,67 2,93 11,41 8,72 -4,0 4,44

0,69 3,64 2,73 3,58 11,98 9,43 -3,95 5,3

0,69 4,0 2,78 4,16 12,48 10,05 -3,85 6,09

0,85 0,72 0,39 -7,88 -17,4 -5,46 -50,8 7,91

0,85 1,08 2,64 3,6 9,81 3,32 8,91 -2,6

0,85 1,45 3,53 5,16 10,6 5,19 5,82 -1,82

0,85 1,81 3,97 5,63 10,67 5,89 3,27 -1,37

0,85 2,18 4,2 5,89 10,93 6,36 1,64 -0,87

0,85 2,54 4,32 6,17 11,33 6,83 0,69 -0,22

0,85 2,91 4,37 6,49 11,78 7,32 0,17 0,51

0,85 3,27 4,4 6,82 12,25 7,81 -0,08 1,28

0,85 3,64 4,4 7,15 12,71 8,3 -0,16 2,05

0,85 4,0 4,38 7,46 13,13 8,76 -0,14 2,81

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

Литература

1. J. Johnson, J. Zollweg, K. Gubbins, Mol. Phys., 78, 3, 591-б18 (1993);

2. М. Ф. Сарры, Журнал технической физики, 68, 10, 1-9(1998)

3. А.В. Клинов, С.А.Казанцев., Г.С. Дьяконов Г.С. Дьяконов С.Г., Вест. Казан. технол. ун-та, 12, 1, 10-1б (2010);

4. Э. Камке, Справочник по дифференциальным уравнениям в частных производных первого порядка. Наука, Москва,19бб. 200c;

5. E.M. Holleran, J.Phys.Chem.Physics, 47, 12, 5318 (19б7);

6. С.А. Казанцев, А.В. Клинов, Г.С. Дьяконов, С.Г Дьяконов, Вест. Казан. технол. ун-та, 12, 1, 17-21 (2010).

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