Доклады Российской академии наук. Математика, информатика, процессы управления, 2022, T. 507, № 1, стр. 40-45

1. Задачи электростатики в пространственно периодических средах возникают естественным образом в ряде актуальных исследований в физической химии и материаловедении, где электрическое поле применяется в качестве способа управления состоянием системы [1]. Численное решение задач электростатики в основном осуществляется тремя группами методов, куда входят (а) разностные методы; (б) проекционные методы, такие как различные версии метода Галеркина; (в) метод интегральных уравнений [2–5]. Наряду с этими методами в ряде специфических задач используются спектральные методы, в которых искомое решение представляется в виде разложения по системе собственных функций соответствующего дифференциального оператора [4]. Спектральный подход эффективен в задачах с периодическими граничными условиями в однородной среде, поскольку в этом случае аналитическое решение представимо в виде разложения по тригонометрической системе функций. Однако дело обстоит иначе в случае неоднородных сред, так как тригонометрические функции, вообще говоря, не являются собственными функциями дифференциального оператора задачи. Тем не менее тригонометрический базис в пространственно периодических задачах все еще является достаточно привлекательным для представления решения в виде линейной комбинации тригонометрических функций, когда характеристики среды являются плавно меняющимися функциями. В настоящей работе предложен аналог псевдоспектрального подхода, примененного к задачам теории полимеров [6], для решения задач электростатики в пространственно периодических диэлектрических средах с плавно меняющейся диэлектрической проницаемостью, на основе которого создан высокоточный экономичный итерационный метод вычисления потенциала и напряженности электрического поля с использованием быстрого преобразования Фурье.

2. Потенциал $\Phi ({\mathbf{r}})$ электрического поля в среде с достаточно гладкой положительной диэлектрической проницаемостью $\varepsilon ({\mathbf{r}})$, где ${\mathbf{r}} = (x,y,z)$ – радиус-вектор, удовлетворяет уравнению в частных производных эллиптического типа [1]

В случае неоднородного диэлектрика, помещенного в первоначально однородное электрическое поле, удобно представить искомый потенциал $\Phi ({\mathbf{r}})$ в виде

где E₀ – вектор напряженности электрического поля в однородной диэлектрической среде с диэлектрической проницаемостью, равной среднему значению $\bar {\varepsilon }$ диэлектрической проницаемости среды $\varepsilon ({\mathbf{r}})$, а $\varphi ({\mathbf{r}})$ – добавочный потенциал. Подстановка правой части (2) в уравнение (1) приводит к следующему уравнению

Искомый добавочный потенциал $\varphi ({\mathbf{r}})$ в правой части (2) удовлетворяет периодическим граничным условиям в пространственно периодической диэлектрической среде

где в качестве функции f используются функции $\varepsilon ({\mathbf{r}})$ и $\varphi ({\mathbf{r}})$, Z – множество целых чисел, $\{ {{{\mathbf{a}}}_{k}}\} $ – базис в E³.

Отметим, что уравнение Пуассона относительно функции $\varphi ({\mathbf{r}})$

в котором функции $\varphi ({\mathbf{r}})$ и $\rho ({\mathbf{r}})$ удовлетворяют условиям периодичности (4), определяет решение с точностью до постоянной составляющей. Поэтому в подпространстве потенциалов $\varphi ({\mathbf{r}})$ с нулевым средним значением уравнение (5) имеет единственное решение. Для такого решения задачи (5) с нулевой постоянной составляющей и граничными условиями (4) будем использовать обозначение

Уравнение (3) определяет добавочный потенциал $\varphi ({\mathbf{r}})$ также с точностью до постоянной составляющей, значение которой без ограничения физической общности можно положить равным нулю. Подставляя правую часть уравнения (3) в уравнение (6) вместо функции $\rho ({\mathbf{r}})$, получаем еще одну форму записи уравнения (3)

В отличие от уравнения (3) уравнение (7) подразумевает, что постоянная составляющая потенциала $\varphi ({\mathbf{r}})$ равна нулю.

Для решения уравнения (7) с периодическими граничными условиями можно применить итерационный метод Пикара с положительным параметром τ

или другие, более быстро сходящиеся итерационные методы, в том числе многослойные [7, 8].

3. Рассмотрим экономичный способ вычисления оператора $F[\varphi ({\mathbf{r}});\varepsilon ({\mathbf{r}})]$ с помощью быстрого преобразования Фурье (БПФ). С этой целью введем взаимный (обратный) базис $\{ {{{\mathbf{b}}}_{j}}\} $, определяемый условиями

где ${{\delta }_{{kj}}}$ – символ Кронекера. Каждая достаточно гладкая действительная периодическая функция $f({\mathbf{r}})$ может быть разложена в равномерно сходящийся тригонометрический ряд Фурье, допускающий почленное дифференцирование (т.е. ряды Фурье для производных также сходятся равномерно), где ${\mathbf{m}} = ({{m}_{1}},{{m}_{2}},{{m}_{3}})$ – целочисленный вектор, Z – множество целых чисел, а коэффициенты ряда (10) определяются равенствами

Следующие дифференциальные операторы (градиент и лапласиан), примененные к достаточно гладкой функции $f({\mathbf{r}})$, также выражаются в виде равномерно сходящихся рядов Фурье

Обращение лапласиана на множестве функций с нулевым средним значением происходит с помощью равенства

Соотношения (10)–(13) лежат в основе высокоточного вычисления градиента и лапласиана функции на равномерной сетке с помощью дискретного преобразования Фурье, когда высокочастотными гармониками функции можно пренебречь.

Для вычисления коэффициентов Фурье (12) воспользуемся формулой прямоугольников на равномерной сетке в параллелепипеде Π

где полагая где N обозначает конечное множество целочисленных векторов m и n, компоненты которых принимают целые неотрицательные значения, описанные соотношениями (15). Отметим тесную связь формулы прямоугольников на отрезке прямой для тригонометрических полиномов и квадратурных формул Гаусса для полиномов Чебышева [9], а именно: формула прямоугольников на равномерной сетке с числом узлов N является точной для тригонометрических полиномов степени n, где $n < N{\text{/}}2$ и, тем самым, выполняет роль высокоточной формулы Гаусса на множестве периодических функций.

Равенство (16) представляет собой с точностью до множителя дискретное преобразование Фурье сеточной функции $f({{{\mathbf{r}}}_{{\mathbf{n}}}})$. Восстановление сеточной функции $f({{{\mathbf{r}}}_{{\mathbf{n}}}})$ по ее коэффициентам (16) происходит с помощью обратного дискретного преобразования Фурье

Соотношения (16) и (17) позволяют использовать алгоритм БПФ для высокоточного и экономичного вычисления коэффициентов Фурье и восстановления по ним исходной сеточной функции. Для вычисления коэффициентов Фурье дискретных аналогов дифференциальных операторов (12) градиента $\tilde {\nabla }{{f}_{{\mathbf{m}}}}$, лапласиана ${{\tilde {\nabla }}^{2}}{{f}_{{\mathbf{m}}}}$ и обращения лапласиана (13) ${{\tilde {\nabla }}^{{ - 2}}}{{f}_{{\mathbf{m}}}}$ сеточной функции $f({{{\mathbf{r}}}_{{\mathbf{n}}}})$ применяются соответственно соотношения

и где

Дискретные аналоги градиента, лапласиана и обращения лапласиана сеточной функции $f({{{\mathbf{r}}}_{{\mathbf{n}}}})$ вычисляются по формуле (17) обратного дискретного преобразования Фурье с заменой коэффициентов  f_m соответственно на коэффициенты $\tilde {\nabla }{{f}_{{\mathbf{m}}}}$, ${{\tilde {\nabla }}^{2}}{{f}_{{\mathbf{m}}}}$ и ${{\tilde {\nabla }}^{{ - 2}}}{{f}_{{\mathbf{m}}}}$, определенные выше в (18)–(21).

Отметим, что в сравнении с широко распространенными разностными аппроксимациями дифференциальных операторов [2], введенные здесь их дискретные аналоги (18)–(21) обладают определенным преимуществом. Они совпадают с исходными операторами (12) и (13) в точках сетки (14), (15) на множестве конечных линейных комбинаций тригонометрических функций вида (10), если компоненты индексов гармоник удовлетворяют неравенствам $2{{m}_{k}} < {{N}_{k}}$ ($k = 1,\;2,\;3$), причем сами эти линейные комбинации могут быть использованы для вычисления функций и соответствующих операторов (12) и (13) во всех точках параллелепипеда Π, а не только в узлах сетки.

Замечание. Стремление к нулю коэффициентов разложения в тригонометрический ряд Фурье функции lnε(r) в уравнении (7) должно быть достаточно быстрым. Для простоты понимания условий, накладываемых на спектр функции lnε(r), рассмотрим одномерное уравнение (7). Если шаг равномерной сетки по пространственной переменной равен h, то пространственная частота Найквиста равна ${{(2h)}^{{ - 1}}}$. Следовательно, выбор шага сетки h должен производиться так, чтобы частичная тригонометрическая сумма ряда Фурье векторной функции $\nabla \ln \varepsilon ({\mathbf{r}})$ с гармониками, включающими волновые числа q, удовлетворяющие неравенству ${\text{|}}q{\text{|}} \leqslant \pi {{h}^{{ - 1}}}$, аппроксимировала исходный ряд и его производную с достаточной точностью, определяемой конкретным приложением. Отметим, что в правую часть уравнения (7) входит скалярное произведение градиентов функции $\ln \varepsilon ({\mathbf{r}})$ и искомого добавочного потенциала $\varphi ({\mathbf{r}})$. Если каждая из этих функций содержит гармоники с модулем волнового числа $\pi {{h}^{{ - 1}}}$, то их скалярное произведение может содержать гармоники с модулем волнового числа $2\pi {{h}^{{ - 1}}}$, т.е. включать пространственные частоты, превышающие частоту Найквиста. Поэтому ограничение на шаг сетки должно быть усилено по меньшей мере вдвое: значимые гармоники функции $\ln \varepsilon ({\mathbf{r}})$ (т.е. те, которые передают существенную информацию о функции) должны включать только волновые вектора с длинами ${\text{|}}q{\text{|}} \leqslant \pi {{(2h)}^{{ - 1}}}$. В случае многомерного уравнения (7) под q следует понимать каждую из независимых компонент волнового вектора ${{{\mathbf{q}}}_{{\mathbf{m}}}}$.

Кратко опишем процедуру вычисления правой части уравнения (7):

1. Последовательно полагаем функцию $f({{{\mathbf{r}}}_{{\mathbf{n}}}})$ равной $\varphi ({{{\mathbf{r}}}_{{\mathbf{n}}}})$ и $\ln \varepsilon ({{{\mathbf{r}}}_{{\mathbf{n}}}})$, применяя затем к ней каждый раз дискретное преобразование Фурье (16) для вычисления коэффициентов Фурье f_m, находя векторные коэффициенты $\tilde {\nabla }{{f}_{{\mathbf{m}}}}$ в соответствии с левой формулой (18) и вычисляя дискретный аналог градиента $\tilde {\nabla }f({{{\mathbf{r}}}_{{\mathbf{n}}}})$ по формуле (17), в которой скалярные коэффициенты f_m заменены на векторные коэффициенты $\tilde {\nabla }{{f}_{{\mathbf{m}}}}$.

2. Вычисляем скалярное произведение ρ(r_n) = =  $({{{\mathbf{E}}}_{0}}\, - \,\tilde {\nabla }\varphi ({{{\mathbf{r}}}_{{\mathbf{n}}}}),\tilde {\nabla }{\text{ln}}\varepsilon ({{{\mathbf{r}}}_{{\mathbf{n}}}}))$, где вектора $\tilde {\nabla }\varphi ({{{\mathbf{r}}}_{{\mathbf{n}}}})$ и $\tilde {\nabla }{\text{ln}}\varepsilon ({{{\mathbf{r}}}_{{\mathbf{n}}}})$ уже найдены на предыдущем шаге, и полагаем функцию $f({{{\mathbf{r}}}_{{\mathbf{n}}}})$ равной $\rho ({{{\mathbf{r}}}_{{\mathbf{n}}}})$, применяя к ней дискретное преобразование Фурье (16) для вычисления коэффициентов Фурье  f_m, находя величины ${{\tilde {\nabla }}^{{ - 2}}}{{f}_{{\mathbf{m}}}}$ в соответствии с формулой (19) и вычисляя дискретный аналог обращения лапласиана ${{\tilde {\nabla }}^{{ - 2}}}f({{{\mathbf{r}}}_{{\mathbf{n}}}})$ = = ${{\tilde {\nabla }}^{{ - 2}}}\rho ({{{\mathbf{r}}}_{{\mathbf{n}}}})$ по формуле (17), в которой коэффициенты  f_m заменены на коэффициенты ${{\tilde {\nabla }}^{{ - 2}}}{{f}_{{\mathbf{m}}}}$.

Результатом применения указанных выше двух шагов является вычисление значений правой части $F[\varphi ({\mathbf{r}});\varepsilon ({\mathbf{r}})]$ уравнения (7) в узлах сетки (14). Далее для численного решения уравнения (7) следует применить итерационный процесс (8) или какой-либо иной метод, например, из работ [7, 8].

Вычислительная трудоемкость прямого (16) и обратного (17) дискретных преобразований Фурье методом БПФ оценивается величиной порядка $O(M\log M)$ операций, где $M = {{N}_{1}}{{N}_{2}}{{N}_{3}}$. Указанная оценка и возможность эффективного распараллеливания алгоритмов БПФ [10] предоставляют возможность численного решения трехмерных задач электростатики в периодических диэлектрических средах итерационными методами с использованием большого числа узлов сетки в параллелепипеде Π за реальное время.

4. Применим предложенный выше алгоритм для расчета электрического поля в модельных одномерной и трехмерной пространственно периодических диэлектрических средах, состоящих из компонентов А и В, каждый из которых характеризуется значением диэлектрической проницаемости ${{\varepsilon }_{A}}$ и ${{\varepsilon }_{B}}$ соответственно.

Сначала рассмотрим одномерную среду с ламеллярной морфологией, положив объемные доли компонентов А и В соответственно равными

Локальную диэлектрическую проницаемость среды определим следующим образом

Однородное электрическое поле ${{E}_{0}} = 1$ считаем направленным вдоль оси x. Уравнение (3) относительно добавочного потенциала $\varphi (x)$ в этом случае принимает вид

с периодическим граничным условием $\varphi (x + 1)$ = = φ(x).

Для применения вычислительного алгоритма, описанного в пп. 2 и 3, рассмотрим равномерную сетку ${{x}_{n}} = n{\text{/}}N$ ($n = 0,1,...,N - 1$) на отрезке [0, 1], где $N\, = \,{{2}^{k}}$ ($k\, = \,1,2,...,10$). Отметим, что коэффициенты разложения функции (lnε(x))_x = ${{\varepsilon }_{x}}(x){\text{/}}\varepsilon (x)$ в тригонометрический ряд Фурье (по функциям $\exp (i2\pi nx)$) достаточно быстро стремятся к нулю, так что при выполнении неравенства ${\text{|}}n{\text{|}} \leqslant 20$ частичная сумма этого ряда аппроксимирует функцию ${{(\ln \varepsilon (x))}_{x}}$ с точностью 12–13 значащих цифр. Дальнейшее увеличение числа членов ряда при стандартной длине мантиссы 16 действительного числа нецелесообразно, поскольку оно не приводит к увеличению точности вычислений, так как гармоники с большими значениями номеров содержат только шум с характерными значениями амплитуд высших гармоник порядка ~10^–15. Суммирование таких рядов представляет собой хорошо известную некорректно поставленную задачу [11]. Анализ гармоник численного решения $\varphi (x)$ и его первой производной ${{\varphi }_{x}}(x)$ показывает, что максимальная достигаемая точность вычислений этих величин составляет 12–14 значащих цифр. Поэтому для вычисления коэффициентов Фурье данных функций достаточно использовать сетки с числом узлов $N = 32,64$ и использовать частичные суммы их тригонометрических рядов, включающие гармоники с номерами ${\text{|}}n{\text{|}} \leqslant 20$. Сетки с большим числом N узлов не улучшают точность вычислений. В качестве иллюстраций на рис. 1 приведены графики напряженности электрического поля (рис. 1а)

и плотности энергии электрического поля (рис. 1б)

Перейдем к расчету электрического поля в трехмерной пространственно периодической диэлектрической среде, имеющей морфологию алмаза [12]. Объемные доли компонентов А и В, использованные в расчетах, имеют вид

а диэлектрическая проницаемость среды определяетя равенством

С помощью предложенного выше численного алгоритма вычисляем добавочный потенциал $\varphi ({\mathbf{r}})$, удовлетворяющий уравнению (3) и условиям периодичности

Для определения численного решения использовалась равномерная сетка в единичном кубе с числом узлов 32³ и 64³. Точность вычисления напряженности электрического поля в любой точке куба составляла 6 и 10 значащих цифр соответственно. Для численной сходимости итераций с относительной точностью 10^–13 потребовалось 9 итераций по методу Ng [7] или 12 итераций метода последовательных приближений (8) при τ = 1.

На рис. 2 представлена поверхность уровня плотности энергии электрического поля

где $\overline {{{\rho }_{E}}({\mathbf{r}})} $ обозначает среднее значение ${{\rho }_{E}}({\mathbf{r}})$ в единичном кубе.

В заключение еще раз отметим, что при фиксированном числе узлов сетки точность численного решения существенно зависит от скорости убывания коэффициентов Фурье функции $\ln \varepsilon ({\mathbf{r}})$ с ростом номеров гармоник. Чем выше скорость убывания коэффициентов Фурье, тем более точно вычисляется решение.

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

КОНФЛИКТ ИНТЕРЕСОВ

Авторы заявляют об отсутствии конфликта интересов.