CC BY
43
А. В. Клинов
ЧИСЛЕННЫЙ МЕТОД РЕШЕНИЯ УРАВНЕНИЯ ОРНШТЕЙНА-ЦЕРНИКЕ ДЛЯ МНОГОКОМПОНЕНТНЫХ СИСТЕМ С ДЛИННОДЕЙСТВУЮЩИМИ КУЛОНОВСКИМИ МЕЖМОЛЕКУЛЯРНЫМИ ВЗАИМОДЕЙСТВИЯМИ
Ключевые слова: уравнение Орнштейна-Цернике, радиальная функция распределения, кулоновские взаимодействия.
Предлагается вариант переформулирования уравнения Орнштейна-Цернике для многокомпонентных систем с межмолекулярным дальнодействующим кулоновским взаимодействием. Для расчета радиальных функций распределения и их производных по термодинамическим параметрам адаптирован эффективный алгоритм Лабика-Малиевского.
Keywords: Ornstein-Zernike equation, radial distribution function, the Coulomb interaction.
A version of the reformulation of the Ornstein-Zernike equation for multicomponent systems with intermolecular long-range Coulomb interaction was suggested. For the calculation of the radial distribution functions and their derivatives with respect to the thermodynamic parameters adapted efficient Labik-Malijevsky algorithm.
Теория интегральных уравнений для частичных функций распределения является одним из перспективных направлений современной молекулярно-статистической теории [1]. В ее основе лежит факт, что многие термодинамические свойства газов и жидкостей определяются локальным молекулярным окружением,
возникающим вокруг одной частицы в радиусе порядка 10 А0. Достоверно описать такое локальное окружение в большинстве случаев оказывается возможным на основе функций распределения, зависящих от координат небольшого числа молекул, например радиальной функции распределения (РФР). В этом случае теория интегральных уравнений для частичных функций оказывается значительно эффективней в вычислительном плане по сравнению с методами молекулярной динамики или Монте-Карло, и в тоже время фундаментально строгой. Основным уравнением этой теории является интегральное уравнение Орнштейна-Цернике (ОЦ), которое в случае т-компонентной смеси представляет следующую систему интегральных уравнений:
(1)
т
Ь ар (Г) = Сар (Г) + 2 Р* IС(г)Ь Р* (г)а^
Я.=1 у
а = 1... т; р = 1... т
где И(г) - полная корреляционная функция; с(г) -прямая корреляционная функция; р-числовая
плотность; нижние греческие индексы означают сорт молекулы. Функции, входящие в уравнение ОЦ, имеют следующую связь с РФР g и потенциалом межмолекулярного взаимодействия ф:
(
h
aß (r) = gaß (Г) -1 = eXp
Фаъ(г)
квТ
Л
®aß (Г)
-1 (2)
Cxß (Г) = haß (г) -юар (г) + Baß (г) (3)
где kB - константа Больцмана; Т - температура; ю(г) - термический потенциал; В(г) - бридж-функционал, который представляет собой бесконечный ряд неприводимых диаграмм.
В работе рассматриваются системы с центральными межмолекулярным взаимодействием вида:
фар(г)=фар(г)+Фар(г) (4)
здесь фя - короткодействующий потенциал типа Леннард-Джонса (ЛД):
Фaß (r) = 4saß (aß / г)12 - (aß /r)
(5)
или какой-либо другой (Букингема, твердой сферы и т.д.); ф1 - длиннодействующая часть, описываемая законом Кулона:
Ф aß(Г) =
Za Zß
(6)
В выражениях межмолекулярного
(5), (6) для потенциалов
взаимодействия ст -эффективный молекулярный диаметр, определяется из условия ф8(ст)=0, е -глубина потенциальной ямы, ъ - заряд.
Несмотря на то, что молекулы реальных веществ в основном не являются сферическими, рассматриваемые в работе системы с потенциалом типа (4), представляют большой научный интерес в плане изучения возможностей теории интегральных уравнений, а так же практический интерес при использовании их в качестве грубой модели для описания термодинамического поведения широко класса веществ [2], [3], в том числе ионных жидкостей.
Особенностью рассматриваемой системы, является наличие прямого электростатического взаимодействия (6) между частицами.
Термодинамически устойчивыми такие системы могут быть только при некоторых условиях, например: если кинетическая энергия молекул
намного больше потенциальной или если плотность частиц достаточна для того что бы экранировать длиннодействующее взаимодействие между двумя молекулами за счет коллективных эффектов. Для возможности численного решения уравнения ОЦ в случае систем с кулоновским взаимодействием,
проведем его переформулировку. Учитывая, что для термодинамически устойчивого состояния полная РФР должна быть короткодействующей, запишем:
(
Ь(г) = ехр
-ф '(г) -ф ‘(г)
Л
= ехр
V
^-ф5(г) квт
квт
квТ
+ ю(г)
-1 =
(7)
+ ю8(г)
V
-1
Соответственно короткодействующая прямая корреляционная функция будет иметь вид:
С (г) = с(г) + ф1 (г) = Ь(г) - ю8 (г) + Б(г) (8)
Запишем систему уравнений ОЦ в Фурье пространстве в матричном виде:
Ь = с + Ь рс (9)
здесь и далее жирным обозначаются матрицы, волной Фурье образ, р-диагональная матрица плотностей компонентов. Выполнив замену (7), (8) получим:
Ь (I+ Р91 ) = с8-ф1 +11 рс8 (10)
Введем следующие определения:
С = <ф1 (I + р(ф1 )1 (11)
(1) = I -ря
(12)
В итоге получим следующую систему переформулированных уравнений ОЦ:
~ = г 8<3(1) + ьр?8 <3(1) - ~ (13)
Здесь нужно отметить, что решением системы уравнений (11), является эффективный потенциал Дебая-Хюккеля:
Яар(г) = ЪаЪр5ХР<-ХЕ1) (14)
где Хо =-
4л
рЛ - обратный радиус Дебая.
квТ' л
Численное решения уравнения (13) осуществлялось с использованием адаптированного гибридного метода Лабика-Малиевсого [4], в котором итерации строятся по короткодействующей непрямой корреляционной функции:
Г(г) = гу8 (г) = г(Ь(г) - с8 (г) + д(г)) .
Перепишем уравнение (13) в виде:
у8 = р-1(з(1)рс (3(1) (I - рс83(1)) - с8 (15)
Используя очевидное соотношение: р-10 (1)р = I - С р = (3(2), получим:
у8 = (3(2)с8(3(1) (I - рс8<3(1)) - с8 (16)
В таком виде система уравнений ОЦ для кулоновских систем может численно решаться во всем диапазоне концентраций и плотностей, включая их нулевые значения, в отличие от традиционно используемой в литературе записи [5], в которой присутствует операция деления на плотность.
Введем следующие обозначения: I"(1) = 8(1) ,
С8 (1) = й8 (1) и окончательно запишем:
і" = (3 (2)с 83(1)
с
№
- с8
(17)
вычислительного
I - с 8<з(
ч. Ї У
Основные этапы
алгоритма:
1) Задаются термодинамические условия: температура, плотность, состав, а так же потенциалы межмолекулярного взаимодействия.
2) Пространственный интервал, на котором ищется РФР, разбивается на дискретный набор точек і = 1,2,...,К с шагом Дг между ними.
Причем N = 2е, что необходимо для осуществления численной процедуры «быстрого» Фурье-преобразования. Обычно, достаточно
пространственного интервала порядка 10 - 20с , а удовлетворительную точность обеспечивают
Є = 8 -12 и — = 0.025 - 0.01 . На выбранном с
дискретном наборе точек задается начальное приближение значений Г0.
3) Определяется значения С^р = р(гар) в
соответствии с выбранным замыканием уравнения ОЦ и совершается дискретное Фурье-преобразование для функций Г и С8:
N -1
Г = 4 лДгX Гі 8Ш(У л / N , ] = 1,2,... ^ -1
і=1
Координата в обычном пространстве определена как гі = іДг , в Фурье-пространстве - как ^ = _)Д1 ,
между ними должно выполняться следующее
. . л соотношение ДгД1 = — .
N
4) Выбирается некоторое число точек і = 1,2,..., М , на которых решение будет искаться методом Ньютона-Рафсона. Обычно бывает достаточно задать такое М, при котором МДг« 0.5 - 0.7с . Решается система
линеаризованных уравнений (17) относительно ДГ):
X ^аР,і;кЛ.,РГ^ = 8ард (18)
к^,І
где
8арі =^(3(2)с (3(1)|1 - ^ С (3(1) |- с1
а,і
ар,і
Для определения Якобиана в (18) необходимо иметь аналитический вид замыкания в Фурье
пространстве: Сар = Р(Гар) что не представляется
возможным. Поэтому его определяют приближенно
разложив функцию С® в точке Г0 и оставив
линейный член:
•'ар.
і(Гар,і) = Сар,і(Г ар,і) +(с1^аГаР,і )г
ДГ
1 ар,і _і ар,і
Фурье-преобразование этого соотношения
ар,і
даст:
N -1
С8 = С8 (Г0 ) + V С(2) (Г
ар, ар,) ар,) / ^ ар,) ар,п
-Г'
ар,п
).
Коэффициенты Cjjn определяются выражением
'S]. = Ц (d^dГaß,)Г^,=Г:,.l sin ^ ) si" ^
Не сложно заметить, что
N-11 . SC
Z C(2) ( - Г0 ) -
aß,jn y aß,n aß,n J
s
—ЛГ .
Sj^ aß,j
n—1 '-/x aß,j
Таком образом, Якобиан определяется согласно выражения:
Г
-J ßi . . = -^-s SßAS.. =
aß,.;^,. (к ßA ij
1 кЛ,.
: Г A(1) A(2)-s SßJ С2^.-S SßAS..
|_ ак Aß ак ßA J кЛ;у ак ßA ij
где
A(1) = Q(1) (i - p С sQ(1)
A(2) = Q(2) 11 - СsQ(1) p
пока
V
Система уравнений (18) решается до тех пор
>^1
Í M Л
ілг2
V j=1
корректируя на каждом
шаге Г- ГJ = Г° + AГJ.
Значение ^1 «10-7 . Если условие по точности выполняется, переходят к следующему пункту.
5) На интервале М N вычисляется
следующее приближение для Гj по (17):.
6) Совершают обратное Фурье-
преобразование для Гj:
А+ N-1
Г= —Г ЕГ;8ш(ул /N),
2ж2
i = 1,2,...,N -1,
j=1
проверяют условие
І
Z(Гі-Г0) I > ^
проводят корректировку Г¡ = Г. . Если условие точности не выполняется (обычно Е,2 да 10 З ), то
возвращаются к пункту З, в противном случае вычисления заканчивают.
Представленный выше алгоритм тестировался для двухкомпонентных систем с потенциалом межмолекулярного взаимодействия ЛД (З) + (6), для различных значений зарядов, величины которых
* z
задавались в приведенных единицах z = ,— . В
Vas
работе использовалось гиперцепное замыкание уравнения ОЦ B(r) = 0 . Расчеты показали, что для большей части фазовой диаграммы для систем с z < 3 , решение сходиться нулевого
приближения Г0 = 0 . Как правило, для получения решения хватает порядка семи прямых итераций, а
внутри одной прямой итерации совершается в среднем пять «ньютоновских» итераций.
* о
В случае ъ > 3 , происходит значительная
перестройка РФР по сравнению с функцией
exp
í-фЧг)^
kBT
(рис.1), которая используется в
V “Б- /
качестве начального приближения, поэтому рекомендуется выполнить сначала несколько десятков прямых итераций, а затем включать «ньютоновские». Результаты расчетов РФР для
кБТ * 3
T =
=З
p = pa = 0.9 и различных
значений зарядов представлены на рис. 1, 2.
Данные по РФР позволяют рассчитать давление и внутреннюю энергию [6]. Другие термодинамические величины можно определить на основе дифференциальных соотношений Максвелла. Наиболее просто и точно определить производные от давления и внутренней энергии по различным термодинамическим параметрам, можно на основе соответствующих производных РФР [7]. Эти производные находятся из решения продифференцированного уравнение ОЦ по требуемым термодинамическим параметрам. Продифференцируем уравнение (17) по некоторому параметру п:
Ґ Л -ч^С)
SSL
Sn ^
где
M(1) =
M(1) + M(2) — Q(1)
Sn
M(З) -
ge0
Sn
(19)
SQ
(2)
(К Sp 1
Sn
- + (Г + С °)-^-Sn t
С^(1) +
+м(2)C0
SQ
(1)
Sn
m (2) = Q(2) + (г+С0)p
M(3) =| I + ■£ С^
SP
Ели n - плотность, то — - матрица мольных долей, Sn
если n-температура
Sp
Sn
= 0 и т.д. Решение
уравнения (19) осуществляется после решения уравнения (17), используя приведенный выше алгоритм, где Якобиан определяется как:
аг „
aß,і
Sn
aß^^AJ
Sr
'-SaicSßASiJ =
кЛ, j
Sn
m (2) [Q(1) m (3)
Aß SaicSßA
С к2);і] - S (к S ßA S ij
~ (2)
С - записывается через соответствующие производные с учетом используемого замыкания.
Рис. 1 - РФР между разноименно заряженными молекулами; сплошная линя - z*=0, короткий пунктир- z*=3, длинный пунктир z*=10
Рис. 2 - РФР между одноименно заряженными молекулами; сплошная линя - z*=0, короткий пунктир- z*=3, длинный пунктир z*=10
Работа выполнена при поддержке РФФИ, грант10-08-01024-а
Литература
1. Мартынов Г.А. Проблема фазовых переходов в статистической механике/ Г.А. Мартынов // Успехи физических наук.- 1999.- Т.-169.- №6
2. Клинов А.В. Термодинамические свойства веществ на основе сферически-симметричного потенциала межмолекулярного взаимодействия / А.В. Клинов , А.В Малыгин, Л.Р. Минибаева, Г.С Дьяконов. // Вестник Казанского технологического университета.- 2010. - №1. - С.105-109
3. Клинов А.В. Использование модели центральных сил для описания двухцетровых молекул / А.В. Клинов , И.П. Анашкин // Вестник Казанского технологического университета.- 2011. - №3. - С.272-274
4. Labik, S. A rapidly convergent method of solving the OZ equation / S. Labik, A. Malijevsky, P. Vonka // Mol. Phys. -1985. - V.56. - № 3. - P. 709 - 715.
5. Toshiko I. Accurate integral equation theory for the central force model of liquid water and ionic solutions / Haymet A.D.J. // J.Chem.Phys. - V89.- №1. - P. -4315
6. Балеску, Р. Равновесная и неравновесная статистическая механика / Р. Балеску: пер. с англ. - М.: Мир, 1978. - 808 с
7. Вомпе, А. Г. Проблема термодинамической согласованности решений уравнения Орнштейна-Цернике / А. Г. Вомпе, Г. А. Мартынов // Журн. физ. химии. - 1994. - Т. 68. - № 3. - С. 41-46.
© А. В. Клинов - д-р техн. наук, проф., зав. каф. процессов и аппаратов химической технологии КНИТУ, а1к1т@к81и.гц.
