ХИМИЧЕСКАЯ ФИЗИКА, 2014, том 33, № 1, с. 20-24
ГОРЕНИЕ, ВЗРЫВ ^^^^^^^^^^^^ И УДАРНЫЕ ВОЛНЫ
УДК 635.715
УРАВНЕНИЕ СОСТОЯНИЯ СВЕРХКРИТИЧЕСКОГО ФЛЮИДА НА ОСНОВЕ УРАВНЕНИЯ ОРНШТЕЙНА-ЦЕРНИКЕ © 2014 г. А. А. Аникеев*, С. Б. Викторов, С. А. Губин
Национальный исследовательский ядерный университет "МИФИ" *Е-таП: anikeev_aa@mail.ru Поступила в редакцию 27.04.2012; после доработки 20.09.2013
Выполнено численное моделирование термодинамических свойств флюида методом интегральных уравнений. Проведено сравнение полученных данных с результатами МК- и МД-моделирования. Изучена проблема устойчивости численного решения. Предложены методы коррекции корреляционных функций и анализа погрешностей их расчета.
Ключевые слова: уравнение состояния, сверхкритическое состояние, флюид, корреляционные функции, функция распределения, уравнение Орнштейна—Цернике, интегральные уравнения, высокие давления и температуры.
Б01: 10.7868/80207401X14010038
Для расчетов давления и избыточной внутренней энергии газообразных систем при высоких плотностях и температурах необходимо знание уравнения состояния вещества (УРС) в диапазоне сверхкритической газожидкостной фазы, которую принято называть сверхкритическим флюидом (СКФ). Термодинамическое моделирование на основе УРС необходимо при решении задач по физике и химии ударных и детонационных волн [1], астрофизических задач о внутренних слоях планет-гигантов и многих других. Использование эмпирических уравнений не позволяет правильно прогнозировать параметры состояния за пределами применимости экспериментальной области определения таких УРС [2]. Поэтому в современном термодинамическом моделировании все больше используются теоретические УРС, основанные на методах статистической механики и потенциалах взаимодействия молекул.
Существует два подхода для расчета параметров состояния. Первый — проведение численных экспериментов методами Монте-Карло (МК) или молекулярной динамики (МД), основанных на получении термодинамических макропараметров из результатов моделирования движения отдельных микрочастиц. Такой подход позволяет получить теоретически обоснованные данные с точно учитываемой статистической погрешностью. Однако прямое моделирование систем требует огромных вычислительных ресурсов и принципиально не может быть проведено за короткое время. Эти недостатки предопределили основное применение методов МК и МД как источников
референсных данных, по которым проводится калибровка новых УРС.
Второй подход — нахождение новых УРС на основе методов термодинамической теории возмущений (ТТВ) и интегральных уравнений функций распределения молекул (ИУР). Термодинамические расчеты с использованием этих УРС проводятся намного быстрее расчетов на основе методов МК и МД. Чаще всего для флюида применяются приближение сферически-симметричных частиц и модель ближнего порядка.
Метод ТТВ позволяет получать результаты с хорошей точностью вплоть до давлений порядка миллионов атмосфер для однокомпонентного флюида. При этом для расчетов может использоваться персональный компьютер. Однако при попытке расширить модель на многокомпонентные составы возникают значительные сложности [2]. Тем не мене на основе метода ТТВ построены современные термодинамические коды [2—7].
Кроме метода ТТВ, была построена модель для метода ИУР [8], который имеет существенные отличия от ТТВ. С одной стороны, этот метод требует больше машинного времени и не позволяет рассчитать свободную энергию системы. С другой стороны, теория интегральных уравнений проще распространяется на многокомпонентные составы. В ней не используются эффективный радиус твердых сфер и другие понятия модели ТТВ, не имеющие однозначной интерпретации для многокомпонентных составов.
g
Рис. 1. Радиальная функция распределения. Увеличено изображение шумов на плато.
Цель настоящей работы заключалась в реализации собственной модели уравнения состояния на основе ИУР и оценке возможностей этой модели.
МОДЕЛЬ БЛИЖНЕГО ПОРЯДКА
Рассмотрим макроскопическую систему, состоящую из большого числа одинаковых микрочастиц, взаимодействующих между собой. Зададим парный потенциал взаимодействия ф(г), где г — расстояние между частицами. Проведем осреднение положения частиц по времени и введем парную функцию распределения частиц:
/ = ЛЫ/^кОъ), Г12 = И - /2|,
где /1(г) — одночастичная функция распределения, g(г) — радиальная функция распределения (рис. 1).
Значение радиальной функции распределения при малых радиусах близко к нулю, поскольку вероятность того, что две частицы находятся очень близко друг к другу, мала вследствие сил отталкивания. Вблизи минимума потенциала взаимодействия нахождение частиц наиболее вероятно, что выражается в появлении максимума функции g(r). Дальнейшее поведение радиальной функции распределения отражает характер упаковки частиц. При больших расстояниях ближний порядок пропадает и положение частиц становится независимым ^(г) = 1), т.е. их координаты не коррелируют. Для большего удобства введем полную корреляционную функцию
Н(г) = g(r) -1.
Влияние положения одной частицы на положение другой может происходить двумя способами: непосредственно или через положение некоторой третьей частицы. Можно показать [9], что при этом выполняется уравнение Орнштейна— Цернике:
h(rn) = c(rn) + рN |й(г1з)с(г2з)^гз, (!)
где c(r) — парная корреляционная функция, pN — численная плотность.
Для того чтобы замкнуть систему уравнений, необходимо еще одно уравнение, которое будет отражать некоторую корреляционную модель. В настоящей работе применялось ИУР HMSA. Оно представляет собой экспоненциальную интерполяцию между двумя аналитическими выражениями, построенными в приближении очень маленьких и очень больших расстояний [8]:
t(r) = h(r) - c(r), квТ
c(r) = -1 - t(r) + exp [-Рф0(г)] x 1 + exp[f (r)(t(r) -p<pi(r))]- 1
f (r)
r > r*
c(r) = y[t(r)] = -1 - t(r), r < r*, f (r) = 1 - exp(a r), Фо(г) = ф(г) - фО; Ф1(г) = ф(0, r < rm, Фо(г) = 0; Ф1(г) = ф(г), r > rm,
ф(0 = ШШ ф(г),
Г
ф(г < Г*) = «• = о,
\ йг )г=Г*
(2)
где ф0 — притягивающая часть потенциала, ф1 — отталкивающая часть потенциала, кв — постоянная Больцмана, гт — радиус минимума, г* — радиус максимума для внутренней границы обрезания потенциала, а — интерполяционный параметр, / (г) — интерполяционная функция.
Функция /(г) называется мостовой и обеспечивает однопараметрическую интерполяцию парной корреляционной функции. Параметр а неизвестен и подлежит определению из условия выполнения термодинамического равенства прямого и вириального выражения для изотермического модуля сжатия. Прямой изотермический модуль сжатия вычисляется непосредственным интегрированием радиальной функции распределения, а вириальный — дифференцированием фактора сжимаемости:
К* = дР
V* — Кт —
.дР$ Jт
да
1 + 4пр$ |[(г) - 1] г2йг
Z = в— = 1
Р N
ФРN |
К* - р' дР
й Ф(г)я(г )г Зйг,
)т V
(3)
Фх.(г)
N п /
^(г - гт)
п!
йпф(г ) йг
п=0 4 "" / г =гт
где N — порядок полинома. Применение такого потенциала не влияет на параметры модельной системы, поскольку суммарный потенциал оста-
ется равным заданному потенциалу парного взаимодействия частиц. Однако уравнение НМ8Л с таким разбиением потенциалов становится более устойчивым для численного решения, что позволило добиться повышенной точности результатов для ТТВ по сравнению с классической формой НМ8Л [2, 10].
Зная радиальную функцию распределения, можно найти значение избыточной внутренней энергии и давления:
и* = 2пвр$ |ф(г^(г)г2йг,
р _ %
в '
(4)
Чдр$л V ;т
где Z — фактор сжимаемости.
Также в работе использовалась улучшенная версия ИУР НМ8Л с применением сглаживающего потенциала ф^(г). Последний необходим для более корректного разбиения заданного потенциала парного взаимодействия частиц на притягивающую и отталкивающую части — ф0(г) и ф1(г) соответственно. Наиболее распространенной формой сглаживающего потенциала является полиномиальное разложение парного потенциала относительно точки минимума:
фо(г) = ф(г) - Фх(г); Ф1(г) = фх(г), г < гт, Фо(г) = 0; ф1(г) = ф(г), г > гт,
Таким образом, система (1)—(4) представляет собой уравнение состояния с учетом ближнего порядка для однокомпонентной смеси из сферически-симметричных частиц в неявном виде.
ВЫЧИСЛИТЕЛЬНАЯ ЧАСТЬ
В данной работе система уравнений решалась с использованием алгоритма, аналогичного алгоритму Лабика [11]. В этом методе применяется фурье-преобразование уравнения (1), что позволяет избавиться от трехмерного интеграла-свертки. Далее система (1)—(4) решается методом простой итерации с применением оптимизации шага по методу Ньютона—Рафсона только для первых
М точек сеточной функции 'Г(д1), представляющей собой дискретное фурье-отображение сеточной функции Т(г) = Ог^г. Таким образом, оптимизация высокого порядка сходимости проводится только для длинных гармоник функции ^г).
Подобный метод решения позволяет добиться компромисса между скоростью решения системы, устойчивостью решения и экономией оперативной памяти вычислительной техники. Как следствие, у этого алгоритма имеются врожденные недостатки — он не самый быстрый и теряет устойчивость при высоких плотностях для потенциалов взаимодействия с очень жесткой отталкивающей частью [12].
Для изучения ограничений устойчивости метода в работе было проведено собственное объемное исследование. Оно приводилось для потенциалов ЕХР-6 с различными параметрами жесткости а ф отталкивающей части. В ходе исследования было выяснено, что при больших значениях плотностей и больших числах М алгоритм теряет устойчивость на этапе оптимизации шага. При использовании малых чисел М или чистого метода прямой итерации (М = 0) алгоритм теряет сходимость. Однако если проводить вариативную подборку числа М, можно значительно расширить
0
0
14000 12000 10000 8000 6000 4000 2000
0 1 2 3 4 5 6 7 Рис. 2. Диапазон полученных значений на примере аргона. Ось абсцисс — температура, К; ось ординат — плотность, г/см3.
диапазон сходимости и, как следствие, диапазон применимости метода.
В качестве еще одного метода улучшения устойчивости была изучена возможность использования модернизированного замыкающег
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.