CC BY
36
УДК 532.593+536.715:546.3
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ТЕПЛОВОЙ СОСТАВЛЯЮЩЕЙ УРАВНЕНИЯ СОСТОЯНИЯ МОЛЕКУЛЯРНЫХ КРИСТАЛЛОВ
Ю.М. Ковалев
Данная работа посвящена построению математической модели уравнения состояния молекулярных кристаллов. Ее практическая ценность заключается в том, что все твердые взрывчатые вещества являются молекулярными кристаллами. Следовательно, разработав математическую модель уравнения состояния молекулярного кристалла, можно будет прогнозировать поведение твердых взрывчатых веществ при высоких давлениях и температурах. Сложность построения уравнения состояния молекулярного кристалла состоит в том, что большое число степеней свободы молекул, входящих в состав кристалла, не позволяет проведение прямых вычислений. В данной работе был предложен подход, который позволил использовать все лучшее, что есть в моделях Дебая и Эйнштейна для описания термодинамики кристаллов. Разделение частот нормальных колебаний в кристалле на высокочастотные и низкочастотные ( деформационные) колебания позволило получить аналитическое выражение для коэффициента Грюнайзена и параметры для тепловой составляющей уравнения состояния молекулярного кристалла.
Ключевые слова: уравнение состояния, давление, температура, энергия, коэффициент Грюнайзена.
Законы сохранения массы импульса и энергии лежат в основе математических моделей механики сплошных сред, термодинамики, электродинамики и т.д. Однако законы сохранения тте являются замкнутой системой. Требуются зависимости между входящими в уравнения сохранения величинами - уравнения состояния.
В настоящее время принято считать, что в уравнения состояния молекулярных кристаллов входит две составляющие: тепловая и холодная. Тепловая составляющая определяется колебательным движением молекул, входящих в состав кристалла, а холодная составляющая - изменением энергии взаимодействия внутри молекулы и между молекулами. Определение зависимости коэффициента Грюнайзена 7 от удельного объема является одной из основных задач при построении уравнений состояния твердых тел. Это обусловлено тем, что данная зависимость является связующей между тепловой и холодной составляющими уравнения состояния в соответствии с формулами Ландау - Слейтера, Дугдала - Мак-Дональда и т.д. [1]. Решению этой задачи посвящено достаточно большое количество как экспериментальных, так и теоретических работ [1-3]. Однако, теоретическое определение зависимостей, характеризующих поведение твердых взрывчатых веществ (ВВ), осложняется тем, что они относятся к молекулярным кристаллам, а молекулы, входящие в состав кристалла, обладают большим числом внутренних степеней свободы. Тем не меттее, и в этом случае можно получить аналитическую зависимость для коэффициента Гртопайзепа от удельного объема твердого тела.
Термодинамические свойства вещества полностью определяются, если известен один из термодинамических потенциалов. Удобно исходить из определения свободной энергии Гельмгольца ^ (У,Т), которая наиболее простым образом связана с моделью строения вещества:
^ = и + Ео + кТ^ 1п(1 - ехр(-ка)), Ео = 2 ^ Нта. (1)
Здесь и - энергия взаимодействия между атомами; V - удельный объем; Т - температура тела; к - постоянная Планка; ,ша - частоты нормальных колебаний; Ео - энергия пулевых колебаний.
Если К(V, Т) задана, то дифференцированием могут быть найдены выражения для всех измеряемых термодинамических величии:
дР ^ - _ (дК \ . (дР \ - (дР \
Б Vдт) у ’ Р \сЖ)т’ \сЖ)т \6Т)у
Известно, что энергия взаимодействия между атомами молекулярного кристалла складывается из внутримолекулярной и межмолекулярпой. Внутримолекулярная энергия состоит из энергии валентных взаимодействий и энергии певалептпых взаимодействий атомов внутри молекулы. Межмолекулярпая энергия является энергией не валентных взаимодействий атомов. В силу того, что внутримолекулярная энергия валентных взаимодействий существенно болытте энергии певалептпых взаимодействий между атомами как внутри молекулы, так и разных молекул, имеет смысл разделить энергию взаимодействия молекулярного кристалла па молекулярную им (энергия валентных взаимодействий) и кристаллическую и к (энергия невалентных взаимодействий). Если энергия и к зависит от пространственного
им
вштбитных связей и Всшентных углов.
Учитывая тот факт, что частоты нормальных колебаний внутри молекулы па порядок больше частот нормальных колебаний молекулы как целого и деформационных частот, можно ввести две характеристические температуры и разделить колебательную составляющую свободной энергии па низкочастотную и высоко частотную. В силу того, что частоты нормальных колебаний молекулы как целого и деформационные частоты определяются из-ик,
малытых колебаний и будут зависеть от объема.
Предположив возможность использования для низкочастотной составляющей свободной энергии подхода Дебая, а для высокочастотной - подхода Эйнштейна, перепишем выражение (1) в виде
вю
К = ик + им + Ео + ЗМЕТ^£21п(1 — ехр(—£)) ^£+ (2)
+ (ЗЫ — М) КТ 1п (1 — ехр (— ,
где К - универсальная газовая постоянная; М - число низкочастотных колебаний; N - число атомов в молекуле; ЗЫ — М - число высокочастотных колебаний; вв - характеристическая температура Дебая; вЕ - характеристическая температура Эйнштейна.
Если для одпоатомпого вещества с физической точки зрения модель Эйнштейна представляется малореалистичпой, то для молекулярных кристаллов, каждая молекула которых имеет свой набор частот, часть спектра, соответствующая оптическим частотам, может быть приближенно описана эйнштейновской моделью [4].
Интегрируя по частям выражение для низкочастотной составляющей свободной энергии К (^Т) (2) и вводя функцию Дебая О (х) по формуле приведенной в монографии [4, 5],
X
О(х)=Хз /
х3 У ехр (£) — 1:
о
получаем
F = UM + UK + MRT ^ln (1 — exp (—xD)) — D ^ + (3N — M) RTln (1 — exp (—xE)), (3)
где Xd = Xe = T'
Используя выражение (3), легко получить выражения для давления P и энтропии S:
P = — (—^ ^ — Ё£к — E — MRTD (xD) d (ln °D) I —
V dVJ T dV dV dV ( D) d (ln V) V
— (3N — M) RTxed1 (/l,nвв}—г---------------------------------------^-r; (4)
v ; E d (ln V) V (exp (xe) — 1)
1п (1 — exp (—xD)) — D (XD)
З
I
I MRT
exp ( xd ) D (xd )
1 — exp (—xD) З
дXD , c)T
I(3N - M) R 1п (1 - exp (-xe)) I (3N - M) RT16xp ( XE\ ^XM. ]
1 — exp (xE) д! )
1п (1 - exp (-XD)) - D (^D)
—MRD (xD) I (3N — M) R 1п (1 — exp (—xE)) —
(3N — M) R(rXlE .} . (5)
exp (Xe) — 1 J (5)
D<x> = expdr—T -x D'<x>-
x
Зная равенства (3) и (5), легко определить выражения для полной энергии E и теплоемкости при постоянном объеме Су.
E = F + TS = Um + Uk + E0 + MRTD (xd) + (3N — M) RTXe
exp (xE — 1)
Су = MR I 4D (xd) - —71 + (3N - M) Rx2E f exp(xE^ ^.
exp (xd ) - 1) (exp (xE) - 1)
Следуя определению коэффициента Гртопайзепа
YD (V ) = - d <‘" в0 > ,
D (V ) d (1п V) ’
выражение (4) можно записать в виде
= _ д£м_ дЦк_ E , MRTyd (V) D (XD) (6)
c)V c)V dV + V ' ()
Последний член в выражении (4) равен нулю, так как при разделении частот было сделано предположение о независимости высоких частот от объема.
Исходя из определения энергии пулевых колебаний и учитывая разделение частот, получаем выражения для функций Eo и
тв
1^, ЗЫ' — М' ЗМ' [ 2, , ЗЫ — М п. З„,„л
Ео = ~2_^ Нша =-----2-----Ьыа + ш кшйш =------- КвЕ + -МКвв (V)
а о
2 ^ 2 2 7 2 8
о
йЕо З МК^в (V) вв (V)
(7)
(IV 8 V
где N - число атомов, М1 - число низкочастотных колебаний в объеме V.
Подставляя выражение для производной от энергии нулевых колебаний по объему (7) в равенство (6), получаем уравнение для определения давления в виде
„ МКТ^в (V) /З Л „ „ йик
Р = —V—\8 х° + О )+ Ру у = —' (8)
Для получения возможной зависимости коэффициента Гргопайзепа для молекулярных кристаллов от объема ВВ может быть использован следующий подход.
Изотермический модуль сжатия вт (изотермическая сжимаемость) определяется выражением
± = Ст = _ у(дР\
вт V т ’ (9)
где Ст обозначает изотермическую скорость звука.
Подставляя выражение (8) в правую часть равенства (9), получаем
— = —V— вт дV
= - V
т
МКТ^в (V) ( З
V2
- Хв + О (хв Л +
МКТ^в (V ) /З . \ МКТ^в (V ) /З г.. Лдхв дРу
+ -----^ -хв + О (хв ) + -----------О (хв И ж + дР
т
= - V
МКТ^в (V) (З . .
-----V2---- V 8 Хв + О (хв ) ) +
МКТ^в (V) (З . \ МКТ^в (V) (З - , . Л дРу
+ V ( 8 Хв + О (хв ) )--------------V2 ( 8 Хв — (хв) + О (хв П +
МКТ
V
Ьв ^) + 1в(V)) (^8хв + ° (хв^ — y
дРу МКТ^в (V) Су в (хв)
дУ
V
—МКТ^'в (V) (8хв + о (хв))
где ^'в (V) _ производная по V от коэффициента Грюнайзена:
Су в (хв) = 4О (хв) —
3хв
дхв хв
= —-тт- 1в (V) .
ехр (хв) — 1 ’ дV V
Определим значение изотермической сжимаемости при Т ^ 0. С этой целью переходим в последнем выражении к пределу:
З З
Ит МКрТ (72 (V) + 7в ^^ - хв = д вв МКР Ьв (У) + 1в (У^ ; т ——о 8 8
— МКрТ (7в (V) + ^в (V)) о (хв) = 0;
]\т^МКрТ^в (V) Сув (хв) = 0;
ЗЗ
11ш МКТ^в (V) - хв = - МКвв 1'в (V) т— 0 8 8
Ито МКТ^в (V) О (хв) = 0.
Следовательно, при Т ^ 0 изотермическая сжимаемость определяется выражением вида
' 1 \ з з дР
Т) = дввМКр№ (V) + 1в (V)) — 8МКвв 1в (V) — VдУ•
вт
звука Ст формулой (9), т.е.
С2 3 3 дР
V = 8ввМКр № (V) + 1в (V)) — 8МКвв!'в (V) — VдpУy.
Полагая далее, что скорость звука определяется только упругими свойствами кристаллов, получаем следующее дифференциальное уравнение типа уравнения Бернулли для определения зависимости коэффициента Гргопайзепа от плотности:
1в ^) — V1в ^) — V 1в ^) = 0-
Данное уравнение заменой г = 7д(у) приводится к линейному дифференциальному уравнению
, г 1
г +V = —V ■
которое легко интегрируется. Зависимость коэффициента Гргопайзепа от плотности описывается выражением следующего вида:
V
7в
С -V'
где константа С определяется из условия 7в (^) = 7°- Похожее выражение для коэффициента Гргопайзепа было получено из других соображений в работах А.М. Молодца [3, 6]. В результате получим
О 1
7в = м у)----------гт-----• (10)
Vуo/ 1+ т“ (1 — у)
Формула (10) в случае слабого сжатия переходит в известное эмпирическое выражение
о ( V
7в=7Ч ^
которое пшроко используется при обработке экспериментальных.
Т=0
определяется только упругими свойствами кристаллов, позволило получить аналитическую зависимость коэффициента Гргопайзепа от объема.
Для того чтобы воспользоваться выражением для коэффициента Грюнайзена 7в (Ю); необходимо определить его значение 7°, соответствующее начальному удельному объему
Уо- С этой целью получим аналог коэффициента Грюнайзена для молекулярных кристаллов. Для этого раскроем выражение отношения коэффициента теплового расширения а к изотермическому модулю объемного сжатия вт [4]
а
вт
д (д£\ дУ [от)у\
(11)
т
т
У, получим
вт -
+МКТ
Ь ( 1 - ехр ( - Т
- 3 °(¥.
д (во дТ\Т
+
= 1В (Уо) Р°МЕ
40(хо) -
3хо
(ехр(хо) - 1)
(12)
Вводя обозначение Су о для выражения, стоящего в квадратных скобках в правой части равенства (12), получим аналог уравнения Грюнайзена в традиционной форме:
а _ (У)Суо
вт _ У
(13)
вт
Ст
_ С!
вт _ У,
равенство (13) запишется в виде
аСТ _ (У) Суо ■
(14)
Следовательно, если из эксперимента известны значения коэффициента объемного расширения а и изотермической скорости звука Ст и определен способ нахождения Суо, из уравнения (14) можно найти значение коэффициента Грюнайзена при начальном значении удельного объема. Это дает возможность использования выражения (10) для определения зависимостей коэффициента Грюнайзена молекулярных кристаллов от удельного объема.
При определении теплоемкости при постоянном объеме, относящейся к деформационной часта, проведем анализ частот внутримолекулярных колебаний атомов па примере кристалла тэна (СбНв^О^)- Силовые постоянные для расчета спектров нормальных колебаний были определены с помощью квантово-химического метода РМ-3, подробно описанного в монографии [7]. Этот метод был выбран по тому, дает хорошее совпадение с экспериментальными данными, приведенными в работе [8]. Автор выражает благодарность профессору А.В. Велику за предоставленные ПК - спектры для тэтта, приведение в таблице. Анализ конформации молекулы тэна показывает, что в качестве колебаний, не зависящих от кристаллического окружения, с уверенностью можно выбрать 28 колебаний валентных связей, 4 колебания валентных углов ттитрогрупп и 3 колебания валентных углов центрального атома углерода.
Следовательно, начальное значение величины (3Ж - М) для молекулы тэна равно 35. Зная частоты нормальных колебаний и число колебаний, относящихся к не деформационной части (3Ж - М), легко определить величину Сум по формуле
1
Таблица
Частоты нормальных колебаний молекулы тэна
№ 1 СМ 1 № 1 і^,СМ 1 № 1 і^,СМ 1 № 1 і^,СМ 1
1 0.0 23 228,93 45 635,19 67 1343,03
2 0.0 24 240,66 46 661,46 68 1352,47
3 0.0 25 257,97 47 737,19 69 1355,29
4 0.0 26 263,26 48 823,13 70 1370,46
5 0.0 27 284,28 49 891,99 71 1383,53
6 0.0 28 293,78 50 917,14 72 1519,57
7 19,546 29 349,13 51 933,68 73 1526,12
8 33,861 30 375,19 52 944,77 74 1526,60
9 35,332 31 225,47 53 1097,99 75 1537,44
10 40,181 32 462,83 54 1110,75 76 2125,80
И 44,269 33 475,73 55 1130,86 77 2137,02
12 60,729 34 481,61 56 1146,12 78 2142,39
13 63,889 35 498,15 57 1150,42 79 2164,06
14 65,898 36 515,86 58 1163,09 80 2917,08
15 69,777 37 517,87 59 1178,33 81 2918,24
16 77,231 38 530,22 60 1183,26 82 2922,99
17 89,939 39 585,72 61 1194,98 83 2923,84
18 100,819 40 600,38 62 1250,72 84 2964,46
19 123,008 41 601,65 63 1282,43 85 2990,96
20 174,109 42 616,82 64 1310,43 86 2992,74
21 181,914 43 628,04 65 1314,59 87 2995,86
22 185,601 44 634,25 66 1330,37
3М о /
Сум, = Н £ ( х , (15)
і=М+1 (ЄХР х) - 1)
где Хі = кі, а значения Ші для тэна выбираются из таблицы соответственно. Зная величины Сум и Су, легко определить значение теплоемкости Су о, относящуюся к деформационным колебаниям
Суо = Су — Сум = МЕ^4В (хо) - 6Хр(Х£В — 1) ’ (16)
Легко показать, что выражение, стоящее в скобках правой части формулы (16), есть результат интегрирования по частям выражения Вс (х)
X
Вс(х) = Х3/Є
_ ехр({)
х3 } 4 (ехр (£) — 1)
йх,
которое представляет собой так называемую функцию теплоемкости Дебая. Данная функция при изменении хв от 0,05 до 1,0 изменяется от 0,999875 до 0,951732. Используя это свойство, организуем итерационный процесс следующим образом: изменяя величину М на
о
единицу и вычисляя значения Cvm и Cvd по формулам (15), (16) соответственно, добиваемся того, чтобы значение функции Dc (x) попало в указанный выше коридор. Таким образом, определяем величины M и xd■ Значения xe определяется путем решения следующего уравнения
(3N — M) Rx2e exp(xE)
Cvm =--------------------о----,
(exp(xE) — 1)
а начальное значение коэффициента. Гртопайзепаї из равенства. (14).
В качестве исходных данных для расчета начального значения коэффициента Гртопай-зетта и параметров тепловой части уравнения состояния тэтта используем данные, приведенные в работе [8J: N = 29,cp = 1088 Дж/(кг-К), cv = 1082, 36 Дж/(кг-К), ц = 316 кг/кмоль, а = 0, 2296 - 10_3 1/К, csq = 2320 м/с. В результате применения предложенного в работе алгоритма были получены следующие начальные значения для уравнения состояния кристаллического тэна: M = 23,xd = 0, 16530,$d = 48,43 К, xe = 4,12655, 0e = 1209,079 К,
Y = 2,03. Полученное начальное значение коэффициента Грюнайзена хорошо согласуется с известными экспериментальными данными для аналогичных молекулярных кристаллов [5].
Таким образом, в данной работе построена математическая модель pi подробно описан подход, позволяющие получать уравнения состояния для молекулярных кристаллов рі кррі-сталлріческріх взрывчатых веществ.
Литература
1. Жарков, В.Н. Уравпетшя состоят™я прті вьісокріх температурах pi давлетгаях / В.Н. Жарков, В.А. Калрпгап. - М.: Наука, 1968.
2. Ковалев, К).М. Уравпетшя состоятшя pi температуры ударного сжаттія кррісталлріческріх ВВ / Ю.М. Ковалев /'/ Фрізріка горетгая pi взрыва. - 1984. - Т. 20, № 2. - С. 102-107.
3. Молодец, А.М. Футікцрія Грюнайзена pi пулевая різотерма трех металлов до давлетгай 10 ГПа / A.M. Молодец// Журті. эксперршепт. рі теор. Фрізрікрі. - 1995. - Т. 107, № 3. -С. 824-831.
4. Жрірріфалько, Л. Статрістріческая фрізріка твердого тела/ Л. Жрірріфалько. - М.: Мтір, 1975.
5. Крітайгородскрій, А.И. Молекулярные кррісталльї/ А.И. Ктітайгородскрій. - М.: Наука, 1971.
6. Молодец, А.М. Футікцрія Грюнайзена, определенная па основе закономерностей ударноволнового сжаттія МОТІОЛРІТТЮГО матерріала / A.M. Молодец // Докл. РАН. - 1995. -Т. 341, .№ 6. - С. 753-754.
7. Кларк,Т. Компьютерная хрімрія / Т. Кларк. - М.: Мрір, 1990.
8. Dobrat.s, В.М. LLNL Explosives Handbook. Properties of Chemical Explosives and Explosive Simulants / B.M. Dobrat.s, P.C Crawford. - riverrriore, California: University of California, 1985.
Юррш Мріхайловріч Ковалев, доктор фрізріко-математріческріх паук, профессор, заведуто-щрій кафедрой вьічріслрітельпой мехатірікрі сплошных сред, Южпо-Уральскрш государственный утіріверсрітет (НИУ), (г. Челябрітіск, Росстійская Федерацрш), kov@mail.ru.
Bulletin of the South Ural State University. Series «Mathematical Modelling, Programming & Computer Software:»,
2013, vol. 6, no. 1, pp. 34-42.
MSC 76T25
Mathematical Modelling of the Thermal Component of the Equation of State of Molecular Crystals
Yu.M. Kovalev, South Ural State University, Chelyabinsk, Russian Federation, kov@inail.ru
This work is devoted to the construction of mathematical models of the equation of state of molecular crystals. Its practical value lies in the fact that all the solid explosives are molecular crystals. Therefore, after developing a mathematical model of the equation of state of molecular crystals, it will be possible to predict the behavior of explosives at high pressures and temperatures. The difficulty of the construction of the equation of state molecular crystals is that a large number of degrees of freedom of the molecules in the crystal structure does not allow the use of direct calculations. In this paper, we proposed an approach that allowed the use of all that is best in the Debye and Einstein models to describe the thermodynamics of crystals. Frequency separation of the normal vibrations in the crystal at high-frequency and low-frequency fluctuations allowed to getan analytical expression for the Griineisen and parameters for the thermal component of the equation of state of molecular crystal.
Keywords: equation of state, pressure, temperature, energy, Griineisen coefficient.
References
1. Zharkov V.N., Kalinin V.A. Equation of State at High Temperatures and Pressures. Moscow, Nauka, 1968. (in Russian)
2. Kovalev Yu.M. Equation of State and Temperature of Shock Compression of Crystalline Explosives. Combustion, Explosion, and Shock Waves, 1984, vol. 20, no. 2, pp. 102-107. (in Russian)
3. Molodet.s A.M. Griineisen Function and Zero Isotherm Three Metals to Pressures of 10 GPa. J. of Experimental and Theoretical Physics, 1995, vol. 107, no. 3, pp. 824-831.(in Russian)
4. ZhirifaPko L. Statistical Physics of Solids. Moscow, Mir, 1975. (in Russian)
5. Kitaygorodskiy A.I. Molecular Crystals. Moscow, Nauka, 1971.(in Russian)
6. Molodet.s A.M. Griineisen Function, Defined on the Basis of the Laws of Shock Compression of Solid Material. Doklady RAN, 1995, vol. 341, no. 6, pp.753-754. (in Russian)
7. Clark T. Computer Chemistry. Moscow, Mir, 1990. (in Russian)
8. Dobrat.s B.M., Crawford P.C LLNL Explosives Handbook. Properties of Chemical Explosives and Explosive Simulants. Livermore, California, L'niversit.y of California, 1985.
Поступила в редакцию 28 ноября 2012 г.
