янии на них нелинейности уравнения ПБ. Раздел VI содержит краткое заключение.
II. ОПИСАНИЕ ЗАДАЧИ И ОСНОВНЫЕ УРАВНЕНИЯ
Рассматриваемый в настоящей работе коллоидный кристалл показан на рис. 1. Он состоит из бесконечно длинных цилиндрических коллоидных частиц радиуса Я, образующих двумерную гексагональную кристаллическую решетку с постоянной а. Векторы а1 и а2 - векторы примитивных трансляций решетки. Система частиц погружена в симметричный бинарный одновалентный электролит. Частицы являются абсолютно твердыми диэлектрическими стержнями. Предполагается, что диэлектрическая проницаемость материала частиц гораздо меньше, чем проницаемость электролита, поэтому во всех последующих вычислениях она принята равной нулю. Частицы характеризуются однородной поверхностной плотностью заряда а, которая во всех обстоятельствах остается постоянной (так называемая модель постоянного заряда).
Мы полагаем, что рассматриваемая система подчиняется общему нелинейному уравнению ПБ [3, 4]
V2 Ф = -1
q zrnr
££n
-exp
qeZ-ф " kT
(1)
для электрического потенциала ф в области электролита вне частиц и уравнению Лапласа
V2 ф = О
(2)
О О О О О
3 2
3
внутри частиц, причем для записи уравнений используется рационализированная система единиц СИ. Здесь п1 - равновесная концентрация /-го компонента в объеме электролита, - валентность /-го компонента, де - элементарный заряд, £ - относительная диэлектрическая проницаемость электролита, £0 - электрическая постоянная, к - постоянная Больцмана, Т - абсолютная температура.
Электрическое поле Е = -Уф на границе фаз частица-электролит удовлетворяет общим электростатическим граничным условиям
Рис. 1. Двумерный гексагональный коллоидный кристалл. Частицы - бесконечно длинные стержни, перпендикулярные плоскости рисунка. Область определения численной задачи состоит из семи частиц с номерами от 0 до 6 вместе с соответствующими гексагональными ячейками Вигнера-Зейтца.
В рассматриваемой системе имеются только два компонента электролита (иона) с валентностями 2Х = +1 и г2 = -1. Их равновесные объемные концентрации могут быть обозначены как п (без индекса): п1 = п2 = п. Чтобы перейти к безразмерной форме уравнений и граничных условий, мы
вводим единицы к-1 = (2пд2е /££0кТ)-1/2 (длина Де-бая), кТ/де и (2п££0кТ)1/2 для длины, электрического потенциала и поверхностной плотности заряда соответственно. Начиная с этого места, в оставшейся части статьи подразумевается, что все величины, такие как х, у, ф и пр., являются безразмерными. Безразмерная форма уравнений для нашей системы записывается в виде
Э2ф + Э^ф _ |втЬф (в электролите) 2 л 2 [ 0 (внутри частицы).
д x ду
С учетом соотношения Е = -Уф, граничные условия в безразмерном виде для потенциала на границе фаз частица-электролит в нашей задаче можно записать в терминах электрического поля следующим образом:
р e — р e
n °part n part
(3)
En = a
(6)
Et = Et
part
(4)
Et = Et
part
(7)
где Еп и Ех - нормальная и тангенциальная компоненты электрического поля в электролите, Еп раг( и Е( раг( - соответствующие компоненты поля внутри частицы, £раг( - относительная диэлектрическая проницаемость материала частицы.
Нормальная компонента Еп раг( электрического поля внутри частицы исчезла из (6) благодаря условию £раг( = 0, принятому в нашем рассмотрении. Таким образом, нормальная компонента Еп электрического поля в электролите на границе фаз частица-электролит полностью определяется плотностью заряда а на поверхности частицы.
и
О
и
Следует отметить, что, тем не менее, электрическое поле внутри частицы в общем случае не равно нулю.
Область определения задачи для численного решения состоит из ячеек Вигнера-Зейтца произвольно выбранной центральной частицы 0 и ее шести ближайших соседей с номерами от 1 до 6. На внешней границе области определения задано стандартное граничное условие фон Неймана
(8)
En = 0.
- = 2 х
Э2 V
а, ß, N, M Zа' Nd Zß.N + M
Zа, NZß, N + M, (9)
^m = д2V aß dZ
ГУ n д Z,
(10)
а, NUZi ß, N + M
что отражается в обозначении C\
M
aß ,
вектор N не входит. По этой причине вектор N в (9) может быть выбран произвольно. Общая (бесконечная) матрица С квадратичной формы (9) состоит из подматриц См коэффициентов С^р, относящихся к данному вектору М. В двумерном случае матрицы См имеют размерность 2 х 2. Силовые постоянные обладают общими свойствами симметрии [19]:
Это граничное условие выполняется точно, когда все частицы находятся в своих положениях равновесия, и является только приближением, когда одна или более частиц отклоняются от положений равновесия. Метод определения силовых постоянных [14, 15], используемый в данной работе, требует решения уравнения ПБ для ряда смещений только центральной частицы 0 (см. раздел III). Это смещение не превышает 10% от расстояния a - 2R между ближайшими частицами и, следовательно, не искажает существенно распределение заряда на внешней границе, поэтому можно ожидать, что условие (8) выполняется с хорошей точностью.
Уравнения (5) вместе с граничными условиями (6)-(8) решались численно методом, описанным в [16, 17]. Он включает в себя решение дифференциальных уравнений методом конечных элементов и адаптивное перестроение сеток. Результатом решения является распределение ф = ф(х, у) электрического потенциала в области определения. После этого возможно вычисление других относящихся к делу величин. Нас будут интересовать силовые постоянные и давление.
III. ОПРЕДЕЛЕНИЕ СИЛОВЫХ
ПОСТОЯННЫХ И ДАВЛЕНИЯ
Согласно гармоническому приближению классическая потенциальная энергия кристалла в случае малых смещений частиц из положения равновесия есть квадратичная форма этих смещений [18, 19]:
_ ^ -M
Caß = Cßa >
M -M
Caß = Caß >
X C
= 0.
(11а) (116) (11в)
Из свойств (11) и (11) немедленно следует свойство симметричности
MM Caß = Cßa
(12)
матриц См для любого М.
Определение силовых постоянных основано на наблюдении, что первая производная ЭК/Э2а,№ а = х, у, с точностью до знака есть компонента ^„¿у силы, действующей на частицу N. Компоненты сил можно вычислить без использования численного дифференцирования с помощью тензора напряжений, который для рассматриваемой системы имеет вид
T = E
E - ( 2 E2
cos ф- 1 I
(13)
где Е = -Уф и I - матрица с единичной диагональю. Чтобы получить силовые постоянные, функциональную зависимость силы от смещения остается однократно продифференцировать по соответствующему смещению:
д ,,
а' (14)
.^M _ U1 а, 0
Caß" -dzz:
где у есть а-компонента смещения Z частицы из ее положения равновесия, задаваемого вектором У, а = х, у, в = х, у, N и М - векторы решетки Бравэ. Коэффициенты
квадратичной формы (9) называются силовыми постоянными. Благодаря трансляционной симметрии кристалла они не зависят от вектора N
которое
'Р, м
что прямо следует из определения (10) и где без ограничения общности положено N = 0.
Основываясь на свойствах симметрии, можно показать [14], что (14) эквивалентно следующему выражению для силовых постоянных, которое и использовалось при вычислениях:
дК м (15)
Р, о
где М - произвольный вектор решетки Бравэ. Согласно этому выражению силовая постоянная
м
Сар есть взятая со знаком минус а-компонента силы на единицу длины, действующей на частицу М при бесконечно малом смещении центральной частицы 0 в Р-направлении. Хотя способы определения силовых постоянных (14) и (15) эквива-
CM = _
C(lß д Z,
лентны, выражение (15) гораздо удобнее с практической точки зрения, поскольку требует смещения только частицы 0, а не всех присутствующих в задаче частиц М, как в (14). Помимо значительного сокращения объема вычислений это обстоятельство позволяет избежать смещений периферийных частиц, что делает более обоснованным наложение граничного условия фон Неймана (8) на внешней границе.
Компоненты сил, действующих на различные частицы, вычисляются путем численного интегрирования тензора напряжений (13) по замкнутому контуру, охватывающему частицу:
Е =
Е
Е - (1Е2
С08ф - 1 I
• п АГ, (16)
где п - внешняя нормаль к элементу АГ. Для минимизации численных ошибок в качестве контура использовалась граница Г соответствующей ячейки Вигнера-Зейтца, на которой электрическое поле мало. Для вычисления производной в (15) требуется определение зависимости ГаМ(2р0) в ряде точек 2р0 в окрестности положения равновесия. Из соображений симметрии следует, что достаточно смещений центральной частицы 0 только в положительном х-направлении. Сила определялась в десяти равноотстоящих точках вдоль положительного направления оси х так, что максимальное смещение составляло 10% от кратчайшего расстояния а - 2Я между частицами. По полученным значениям строилась полиномиальная аппроксимация функциональных зависимостей ГаМ(2х0) в окрестности нуля, коэффициент при линейном члене которой брался в качестве производной.
Давление в системе вычисляется с помощью тензора напряжений (13) следующим образом:
Силовые постоянные 0.8
(а)
0 3.5
3.0
2.5
2.0
1.5
1.0
0.5
0
70 60 50 40 30 20 10
0
Рис. 2. Силовые постоянные С.0, (1), С.1. 0) (2) и
Л1, 0) /оч
Су у (3) двумерного кристалла с гексагональной решеткой как функции плотности заряда о при фиксированной постоянной решетки а = 5.0, для частиц радиусом Я = 0.5 (а), 1.0 (б) и 2.0 (в).
а для частицы 1 с вектором (1, 0) вид
7 ° п
Г
Е
Е - ( 2 Е2
С08Ьф - 1 I
• пА1,(17)
где Ь - длина границы Г ячейки Вигнера-Зейтца. Можно показать [14], что выражение (17) в условиях данной задачи эквивалентно выражению для давления р, следующему из термодинамического соотношения р = -(ЭЕ/д¥)^т, где Г есть свободная энергия Гельмгольца, V - объем системы, а производная берется при постоянном числе N макроионов и постоянной температуре Т.
Матрица силовых постоянных для частицы 0, описываемой вектором (0, 0), имеет вид [14]
С0'0 = С(0'0)
С Схх
/ \
1 0 V 0 1 у
(18)
С
(1,0) _
С
(1,0)
0С
(1,0) УУ У
(19)
Матрицы остальных пяти ближайших соседей не являются независимыми; они связаны с матрицей С(1'0) законом преобразования матриц квадратичных форм при поворотах:
СМ = Ят (ф) СМ Я (ф),
(20)
где Я(ф) - матрица поворота на угол ф, к
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.