CC BY
6
УДК 531.19
Р. Н. Ласовский, ассистент (БГТУ);
Г. С. Бокун, кандидат физико-математических наук, доцент (БГТУ)
ПЕРЕХОДНЫЙ СЛОЙ НЕОДНОРОДНОЙ РАВНОВЕСНОЙ РЕШЕТОЧНОЙ СИСТЕМЫ В ДИАГРАММНОМ ПРИБЛИЖЕНИИ
Предложено комбинированное использование аналитических компьютерных преобразований и численных методов для расчета плотности свободной энергии системы и локальных значений химического потенциала неоднородной системы в диаграммном приближении. Основываясь на данном комбинировании, предложен усовершенствованный метод Ньютона для получения профиля распределения плотности в переходном слое для случая сильной неоднородности. Исследована температурная зависимость профиля плотности.
Computer analytical transformations are combined with computer numerical calculations for investigating the free energy density distribution and the local chemical potentials of inhomogeneous systems in the diagram approximation. Newton's advanced method for the description of the interfacial transition layer in the case of strong heterogeneity is offered. Temperature dependence of the particle density distribution in the transition layer is considered.
Введение. В работе [1] разработан устойчивый алгоритм интегрирования уравнения межфазного переходного слоя в равновесной решеточной системе в приближениях среднего поля и квазихимическом.
Однако известно, что оба эти приближения приводят к значительным погрешностям [2]. Поэтому ниже используем диаграммное приближение [2, 3], позволяющее воспроизводить термодинамические характеристики решеточных систем с высокой точностью и обобщенное на неоднородные системы в работе [4].
Для вычисления плотности свободной энергии системы и локальных значений химического потенциала неоднородной системы в диаграммном приближении приходится выполнять суммирование по узлам всей системы, т. е. фактически иметь дело с разностным представлением функционала свободной энергии и его производных [4].
В результате требуемые компьютерные ресурсы увеличиваются на порядки величины по сравнению с приближениями среднего поля и квазихимическим. Вместе с тем существует возможность комбинированного использования аналитических компьютерных преобразований и численных методов для системы решаемых уравнений большой размерности. Ниже показана реализация такой возможности.
Вычисление функционала свободной энергии и его производных в диаграммном приближении. Свободная энергия неоднородной системы представляется в форме [3]
ßF = ßFKX - D,
(1)
F„ -
где р = 1/ (квТ) - обратная температура; свободная энергия неоднородной системы в квазихимическом приближении; В - диаграммный вклад, равный
D
(2)
= C C C C
C = (W - 1)(PliPo,Pl jPoj ) 4 (Pl
1/2
П = *
Pli
П Po j )(Pi j +nfi Po,) ' ■Poj -W(Pij -Pii) ,
(3)
(4)
2Po
(Pli -Poj - W(Pl j -Pli )) , WPii
4p2, Poi
(5)
Здесь р17 (рш = 1 -р17) - плотность частиц (вакансий) в узле 7.
Конфигурация узлов 7, ], к, I при вычислении свободной энергии (1) по соотношениям (2)-(5) на плоской квадратной решетке представлено на рис. 1.
'3
,2
l
'1
k
j
'5
Рис. 1. Конфигурация узлов для расчета диаграммного вклада в характеристики системы
Далее будем рассматривать случай неоднородности с плоской границей раздела фаз. Тогда узлы, расположенные вдоль вертикальной оси в каждом столбце, одинаковы и в силу симметрии тождественны соответствующие мно-гоузельные функции (3) (ё. = di = ). Двух-узельная функция (4) симметрична по переменным р17 и р1 ^. Отмеченную симметрию учтем
Механика
25
при расчете производных функционала свободной энергии (1).
Дифференцируя свободную энергию по плотности числа частиц ри, получим локальную величину химического потенциала, которую в соответствии с уравнением (1) можно записать в виде
h = цГ + дц<:
(6)
где - значение химического потенциала в квазихимическом приближении; Ац.- - диаграммная добавка [3], равная
Дц,. = -А,-
ds
s=1 dp1. ■
(7)
В свою очередь, из выражения (4) и рис. 1 следует, что переменная ри содержится в функциях С р, С. , С ^ и Са. Для унификации записи последующих выражений введем обозначение
DpCv =
1 dCj
1 C
j Cj dp,. C. фи •
(8)
Из изложенного следует, что в сумме (7) остаются только те слагаемые, которые содержат множитель Су. Как видно из рис. 1, это будут
слагаемые а., а., а., а., так что
. 1 h Ъ
дц , = а
' ФиV
,)• (9)
С учетом обозначений (8) выражение (9) перепишем в следующем виде:
Ац . = -Х(ё г (БрСу + БрС а) + ^ (БрСу + БрС 1) +
+ ёчфрСгч + БрСщ) + ^(ЯрС. + БрСй)). (10) Принимая во внимание, что
ё г = ё^ £>рСй = DpCI.i, ё. = ё^
находим
(11)
Дц,=-2Ц d,. (DpCj + DpCa )-
М DpC , .
DpCa)
При рассматриваемой неоднородности вдоль горизонтальной оси выражение (11) далее не упрощается. Соотношение (11) определяет диаграммную добавку к химическому потенциалу неоднородной системы с осью аксиальной симметрии. Полученное выражение удобно тем, что численные процедуры расчета производных можно заменить аналитическими формулами, что существенно сократит объем
компьютерных ресурсов и повышает точность конечных результатов.
Алгоритм численного интегрирования.
Изложенный в работе [1] метод определения профиля распределения плотности применим и для решения задачи для случая сильной неоднородности, когда разложение в форме [1]
Ц .- = £ (Ри) + С (Р11+1 +Р11-1 - 2Ри) = Ц С> 0 (12)
(13)
становятся неприменимыми. При нарушении условия
ц,= const = ц = ц (pG) = h (pL), (14)
где Д - значение химического потенциала, найденное по формуле расчета ц в однородной фазе; pG, pL - плотности соответственно в однородной разреженной и плотной сосуществующих фазах, химические потенциалы в узлах оказываются не постоянными, что с физической точки зрения обусловливает потоки, приводящие к изменению плотности числа частиц по всему объему системы. В этом случае задается пробный профиль и вводятся невязки решения
Д,+1 ^(a-MMX
(15)
а через них определяются вариации плотности в слое
Sp,^ ,,+1 +д «-1 =
(16)
так что в процессе решения плотность р.. варьируется на величину 5р { и выполняется следующий шаг процедуры интегрирования.
Параметр О выбирается так, чтобы обеспечить сходимость решения задачи и его точность, определяемые величиной 5р.. на последнем шаге расчетов.
Предложенный метод является по существу усовершенствованным методом Ньютона. В стандартном методе Ньютона р.. рассматриваются как независимые переменные и при поиске решения оперируют величинами 5ц i / фу. В предложенном подходе независимые переменные р { дополнительно ранжируются в некотором фиктивном х-пространстве. Из названных переменных конструируется функция р( х), вообще говоря, произвольным образом. Далее в новом построенном фиктивном пространстве с помощью выражений (6) и (9) определяем химические потенциалы, через которые вводим невязки (15), и с помощью (16) обеспечиваем процесс, который в силу своей
и
необратимости с неизбежностью приведет к решению системы сколь угодно большого числа уравнений с заданными граничными условиями.
Компьютерные преобразования и результаты вычислений. Для выполнения компьютерных аналитических преобразований вводим оператор-функции, соответствующие выражениям (3)-(5):
n(x1, xlj,W) = -
xl - (1 - xl j) - W(xl j - xl)
2(1 - xl)
г
\
(xl - (l - xl j) - W (xl j - xli)) Wx1i
, (l7)
C (xl, xl j, W) =
4(l - xli)2 l - xli
(W -1)) xli(l - xli) xl j (l - xl j)
xli + (l - xl j )n(xli, xl j, W)
(l8)
x1 j + (1 - x1i)n(x1i, x1 j, W)'
dl (x1, x1 j, W) = C (x1, x1 j, W)2 x
xC(x1, x1, W)C(x1 j, x1 j, W). (19)
Далее производную (8) определяем идентификатором
dC (xl, xl j ,W)
DpC (xl, xl j ,W ) = -
dxli
C (xl, xl j ,W) '
(20)
что позволяет рассчитать диаграммную добавку (9) как
Ац(х1, х1 у, х. 3, Ж) = -2А,((х1, х1 у, Ж)х
х (£>рС (х1, х1 у, Ж) + £>рС (х1, х1, Ж)) +
+ (3, х1, Ж)(£>рС(х1, х. 3, Ж) +
+£>рС(х1, х1, Ж))). (21)
В свою очередь квазихимическую часть химического потенциала и химический потенциал определяем функциями
цкх(х1, х.3, х1 у, Ж) = 1п 1 . - 1пп((1, х1 у, Ж) -
- 21пп(х1, х1, Ж )- 1пп(х1, х. 3, Ж ) , (22) ц(х1, х.3, х1 у, Ж) = цкх(х1, х.3, х1 у, Ж) +
+ Ац(х1, х.3, х1 у,Ж). (23)
Соответственно, вариацию плотности (16) запишем в виде
5р. = °(ц(р,-1, р.-2,р., Ж)- , р.-1, р<+1, Ж) + + Ц(рг+l, р. , Р.+2, Ж)), (24)
которая позволяет наити плотность частиц в слое i:
рг=р+5рг.
(25)
Зависимость концентрации числа частиц от координаты для неоднородной системы при различных температурах в единицах критической, полученная описанным выше методом, показана на рис. 2.
Р l,0
0,6 0,4 0,2 0,0
0
10 20 30 40 Номер слоя
Рис. 2. Равновесное распределение плотности частиц:
1 - t = 0,8; 2 - t = 0,9; 3 -1 = 0,95
Заключение. Как и следовало ожидать, ширина переходного слоя при повышении температуры увеличивается от 10 до 20 постоянных решетки. Таким образом, рост температуры приводит к размытию переходной области, и при температурах выше критической переходная область вовсе отсутствует.
Литература
1. Ласовский, Р. Н. Определение характеристик переходного слоя в неоднородной равновесной решеточной системе / Р. Н. Ласовский, В. С. Вих-ренко // Труды БГТУ. - 2011. - № 6: Физ.-мат. науки и информатика. - С. 21-23.
2. Вихренко, В. С. Равновесные и диффузионные характеристики интеркаляционных систем на основе решеточных моделей: монография / В. С. Вихренко, Я. Г. Грода, Г. С. Бокун. -Минск: БГТУ, 2008. - 326 с.
3. Vikhrenko, V. S. The diagram approximation for lattice systems / V. S. Vikhrenko, Ya. G. Groda, G. S. Bokun // Phys. Let. A. - 2001. -Vol. 286, № 2/3. - P. 127-133.
4. Ласовский, Р. Н. Диаграммное приближение для неравновесных и неоднородных состояний решеточных систем / Р. Н. Ласовский, Г. С. Бокун // Труды БГТУ. Сер. VI, Физ.-мат. науки и информатика. - 2010. - Вып. XVIII. -С. 56-58.
Поступила 02.03.2011
