Расчет параметров состояния реальных криогенных рабочих веществ

Расчет параметров состояния реальных криогенных рабочих веществ

Контрольная работа

РАСЧЕТ ПАРАМЕТРОВ СОСТОЯНИЯ РЕАЛЬНЫХ КРИОГЕННЫХ РАБОЧИХ ВЕЩЕСТВ

1. Уравнение состояния реальных веществ

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

Простое эмпирическое уравнение состояния реального газа, полученное с учетом сжимаемости 1 кг рабочего вещества, имеет вид

, (1.1)

где p - давление; v - удельный объем; z - коэффициент сжимаемости; R - газовая постоянная; T - температура. Широко применяются вириальные формы уравнения (1.1) и уточненные эмпирические зависимости с различными индивидуальными коэффициентами, в определенной мере учитывающими многообразие и сложность молекулярных взаимодействий и другие особенности свойств реальных рабочих веществ.

При обобщении экспериментальных данных часто используют вириальное уравнение Боголюбова-Майера [1], которое базируется на наиболее строгом теоретическом подходе и имеет вид

, (1.2)

где  - вириальные коэффициенты, зависящие от температуры;  - плотность.

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

. (1.3)

Каждый из вириальных коэффициентов, в свою очередь, записывают как

,

где s_i - границы изменения j; b_ij - коэффициенты разложения;  - приведенная температура, здесь T_кр - критическая температура.

Часто уравнение Боголюбова-Майера (1.3) представляют в виде [2]

. (1.4)

По мере накопления экспериментальных данных по термодинамическим свойствам рабочих веществ теплоэнергетических установок и систем в научно-технической литературе появилось множество уравнений состояния, описывающих с определенной полнотой и точностью параметры их состояния в p-v-T-пространстве. Наряду с уравнением Боголюбова-Майера получили применение уравнения Битте-Бриджмана, Бенедикта-Рубина [1], а в последние годы - уравнения Рабиновича, Сычева, Вассермана и др. [3-6], взаимосогласованные уравнения Клецкого, уравнения состояния гелия-4 Тарана, Мак-Карти и т. д.

. Расчет параметров состояния криогенных рабочих веществ

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

Единое уравнение состояния криогенных рабочих веществ. При разработке системы расчета термодинамических и теплофизических свойств криогенных рабочих веществ в основу было положено уравнение состояния Сычева-Вассермана, которое имеет единую форму для основных криоагентов и обеспечивает относительно высокую точность расчета p-v-T-данных [7].

При аппроксимации экспериментальных данных по основным криогенным веществам было получено уравнение состояния (1.4) в следующем виде:

, (1.5)

где - приведенная плотность, здесь ρ_кр - критическая плотность.

Для повышения точности описания p-v-T-данных водорода и гелия-4 используют два набора коэффициентов для определенных температурных зон. В связи с тем что характер изменения термических свойств гелия-4 в области температуры от λ-линии до 25 К значительно отличается от характера изменения их в диапазоне температур от 25 до 1500 К, одно уравнение действует в зоне температур от λ-линии до 25 К, а второе - при значениях T =-15…1500 К. Диапазон 15…1500 К частично охватывает каждое из двух уравнений, поэтому переход от одной группы коэффициентов разложения b_ij к другой осуществляется при температуре 20 К. Однако необходимо отметить, что уравнение состояния в вириальной форме (1.5) не отображает достоверно поведение термодинамических функций в области критической точки.

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

, (1.6)

где f_i - коэффициенты уравнения; - относительная температура.

Калорические уравнения для криогенных рабочих веществ. Калорические уравнения, связывающие энтальпию, энтропию и теплоемкость с давлением и температурой, могут быть получены на основе обобщения экспериментальных данных либо путем интегрирования p-v-T-зависимости. В работах [3-6, 9] используется второй путь и предлагаются формулы для определения следующих значений: энтальпии

; (1.7)

энтропии

; (1.8)

изобарной теплоемкости

, (1.9)

где  - энтальпия, энтропия и изохорная теплоемкость в идеально-газовом состоянии.

Энтальпию h₀ и энтропию s₀ определяют из соотношений

; (1.10)

, (1.11)

где h₀₀ и s₀₀ - энтальпия и энтропия при температуре T₀;  - теплота сублимации при T1=10 К;  - константа.

Значения констант уравнений (1.7)-(1.11) приводятся в таблице.

Изобарную теплоемкость в идеально-газовом состоянии рассчитывают с помощью полинома

, (1.12)

где α _j и β _j - коэффициенты полинома.

При аппроксимации значений изобарной теплоемкости в идеально-газовом состоянии для азота в полиноме (1.12) вместо Θ_T используется значение приведенной температуры .

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

. (1.13)

Параметры криогенных веществ и константы уравнений (1.7)-(1.11) для T₀ = 100 K [3-6, 9]

криогенный газовый энтропия смесь

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

На базе уравнений (1.5)-(1.13) разработана единая подсистема расчета термодинамических и калорических параметров основных рабочих веществ [10, 11] криогенных установок и систем. Значения параметров этих веществ можно рассчитывать с помощью отдельных подпрограмм-функций, структурно входящих в данную подсистему.

Коэффициент сжимаемости z и давление p реального криогенного рабочего вещества по известной плотности ρ и температуре T определяем с помощью подпрограмм, изображенных на рис. 1.1 и 1.2.

Рис. 1.1. Подпрограмма вычисления коэффициента сжимаемости

В подпрограммах на рис. 1.1 и 1.2 приняты следующие идентификаторы: W - приведенная плотность; RO, ROKR - плотность и критическая плотность; TAU - приведенная температура; Т, TKR - температура и критическая температура; Р, PKR - давление и критическое давление; R - специфическая газовая постоянная; B(I,J) - коэффициенты разложения уравнения (1.5).

Рис. 1.2. Подпрограмма вычисления давления

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

Рис. 1.3. Вспомогательная подпрограмма для определения плотности рабочего вещества

Подпрограмма на рис. 1.3, в свою очередь, используется в основной расчетной подпрограмме CRRO(P,T,N) при определении плотности рабочего вещества по заданным значениям давления и температуры.

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

Рис. 1.4. Подпрограмма расчета плотности рабочего вещества в зависимости от фазового состояния

Как видим, сходимость решения осуществляется по заданной точности EPS, а число приближений задается величиной KN. Область состояния рабочего вещества задается значением переменной N: для жидкости N = 0 , для парогазовой смеси N =1. В соответствии с этим подпрограмма (см. рис. 1.4) для области жидкости записывается как CRRO(P,T,0), а для парогазовой части p-v-T-пространства - как CRRO(P,T,1).

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

Рис. 1.5. Универсальная подпрограмма расчета плотности по давлению и температуре

На рис. 1.5 FTS(P) - подпрограмма расчета температуры насыщения по давлению.

Энтальпию реального рабочего вещества определяем по плотности и температуре [см. уравнение (1.7)]. Подпрограмма расчета энтальпии изображена на рис. 1.6.

Рис. 1.6. Подпрограмма расчета энтальпии реального рабочего вещества

Для вычисления энтальпии рабочего вещества используем две вспомогательные подпрограммы: СР0(Т) - для определения изобарной теплоемкости [см. уравнение (1.12)] и Н0(Т) - для определения энтальпии [см. уравнение (1.10)] в идеально-газовом состоянии. Эти подпрограммы показаны на рис. 1.7 и 1.8.

Рис. 1.7. Подпрограмма расчета изобарной теплоемкости

Рис. 1.8. Подпрограмма расчета энтальпии в идеально-газовом состоянии

Для вычисления определенного интеграла в уравнениях (1.10) и (1.11) применяем стандартную процедуру GAUSS. В подпрограммах (см. рис. 1.6-1.8) приняты следующие идентификаторы: H0 - энтальпия в идеально-газовом состоянии; СР0 - изобарная теплоемкость в идеально-газовом состоянии; Н00 - энтальпия при T₀ (см. таблицу); Н000 - теплота сублимации при T = 0 К.

Энтропию, согласно подпрограмме на рис. 1.9, рассчитываем подобно энтальпии (см. рис. 1.6) по плотности и температуре. При этом используем две вспомогательные подпрограммы (рис. 1.10 и 1.11): CP0T(T), в которых определяем значение подынтегрального отношения СР0/Т, и S0(T) - для вычисления энтропии в идеально-газовом состоянии [см. уравнение (1.11)].

Рис. 1.9. Подпрограмма расчета энтропии рабочего вещества

Рис. 1.10. Вспомогательная подпрограмма для вычисления энтропии

Рис. 1.11. Подпрограмма для расчета энтропии в идеально-газовом состоянии

В подпрограммах (см. рис. 1.9-1.11) использованы следующие идентификаторы: S0 - энтропия в идеально-газовом состоянии; S00 - энтропия при T₀; S000 - константа уравнения (1.11).

Подпрограммы вычисления теплоемкости при постоянном давлении или объеме по плотности и температуре (рис. 1.12 и 1.13) разработаны на базе соотношения (1.9).

Рис. 1.12. Подпрограмма расчета теплоемкости при постоянном давлении

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

Рис. 1.13. Подпрограмма расчета теплоемкости при постоянном объеме

Рис. 1.14. Подпрограмма расчета давления насыщенных паров

На рис. 1.14 F(I) - коэффициент аппроксимации f_i в уравнении (1.6).

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

Рис. 1.15. Подпрограмма расчета температуры насыщения

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

Рис. 1.16. Подпрограмма расчета теплоты парообразования

При решении ряда задач гидродинамики важно знать скорость звука. Для ее вычисления используем подпрограмму, разработанную на базе соотношений, приведенных в [5] (рис. 1.17). В подпрограмме на рис. 1.17 VSO - скорость звука

.

Рис. 1.17. Подпрограмма расчета скорости звука

При моделировании тепловых и гидродинамических процессов в каналах проточных трактов элементов систем криогенного обеспечения объектов криостатирования часто возникает необходимость в определении калорических параметров рабочих веществ и скорости звука по известным значениям давления и температуры. Для этого на базе подпрограмм (см. рис. 1.5, 1.6, 1.9, 1.12, 1.13 и 1.17) разработаны соответствующие расчетные подпрограммы (рис. 1.18-1.22).

Рис. 1.18. Подпрограмма расчета энтальпии по известным значениям давления и температуры

Базовые подпрограммы-функции (см. рис. 1.1-1.22) структурно входят в основное тело программы подсистемы расчета термических и калорических параметров криогенных рабочих веществ. Массивы коэффициентов расчетных уравнений и характерные параметры для каждого криоагента (см. таблицу) описываются в подпрограммах SUBROUTINE N2 для азота, SUBROUTINE О2 для кислорода и в аналогичных подпрограммах для других газов операторами DIMENSION и DATA, а обмен информацией осуществляется с помощью соответствующего оператора СОМMON. Подпрограммы-функции использованы при разработке программы расчета параметров реальных тепловых и гидродинамических процессов в элементах криогенных систем и установок.

Рис. 1.19. Подпрограмма расчета энтропии по известным значениям давления и температуры

Рис. 1.20. Подпрограмма расчета теплоемкости при постоянном давлении по известным значениям давления и температуры

Рис. 1.21. Подпрограмма расчета теплоемкости при постоянном объеме по известным значениям давления и температуры

Рис. 1.22. Подпрограмма расчета скорости звука по известным значениям давления и температуры

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

Для расчета параметров в двухфазной области применяем процедуру FAZA2(P), которая базируется на уравнении состояния (1.5) (рис. 1.23).

Рис. 1.23. Подпрограмма расчета параметров в двухфазной области

В начале решения находим температуру кипения, применяя подпрограмму CRPS(T) (см. рис. 1.14). Затем для локализации корней уравнения состояния, лежащих на левой и правой пограничных кривых области двухфазного состояния, в процедуре CRRO(P,T,N) принимаем соответствующие значения для N (см. рис. 1.4). Далее вычисляем энтальпию, энтропию жидкости и пара и теплоту парообразования.

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

Рис. 1.24. Схема расчета параметров в конце процесса дросселирования:

В соответствии с этим предусмотрены три процедуры.

.        Процедура DROSG(P1,Т1,Р2), с помощью которой по начальным давлению и температуре и конечному давлению находим конечную температуру газовой фазы. В данной процедуре производится определение энтальпии в точке, характеризующей начало процесса дросселирования; далее методом половинного деления решается уравнение h_н - h_к = 0, из которого определяется температура конца дросселирования. Программа этой процедуры показана на рис. 1.25.

.        Процедура DROS(P1,Т1,Р2) служит для расчета параметров газа в конце процесса дросселирования в двухфазной области (рис. 1.26). В этой процедуре по конечному давлению рассчитывается температура кипения, которая и является температурой конца процесса дросселирования.

По процедуре FAZA2(P) находим параметры двухфазной области, степень сухости и соотношение между потоками пара на входе в дроссель и на выходе из него.

Рис. 1.25. Подпрограмма расчета процесса дросселирования, заканчивающегося в газовой фазе

Рис. 1.26. Подпрограмма расчета процесса дросселирования, заканчивающегося в двухфазной области

3.      Процедура DROSL(P1,T1,P2) (рис. 1.27) служит для вычисления параметров в конце процесса дросселирования насыщенной жидкости, так как в установках часто осуществляется данный процесс (процесс 1¢¢-2¢¢ на рис. 1.23), что приводит к изменению энергетического состояния криоагента после дросселя.

Рис. 1.27. Подпрограмма расчета процесса дросселирования насыщенной жидкости

На рис. 1.28 показаны процессы расширения рабочего вещества в парожидкостном детандере (1-2) и изоэнтропного сжатия (1¢-2¢). В обоих случаях стоит задача определения параметров состояния рабочего вещества в конце процессов для случая, когда заданы начальные давление р_н и температура Т_н, конечное давление р_к и изоэнтропный КПД h_s.

Рис. 1.28. Схема расчета параметров в конце процессов расширения в детандере (1-2) и изоэнтропного сжатия

В первом случае по известным параметрам р_н и Т_н находим энтальпию h_н и энтропию s_н, а по значению р_к с помощью процедуры FAZA2(PK) - температуру кипения Т_кип, энтальпию h_L и h_G, энтропию s_L и s_G на пограничных кривых области двухфазного состояния. Далее по известным значениям s_н, s_L, s_G определяем положение точки 2_s, а по заданному значению КПД (CPD) - истинное состояние рабочего вещества в конце расширения (рис. 1.29).

Рис. 1.29. Подпрограмма расчета параметров после расширения в детандере

Для определения состояния рабочего вещества в конце изоэнтропного сжатия применяем процедуру SCONST(PH,TH,PK,K) - рис. 1.30.

Рис. 1.30. Подпрограмма расчета параметров после изоэнтропного сжатия

По этой процедуре находим значения энтальпии и энтропии в точке 1¢ и определяем температуру

.

Поскольку показатель k в зависимости от области состояния может отличаться от принятого первоначального значения, найденное значение температуры Т_s является первым приближением, которое в последующем уточняется в результате решения методом половинного деления уравнения s_н - s_к = 0.

Часто при проведении вычислительных экспериментов возникает необходимость расчета температуры рабочего вещества по известным значениям давления и энтальпии, для этого может быть использована процедура TH(P,H,N) - рис. 1.31.

Рис. 1.31. Подпрограмма расчета температуры рабочего вещества по известным значениям давления и энтальпии

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

. (1.14)

Реализация решения полинома (1.14) приводится в процедуре FLA(RO,T) - рис. 1.32.

Рис. 1.32. Подпрограмма расчета теплопроводности

На рис. 1.32 FLA - теплопроводность, Вт/(м×К).

Динамическую вязкость рассчитываем по следующему уравнению:

. (1.15)

Слагаемые уравнения (1.15) вычисляем по формулам

; (1.16)

Для определения динамической вязкости по уравнениям (1.15)-(1.17) в программе предусмотрена процедура FMU(RO,T) - рис. 1.33.

Рис. 1.33. Подпрограмма расчета динамической вязкости

На рис. 1.33 FMU - динамическая вязкость, Па×с.

. Расчет параметров состояния смеси

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

Уравнение Редлиха-Квонга для смесей в модификации Соаве имеет вид [12]

, (1.18)

где z_m - коэффициент сжимаемости для смеси; J_m - молярный объем смеси, м³/кмоль; WA и WB - простые числа (WA = 0,4274802327; WB = 0,086640350);

, (1.19)

здесь y_j - молярная доля j-го компонента;

; (1.20)

M_j - молярная масса компонента смеси.

Коэффициент F_m учитывает параметры бинарных взаимодействий компонентов смеси и, согласно модификации Соаве, определяется как

, (1.21)

где i, j (i-й и j-й) - парные компоненты смеси;  - параметр бинарного взаимодействия компонентов; F_i, F_j - коэффициенты компонентов, которые определяются по уравнению

, (1.22)

здесь  - фактор ацентричности молекул вещества.

На базе алгоритма (1.1) и выражений (1.18)-(1.22) разработаны блок-схемы расчета коэффициента сжимаемости, давления и плотности смеси.

На рис. 1.34 показана блок-схема расчета коэффициентов уравнения (1.18).

Рис. 1.34. Блок-схема алгоритма расчета коэффициентов уравнения (1.18)

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

.

В блоке 1 осуществляется ввод исходной информации, характеризующей состав и свойства компонентов смеси.

В блоке 2 рассчитываются индивидуальные коэффициенты F_i (F(I)) по уравнению (1.22), b_i (В(I)) по (1.20), а также суммарный коэффициент b_m (ВМ) по (1.19).

В блоках 4 и 5, которые структурно входят в циклы по I и по J, рассчитываются знаменатель и числитель уравнения (1.21), а в блоке 8 завершается его решение определением коэффициента F_m (FM), учитывающего параметры бинарных взаимодействий компонентов смеси.

Блоки 3, 6 и 7 организуют расчетные циклы по I и J по числу компонентов смеси от 1 до N.

Вывод и передачу рассчитанных величин в соответствующие подпрограммы и процедуры осуществляет блок 9.

На основании приведенной блок-схемы разработана и реализована подпрограмма MIXTURE.

Для расчета давления смеси [см. уравнение (1.1)] по известной плотности и температуре может быть использована блок-схема (рис. 1.35).

Рис. 1.35. Блок-схема алгоритма расчета давления

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

В блоке 2 по уравнению (1.18) определяется значение коэффициента сжимаемости смеси, а в блоке 3 решается уравнение состояния реального вещества относительно давления.

Вычисление плотности с использованием уравнения (1.18) производится в соответствии с блок-схемой (рис. 1.36).

Рис. 1.36. Блок-схема алгоритма расчета плотности смеси

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

. (1.23)

Для локализации корней уравнения состояния принимаем первоначальные значения плотности пара и жидкости, которые идентифицируются целочисленной переменной М: для жидкости М = 0, для пара М = 1 (см. рис. 1.36, блок 2).

В блоках 4 и 6 обеспечивается реализация метода Ньютона, для чего при заданных возможных значениях плотности при постоянной температуре с помощью уравнения (1.23) находим частную производную , уточненное значение плотности и степень сходимости DELTA решения уравнения (1.23).

После сравнения (в блоке 7) DELTA с заданной точностью решения EPS принимаем решение либо об обращении в блок 3 для уточнения исходного значения плотности и повторения расчетов, либо, при достижении необходимой сходимости, о передаче информации в блок 8 вывода информации, на чем решение задачи заканчивается.

Литература

1.Архаров A.M., Марфенина И.В., Микулин Е.И. Криогенные системы. Т. 1. Основы теории и расчета. - М.: Машиностроение, 2006. - 576 с.

.Теплофизические свойства неона, аргона, криптона и ксенона /Под ред. В.А. Рабиновича. - М.: Изд-во стандартов, 1976. - 636 с.

.Термодинамические свойства азота /В.В. Сычев, А.А. Вассерман, А.Д. Козлов и др. - М.: Изд-во стандартов, 2007. - 352 с.

.Термодинамические свойства воздуха /В.В. Сычев, А.А. Вассерман, А.Д. Козлов и др. - М.: Изд-во стандартов,2008. - 346 с.

.Термодинамические свойства гелия-4 /В.В. Сычев, А.А. Вассерман, А.Д. Козлов и др. - М.: Изд-во стандартов,2004. - 307 с.

.Термодинамические свойства кислорода /В.В. Сычев, А.А. Вассерман, А.Д. Козлов и др. - М.: Изд-во стандартов,2001. - 312 с.

.Теплофизические свойства криопродуктов /Л.А. Акулов, Е.И. Борзенко, В.Н. Новотельнов, А.В. Зайцев. - СПб.: Политехника, 2001. - 243 с.

.Справочник по физико-техническим основам криогеники /М.П. Малков, И.Б. Данилов, Л.Г. Зельдович, А.Б. Фрадков; Под ред. М.П. Малкова. - 3-е изд., перераб. и доп. - М.: Энергоатомиздат, 2005. - 432 с.