Algorithm for finding parameters of equations of state based on the particle swarm optimization

Толық мәтін

Введение

Исследование посвящено разработке алгоритма для поиска параметров (констант) уравнений состояния вещества. Константы рассматриваемых моделей определяются с помощью численной процедуры на основе метода роя частиц [1, 2]. Это стохастический алгоритм, в котором группа частиц перемещается в пространстве решений и осуществляет поиск оптимума. Весь рой после смещения обменивается данными о наилучшем положении, в котором каждому отдельному элементу (частице) удалось побывать к настоящему моменту (итерации). Следующий шаг каждой частицы делается с учетом полученной информации. Поиск оптимального положения заканчивается тогда, когда весь рой смещается к некоторой наилучшей позиции с требуемой точностью. Скорость сходимости зависит от значений параметров и используемой модификации метода [3, 4]. В этой работе параметры численного алгоритма подобраны для поиска констант уравнения состояния вольфрама в жидкой фазе при высоких давлениях и температурах.

Описание термодинамических свойств исследуемого металла осуществляется в рамках модели на основе уравнения Ван-дер-Ваальса [5, 6], уравнения Ликальтера [7], которое получено при рассмотрении модельной задачи заряженных твердых сфер, и уравнения Каплуна–Мешалкина [8], модифицированного Куперштохом [9]. Эти модели содержат малое количество свободных констант (не более четырех), что обеспечивает максимальную простоту разрабатываемых на их основе уравнений состояния.

Ранее теоретическое описание термодинамических свойств вольфрама при высоких давлениях и температурах проводилось в рамках более сложных моделей [10–18], а также в рамках метода квантовой молекулярной динамики [19]. Помимо этого, жидкий вольфрам исследовался экспериментально, в частности, методом импульсного нагрева [20, 21]. Получены данные по изобарическому расширению жидкой фазы металла при давлениях 0.1 МПа [22], 0.1 [23], 0.2 [24] и 0.3 ГПа [25]. Состояния жидкой фазы вольфрама при более высоких давлениях исследованы в ударно-волновых экспериментах [26–31]. С помощью данных по ударному сжатию [26–30] и изоэнтропическому расширению [31] ударно-сжатых исходно пористых образцов строится калорическое уравнение состояния и определяется область применимости рассматриваемых моделей.

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

Модели уравнений состояния

Как уже было отмечено выше, в работе рассмотрены три модели, первая из которых (M1) определяется уравнением Ван-дер-Ваальса [5, 6], задающим следующую зависимость давления (P) от удельного объема (V) и температуры (T):

$P\left(V,T\right)=\frac{RT}{V-b}-\frac{a}{V^2}$, (1)

где R – удельная газовая постоянная (отношение универсальной газовой постоянной к молярной массе рассматриваемого вещества); a и b – константы.

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

$P\left(V,T\right)=\frac{RT}{V-b}-\frac{a}{V^{4/3}}$. (2)

Третья модельная зависимость (M3) давления от удельного объема и температуры предложена в [8] и впоследствии модифицирована в [9]. Она может быть представлена в виде

$P\left(V,T\right)=\frac{RT}{V-b}\left(1+\frac{c-b}{V}\right)-\frac{a}{V^2}$, (3)

где c – третий параметр модели. Было показано [8], что выражение, определяющее фактор сжимаемости, должно содержать слагаемое, описывающее поведение вещества при предельно высоких температурах и зависящее только от удельного объема. Также установлена [8] структура трехпараметрического термического уравнения состояния (3) для сверхкритического флюида, жидкости и газа, позволяющего более точно количественно описать свойства ряда веществ, чем некоторые другие известные простые модели.

Для построения термодинамически полного уравнения состояния на основе выражения P = P(V, T) необходимо сформулировать вид зависимости удельной внутренней энергии (E) от удельного объема и температуры. Для уравнений (1), (2) или (3) такая функциональная связь предлагается в едином виде

$E\left(V,T\right)=\chi RT+f\left(V\right)$, (4)

где χ – безразмерный параметр; f(V) = =∫ [T(∂P/∂T)_V − P] dV. Нетрудно заметить, что тогда удельная теплоемкость при постоянном объеме является константой: c_V = (∂E / ∂T)_V = χR.

Параметры термического уравнения состояния в каждой модели определялись из условия оптимального описания данных по изобарическому расширению жидкой фазы [22] с помощью метода роя частиц. Предпочтение эксперименту [22], проведенному при атмосферном давлении в воде, было отдано в связи с результатами расчетов [13] и выводами, сделанными о данных по изобарным измерениям [22–25] в [17], где сформулировано предположение о возможном отклонении от условий равновесия в опытах [23–25], проведенных при повышенном давлении в газовой среде (аргон). Значение параметра χ, входящего только в выражение для внутренней энергии (4), находилось из условия наилучшего описания данных по ударной сжимаемости исходно пористых образцов вольфрама [27–30]. Область применимости полученных уравнений состояния определялось на основе информации, полученной в экспериментах с волнами ударного сжатия [27–30] и изоэнтропического расширения [31].

Метод роя частиц

В последнее время предлагаются различные способы обработки и анализа экспериментальной информации [32, 33]. В настоящей работе предлагается способ оптимизации описания данных по изобарическому расширению жидкого металла с помощью метода роя частиц, представляющего собой стохастический алгоритм, в основе которого лежит использование набора элементов (частиц). Каждая частица является вектором в некотором многомерном пространстве решений, координаты которого в начальный момент задаются случайным образом. Крайне желательно, однако, чтобы частицы исходно были равномерно распределены по исследуемой области пространства [1]. Это можно обеспечить, например, путем задания достаточно большого количества элементов роя. Есть и другие способы, связанные с некоторой упорядоченностью задания начального положения [34]. В настоящей работе исходные координаты каждой частицы определяются полностью случайным образом.

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

${\delta }_{\rho }=\frac{1}{\sqrt{N}}\times \sqrt{\sum _{n=1}^N{\left[\left(\frac{{\rho }_{k-1}}{{\rho }_n}\right)\frac{T_k-T_n}{T_k-T_{k-1}}+\left(\frac{{\rho }_k}{{\rho }_n}\right)\frac{T_n-T_{k-1}}{T_k-T_{k-1}}-1\right]}^2}.$ (5)

Здесь ρ – плотность, ρ = V⁻¹; индекс n соответствует значениям температуры и плотности на изобаре в экспериментальных точках; N – количество экспериментальных точек; индексы k и k−1 (такие, что T_k−1 ≤ T_n ≤ T_k) относятся к результатам расчетов по уравнению состояния с набором параметров, соответствующим координате рассматриваемой частицы роя. Количество N_K температурных интервалов (T_k−1;T_k) выбирается таким, чтобы результат расчета (5) практически не зависел от N_K.

Далее происходит смещение частиц роя, которое реализуется путем задания скорости каждого элемента на следующем временном шаге (на следующей итерации):

$v_{\mathrm{i}}^{\mathrm{t}+1}=w{v}_{\mathrm{i}}^{\mathrm{t}}+c_1{r}_1\left(p_{\mathrm{i}}^{\mathrm{t}}-x_{\mathrm{i}}^{\mathrm{t}}\right)+c_2{r}_2\left(g^{\mathrm{t}}-x_{\mathrm{i}}^{\mathrm{t}}\right).$ (6)

Здесь нижний индекс определяет номер частицы, а верхний – номер итерации. Переменные x и v с верхними и нижними индексами указывают положение и скорость частицы с номером i в момент, соответствующий верхнему индексу (под моментом времени подразумевается номер итерации); p_i^t – наилучшее положение (в котором функционал принимает наименьшее значение), достигнутое к итерации t частицей с номером i; g^t – глобальное наилучшее положение среди всех элементов роя за все прошедшие итерации до t; w, c₁ и c₂ – параметры метода роя частиц; r₁ и r₂ – случайные числа из промежутка (0; 1) для каждой координаты векторов (p_i – x_i) и (g – x_i), причем r₁ и r₂ различны для каждой координаты каждой частицы на каждой итерации. В этой работе полагается v_i⁰ = 0.

После задания скорости по формуле (6) осуществляется смещение роя:

$x_{\mathrm{i}}^{\mathrm{t}+1}=x_{\mathrm{i}}^{\mathrm{t}}+v_{\mathrm{i}}^{\mathrm{t}+1}$. (7)

Затем проверяется условие выхода из итерационной процедуры, которое сводится либо к сходимости метода, либо к достижению максимального числа итераций. Сходимость метода определяется по текущему положению всех частиц: если они локализуются в достаточно малой области (здесь малой областью считается окрестность глобального минимума с отклонением от его положения не более чем на тысячную долю процента), то делается вывод, что оптимум достигнут. Принцип работы алгоритма проиллюстрирован на рис. 1.

 

Рис. 1. Иллюстрация работы алгоритма на основе метода роя частиц. По осям отложены параметры уравнения M1. Черные кружочки показывают положение частиц при t = 0, красный квадратик – при t = 728 (частицы сошлись в одну точку с требуемой точностью). Зеленый крестик определяет положение глобального минимума функционала при t = 728. Проиллюстрирован тот факт, что за время работы алгоритма оптимум достигнут.

 

Для качественной работы метода необходимо удачно задать параметры w, c₁ и c₂. Коэффициенты выполняют следующие функции: w (инерционный вес) позволяет сгладить изменение скорости на следующем временном шаге, он может зависеть от времени или быть постоянным; c₁ обеспечивает смещение частицы к собственному наилучшему положению; c₂ перемещает элемент роя к глобальному лучшему положению. Для каждой задачи константы c₁ и c₂ уникальны. В данном случае инерционный вес считается постоянным в течение всего времени работы алгоритма, но задается случайным образом в диапазоне 0.5 < w < 1. Остальным константам присваиваются значения с₁ = 1.5, с₂ = 2.6.

Помимо указанных параметров метода, важную роль играет количество частиц в рое (N_I). На рис. 2 показано количество итераций, за которое достигается оптимум, и отклонение от экспериментальных данных для фиксированного числа векторов, участвующих в расчетах. На этом рисунке видно, что, начиная с определенного числа частиц, оптимум может быть достигнут. В настоящей работе для моделей M1 и M2 было задействовано 200 частиц, а для модели M3 – 1600.

 

Рис. 2. Число итераций N_T (за время работы алгоритма) и достигаемое наименьшее отклонения от экспериментальных данных δ_ρ (5) в зависимости от числа частиц N_I. По оси абсцисс отложено число частиц N_I.

 

Результаты

Параметры уравнений состояния M1, M2 и M3, найденные по представленному алгоритму, приведены в табл. 1. Можно заметить, что значение коэффициента b, имеющего физический смысл исключенного объема, принимает близкие значения для всех рассмотренных моделей.

 

Таблица 1. Параметры уравнений состояния для жидкой фазы вольфрама по моделям M1, M2 и M3

Модель	a, Дж∙см3∙г−ξ	ξ	b, см3∙г−1	c, см3∙г−1	χ
M1	84.19876	2	0.05138	–	3.465134
M2	481.2655	4/3	0.05017	–	3.377452
M3	40.42995	2	0.05322	0.01569	3.427462

 

Как уже было отмечено во введении, границы области применимости полученных уравнений состояния устанавливаются путем сравнения результатов расчетов по ним с данными ударно-волновых экспериментов. На рис. 3 расчетные ударные адиабаты и изоэнтропы разгрузки ударно-сжатых исходно пористых образцов исследуемого металла показаны в сопоставлении с данными работ [27–31]. Продемонстрировано хорошее согласие результатов расчетов, проведенных по моделям M2 и M3, в более широком диапазоне давлений, чем по модели M1, с данными по разгрузке [31].

 

Рис. 3. Ударные адиабаты (H1–H6) и изоэнтропы разгрузки (S1 и S2) для образцов вольфрама с исходной плотностью ρ₀₀ = 4.60 (H1), 5.50, 6.64, 8.87 (S1, S2), 10.59 и 13.36 г∙см⁻³ ρ₀₀ = 13.36 (H1), 10.59, 8.87 (S1, S2), 6.64, 5.50, и 4.60 г∙см⁻³ (H6); маркеры – экспериментальные данные (P0, P3, P5, P6, P9, K3 – [29]; P1, P5 – [30]; P2, P8, K2 – [28]; P7, K1 – [27]; R1, R2 – [31]), из которых часть (K1, K2, K3) выходят за пределы области V > b. Изоэнтропы S1 и S2 соответствуют массовой скорости на ударной адиабате U = 2.7 и 3.11 км∙с⁻¹; лиловые, фиолетовые и голубые линии – результаты расчетов изоэнтроп по моделям M1, M2 и M3 соответственно с учетом образования двухфазной смеси жидкость–пар (сплошные линии) и метастабильных однофазных состояний (штриховые линии).

 

На основании сравнения расчетных давлений и массовых скоростей на ударных адиабатах и изоэнтропах по рассматриваемым уравнениям состояния с экспериментальной информацией можно заметить, что модель M1 применима при давлениях от 100 до 5 ГПа, а M2 и M3 – до 0.03 ГПа.

Заключение

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

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

Авторлар туралы

K. Boyarskikh

Joint Institute for High Temperatures of the Russian Academy of Sciences; Moscow Institute of Physics and Technology

Хат алмасуға жауапты Автор.
Email: boyarskikh.ka@phystech.edu
Ресей, Moscow; Dolgoprudny

K. Khishchenko

Moscow Institute of Physics and Technology; South Ural State University; South Ural State University

Email: boyarskikh.ka@phystech.edu
Ресей, Moscow; Dolgoprudny; Chelyabinsk

Әдебиет тізімі

Қосымша файлдар