Математическое моделирование тепловой составляющей уравнения состояния молекулярных кристаллов

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

В настоящее время принято считать, что в уравнения состояния молекулярных кристаллов входит две составляющие: тепловая и холодная. Тепловая составляющая определяется колебательным движением молекул, входящих в состав кристалла, а. холодная составляющая - изменением энергии взаимодействия внутри молекулы и между молекулами. Определение зависимости коэффициента Грюнайзена у от удельного объема является одной из основных задач при построении уравнений состояния твердых тел. Это обусловлено тем, что данная зависимость является связующей между тепловой и холодной составляющими уравнения состояния в соответствии с формулами Ландау - Слейтера, Дугдала - Мак-Дональда, и т.д. [1]. Решению этой задачи посвящено достаточно большое количество как экспериментальных, так и теоретических работ [1-3]. Однако, теоретическое определение зависимостей, характеризующих поведение твердых взрывчатых веществ (ВВ), осложняется тем, что они относятся к молекулярным кристаллам, а. молекулы, входящие в состав кристалла, обладают большим числом внутренних степеней свободы. Тем не менее, и в этом случае можно получить аналитическую зависимость для коэффициента. Грюнайзена. от удельного объема, твердого тела.

Термодинамические свойства, вещества, полностью определяются, если известен один из термодинамических потенциалов. Удобно исходить из определения свободной энергии Гельмгольца F ( V, T ), которая наиболее простым образом связана с моделью строения вещества:

F = U + E о + kT £ ln(1 - exp( -w)), E о = 1 ^ hw_a.           (1)

Z-^              kl           2 z-^

αα

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ

Здесь U - энергия взаимодействия между атомами: V - хмельный объем: Т - температура, тела: к - постоянная Планка; w _a - частоты норма.тьных колебаний: E о - энергия нулевых колебаний.

Если F ( V, Т ) задана, то дифференцированием могут быть найдены выражения для всех измеряемых термодинамических величин:

S

-

∂F          ∂F     ∂S      ∂P

V дт)V’        V dVT_T ’ V dVТ_т ддТ)_v•

Известно, что энергия взаимодействия между атомами молекулярного кристалла, складывается из внутримолекулярной и межмолекулярной. Внутримолекулярная энергия состоит из энергии валентных взаимодействий и энергии невалентных взаимодействий атомов внутри молекулы. Межмолекулярная энергия является энергией не валентных взаимодействий атомов. В силу того, что внутримолекулярная энергия валентных взаимодействий существенно больше энергии невалентных взаимодействий между атомами как внутри молекулы, так и разных молекул, имеет смысл разделить энергию взаимодействия молекулярного кристалла на молекулярную U m (энергия валентных взаимодействий) и кристаллическую U k (энергия невалентных взаимодействий). Если энергия U k зависит от пространственного расположения молекул, т.е. от объема, то энергия U m зависит исключительно от величины валентных связей и валентных углов.

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

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

θTD

F - U k + U m + E о + 3 MRT ( T V ^J ^ ² ln(1 - exp( - )) d^ +

D ₀

+ (3N - M) RTIn (1 - exp (-6-E^ } , где R - универсальная газовая постоянная: M - число низкочастотных колебаний: N - число атомов в молекуле; 3N - M - число высокочастотных колебаний; 6d - характеристическая температура Дебая; 6е - характеристическая температура Эйнштейна.

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

Интегрируя по частям выражение для низкочастотной составляющей свободной энергии F ( V, Т ) (2) и вводя функцию Дебая D ( x ) по формуле приведенной в монографии [4, 5],

x dξ exp(С) - 1’

D ( x ) - | / , 3

получаем

F = U m + U k + MRT (ln(1 — exp( — x d )) — D ⁽ ^D ) ) +(3 N — M ) RT ln(1 — exp( — x e )) , (3)

где x d = ^ D; x e = T E.

Используя выражение (3), легко получить выражения для давления P и энтропии S;

P = —

∂V _T

-

∂U M   ∂U K   dE 0

-

--

∂V ∂V dV

MRTD ( x d ) Д ТПД Г d (ln V ) V

-

— (3 N — M ) RTx e

d (ln 9 e )

d (In V ) V (exp ( x e ) — 1) ’

S = ^—(IT) v = ^—{MR h¹

— exp ( — x d )) —

D ( X D )

+

+ MRT

exp ( — x d )      D'

1 — exp ( — x d )

-

( X D )! dx D

дТ ⁺

+ (3 N — M ) R ln(1 — exp ( — x e )) + (3 N — M ) RT-

exp ( — x e ) dxE^ exp ( x e ) дТ J

-

— ^MR |^ln (1 — exp ( — x d ))-- D^ j — MRD ( x d ) + (3 N — M ) R ln (1 — exp ( — x e )) —

(3 N — M )   R^ V

exp( x e ) — 1

При выводе формулы (5) было использовано свойство функции Дебая

D (x ) =----— xD' (x), exp (x) — 1   3

где штрих обозначает дифференцирование по характеристической температуре x.

Зная равенства (3) и (5), легко определить выражения для полной энергии E и емкости при постоянном объеме Су;

тепло-

E = F + TS = U m + U k + E о + MRTD ( x d ) + (3 N — M ) — R^TxE ₁. ; exp ( x e — 1)

C_v = MR 4DD ( x_D )--3^) + (3 N — M ) Rx2_E exp ^{( x}E ) ₂ .

V        V ^V D    exp ( x d ) — 1) ^V         ’ ^E (exp ( x e ) — 1) 2

Следуя определению коэффициента Грюнайзена

Y D ( V )

-

d (In 9 d ) d (In V ) ,

выражение (4) можно записать в виде

P = —

∂U M   ∂U K

-

∂V ∂V

-

dE о   MRT y d ( V ) D ( x d )

dV ⁺       V

.

Последний член в выражении (4) равен нулю, так как при разделении частот было сделано предположение о независимости высоких частот от объема.

Исходя из определения энергии нулевых колебаний и учитывая разделение частот, получаем выражения для функций E о и dE ⁰:

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ

E о = 2 £ hu a = α

3 N' - М', hw _α

w D

+ ^М j„ 2 hwdw =

³ N - M R9 e + 3 MR9 d ( V ); 2            8

dE о _   3 MRy d ( V ) 9 d ( V )

dF = - 8      V где N' - число атомов. M' - число шгзкочастотш>ix колебании в объеме V.

Подставляя выражение для производной от энергии нулевых колебаний по объему (7) в равенство (6), получаем уравнение для определения давления в виде n MRTyd (V) /3         А n n    dUK

P = V        x ^{D +} D ⁽ x ^D )J ⁺ P y , P y = ^— dV •              ⁽⁸⁾

Для получения возможной зависимости коэффициента Грюнайзена для молекулярных кристаллов от объема ВВ может быть использован следующий подход.

Изотермический модуль сжатия в т (изотермическая сжимаемость) определяется выражением

1    с T вТ = V

∂P

^V ∂V T ^,

где с т обозначает изотермическую скорость звука.

Подставляя выражение (8) в правую часть равенства (9), получаем

1       д \MTTdd(V) /3                         MRTyd (V) /3       . А вт = -VdV [ V     UxD + D (xD)J + Py] T - -V [    V2    (sxD + D (xD)J +

MRTy в ( V ) /3         A MRTyd ( V ) /3          dx_D dP

+ ----V      (sxD + D (xD)) + ----V (s + D (xD)) "VV + дЯ T = т Г MRTyd (V) /3         A

= ^-V      V 2     ^8 x D ⁺ D ⁽ x D ) J ⁺

MRTy d ( V ) /3         A mrty D ( V ) /3                        dP y l

⁺-- V ---- ( 8 x D ⁺ D ⁽ x D ) )-- V^--- ( 8 x D — C VD ⁽ x D ) ⁺ D ⁽ x D ) ) ⁺ dV =

= MVT ( y d ( V ) + y d ( V )) (I X D + D ( X D )) - V⁸A - MRRT D yp DAx D l -

- MRTy' d

^{( V} ) (8 x D ⁺ D

( X D ) ,

где y D ( V ) ~ производная по V от коэффициента Грюнайзена;

x D

^- vy d ⁽ V ) •

C VD ( x D ) = 4 D ( x D ) - _exp( X D - 1_; dV

Определим значение изотермической сжимаемости при Т ^ 0. С этой целью переходим в последнем выражении к пределу:

, .   . .                 . ..   3                3                       , .   . .

ЙШ0 MRpT ( y D ( V ) + y d ( V )) 8 X D = 8 9 d MRp ( y D ( V ) + y d ( V )) ;

lim

T м 0

MRpT (D_D ( V ) + y d ( V )) D ( x d ) = 0;

jim MRpTy D ( V ) C vd ( x d ) = 0;

3      3

lim MRTY d ( V ) - x d = - MR9 d Y d ( V );

T ^ 0               о о

Um MRTy D ( V ) D ( x d ) = 0 .

Следовательно, при T ^ 0 изотермическая сжимаемость определяется выражением вида

/ 1 \           3             , .                    ,     3.

= я 9 d MRP Y D ( V ) + Y D ( V )) - MR9 d Y d ( V ) - V^ .

\ eT ) т ^ о   8                          8

По определению изотермическая сжимаемость в т связана с изотермической скоростью звука, с т формулой (9). т.е.

-T = -9dMRp (yD (V) + yd (V)) - 5MR9dyD (V) - V—У • V о                           о

Полагая далее, что скорость звука, определяется только упругими свойствами кристаллов, получаем следующее дифференциальное уравнение типа, уравнения Бернулли для определения зависимости коэффициента. Грюнайзена от плотности:

Y d ^{( V} ) ^- VY D ^{( V} ) ^- VY D ⁽ V ) = ⁰ .

Данное уравнение заменой z = y d ( v ) приводится к линейному дифференциальному уравнению

z z + V =

V,

которое легко интегрируется. Зависимость коэффициента. Грюнайзена. от плотности описывается выражением следующего вида:

V

Yd = C - V ’ где константа C определяется из условия yd (Vo) = Y^d- Похожее выражение для коэффициента. Грюнайзена. было получено из других соображений в работах А.М. Молодца. [3, 6].

В результате получим

^{Y D} = ^{y D} ® Г+Дф г у •

Формула. (10) в случае слабого сжатия переходит в известное эмпирическое выражение

YD=yD (у) '

которое широко используется при обработке экспериментальных.

Таким образом, предположение о том, что изотермическая скорость звука при T = 0 определяется только упругими свойствами кристаллов, позволило получить аналитическую зависимость коэффициента. Грюнайзена. от объема.

Для того чтобы воспользоваться выражением для коэффициента Грюнайзена y d (Ю), необходимо определить его значение yD, соответствующее начальному удельному объему

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ

V 0. С этой целью получим аналог коэффициента Грюнайзена для молекулярных кристаллов. Для этого раскроем выражение отношения коэффициента теплового расширения а к изотермическому модулю объемного сжатия в т [4]

а - - m 1 -           di»

Подставляя в равенство (11) выражение (3) и дифференцируя последовательно по T и V, получим вт=-tv M tn с- exp (-т))] - зD (т)+

+ MRT

( ^exP ( D ) - 1 )

- ¹ D

=

T ∂T T

= Y D ( V о ) Р ⁰ MR

4 D ( x d ) -

3 x D

(exp( x d ) - 1)

Вводя обозначение C vd для выражения, стоящего в квадратных скобках в правой части равенства (12), получим аналог уравнения Грюнайзена в традиционной форме:

а _ yd ( V ) C vd в т = V

Учитывая связь между изотермическим модулем сжатия в т ^и изотермической скоростью звука, С т

1 _ ст вт = V равенство (13) запишется в виде аСт = yd (V) Cvd • (14)

Следовательно, если из эксперимента, известны значения коэффициента, объемного расширения а и изотермической скорости звука С т и определен способ нахождения C vd, из уравнения (14) можно найти значение коэффициента. Грюнайзена. при начальном значении удельного объема. Это дает возможность использования выражения (10) для определения зависимостей коэффициента. Грюнайзена. молекулярных кристаллов от удельного объема.

При определении теплоемкости при постоянном объеме, относящейся к деформационной части, проведем анализ частот внутримолекулярных колебаний атомов на. примере кристалла тэна (С5H8N4O12). Силовые постоянные для расчета спектров нормальных колебаний были определены с помощью квантово-химического метода. РМ-3, подробно описанного в монографии [7]. Этот метод был выбран по тому, дает хорошее совпадение с экспериментальными данными, приведенными в работе [8]. Автор выражает благодарность профессору А.В. Белику за. предоставленные ИК - спектры для тэна, приведенью в таблице. Анализ конформации молекулы тэна. показывает, что в качестве колебаний, не зависящих от кристаллического окружения, с уверенностью можно выбрать 28 колебаний валентных связей, 4 колебания валентных углов нитрогрупп и 3 колебания валентных углов центрального атома, углерода.

Следовательно, начальное значение величины (3 N - M ) для молекулы тэна равно 35. Зная частоты нормальных колебаний и число колебаний, относящихся к не деформационной части (3 N - M ) , легко определить величину С ум по формуле

Таблица

3 N

Cvm = R £ i=M +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 11 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 x2exp(xi) (exp (Xi) - 1)2 ’

где X i = hkT i’ а значения w i для тэна выбираются из таблицы соответственно. Зная величины C vm ^и Cv ’ легко определить значение теплоемкости C vd, относящуюся к деформационным колебаниям

C vd = C v - C vm = MR ^4 D ( xD )-- -—D—- ’

V          exp( x d ) - 1/

Легко показать, что выражение, стоящее в скобках правой части формулы (16), есть результат интегрирования по частям выражения D e ( x )

x

3 f 44    ^exp⁽ С )    ,

DC (X)   x3J ^ (exp (e) - 1) dX’ которое представляет собой так называемую функцию теплоемкости Дебая. Данная функция при изменении xd от 0,05 до 1,0 изменяется от 0,999875 до 0,951732. Используя это свойство, организуем итерационный процесс следующим образом: изменяя величину M на

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ единицу и вычисляя значения Cvm и Cvd по формулам (15), (16) соответственно, добиваемся того, чтобы значение функции De (x) попало в указанный выше коридор. Таким образом, определяем величины M и xd. Значения хе определяется путем решения следующего уравнения

(3N - M ) Rx2E exp(xE ) CVM= (exp(xE) - 1)2 , а начальное значение коэффициента Грюнайзена из равенства (14).

В качестве исходных данных для расчета начального значения коэффициента Грюнай-зена и параметров тепловой части уравнения состояния тэна используем данные, приведенные в работе [8]: N = 29 , c p = 1088 Дж (кг - К). c_v = 1082 , 36 Дж (кг - К). ц = 316 кг кмоль. а = 0 , 2296 - 10 ~ ³ 1/К, c s о = 2320 м/с. В результате применения предложенного в работе алгоритма были получены следующие начальные значения для уравнения состояния кристаллического тэна: M = 23 ,x d = 0 , 16530 ,9 d = 48 , 43 К. х е = 4.12655. 9 е = 1209.079 К. Y = 2,03. Полученное начальное значение коэффициента Грюнайзена хорошо согласуется с известными экспериментальными данными для аналогичных молекулярных кристаллов [5].

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