Рассматривается двумерная задача для уравнения Пуассона в области, содержащей металло-диэлектрические углы звездного типа, которые могут быть как только металлическими, так и диэлектрическими. В их окрестности решение является ограниченным, а градиент имеет степенную особенность. Предлагается численный алгоритм решения данной задачи, основанный на методе конечных элементов и учитывающий асимптотическое поведение особенности решения в окрестности металло-диэлектрических углов.
Ключевые слова: уравнение Пуассона, металло-диэлектрический угол, метод конечных элементов.
УДК: 519.632.4, 51-73, 537.213. PACS: 02.70.Dh.
Во многих задачах электростатики и электродинамики требуется найти электромагнитное поле в сложных областях, содержащих металлические или диэлектрические входящие углы. В качестве примера подобных задач можно привести расчет распределения поля внутри многозазорной волноведущей системы, заполненной диэлектриками с различными коэффициентами диэлектрической проницаемости.
Наличие особенности электромагнитного поля, заключающееся в неограниченном росте некоторых его компонент вблизи металлического входящего угла, известно как теоретически [bb1,bb2,bb3,bb4,bb5], так и экспериментально [bb6,bb7,bb8]. Особенность поля также имеет место и при наличии диэлектрических углов [bb9]. В настоящей работе рассматривается случай, когда область содержит металло угол звездного типа, также приводящий к возникновению особенности поля.
Для расчета задач такого рода используются различные численные методы [bb10,bb11,bb12,bb13,bb14,bb15]. Однако наличие особенностей существенно ухудшает их скорость сходимости и точность. Одним из самых распространенных способов повышения эффективности численного метода является адаптивное сгущение сетки в окрестности особых точек, но в случае необходимости решения большого количества таких задач, например при решении задачи синтеза, данный подход приводит к значительным вычислительным затратам.
Другой подход заключается в использовании некоторой априорной информации о поведении решения, которая закладывается в численный алгоритм. В качестве примера можно привести метод конечных элементов [bb16,bb17] с сингулярными элементами, описывающими особенность [bb18]. Однако добавление новых элементов к базису может приводить к ухудшению свойств его матрицы Грама [bb16] и снижению эффективности метода.
Численный алгоритм, предложенный в настоящей работе, основан на построении явного вида особенности решения в окрестности металло углов и его дальнейшем использовании в методе конечных элементов. Однако для преодоления описанных выше проблем дополнительные конечные элементы не добавляются, а заменяют часть базисных функций внутри и на границе малых окрестностей углов, при этом задача решается уже вне этих окрестностей на достаточно грубой сетке. Использование дополнительной информации о поведении решения делает алгоритм быстрым и эффективным и позволяет с высокой точностью находить сингулярную часть градиента решения.
Рассмотрим область $\Omega = \bigcup\limits_{i=1}^3\Omega_i$ с металлической границей, заполненную диэлектриками с различными диэлектрическими проницаемостями в каждой из подобластей $\Omega_i$ (рис. 1). В $\Omega$ присутствуют угловые точки различных типов: металлический угол ($M$), внутренние диэлектрические углы ($D$) и металло углы ($C$).
Распределение электрического потенциала $u$ описывается уравнением Пуассона с граничными условиями Дирихле и условиями сопряжения на границе раздела диэлектриков
где $\varepsilon_i=\mathop{\rm const}\nolimits,$ $M\in\Omega_i$, $\partial\Omega$ --- общая внешняя металлическая граница, $\partial\Omega_{ij}$ --- граница между подобластями $\Omega_i$ и $\Omega_j$, $f(M)$ --- достаточно гладкая функция, выражающаяся через плотность зарядов $\rho$: $f(M)=4\pi\rho(M)$.
Для численного решения задачи получим асимптотику решения по гладкости в окрестности всех указанных типов углов.
Построим асимптотику решения задачи (1)-(4) в достаточно малой окрестности $\Pi$ металло угла с вершиной $M_0$ (рис. 2).
Пусть из вершины угла выходит $L$ лучей. Пронумеруем их таким образом, чтобы первый и $L$-й лучи соответствовали металлическим границам. Введем полярную систему координат с центром в этой вершине и направим полярную ось вдоль $L$-го луча. За $\omega_i,$ $i=1,\dots,L,$ $\omega_L=2\pi$ обозначим углы между лучами и полярной осью.
Разложим в окрестности $\Pi$ правую часть $f(M)$ уравнения (1) в ряд Тейлора по $r$. Рассмотрим частичную сумму ряда, состоящую из ${(P+1)}$ слагаемого:
На внешней границе окрестности $\Pi$ в данный момент граничные условия не ставятся. В дальнейшем на ней будут задаваться условия непрерывности решения и его нормальной производной.
Поскольку задача (6)--(9) является линейной и неоднородной, представим ее решение в виде суммы частного решения неоднородной задачи и общего решения однородной задачи, каждое из которых удовлетворяет условиям (7)--(9).
В областях строго внутри диэлектриков нетривиальные решения однородной задачи (6)--(9) в полярных координатах с центром в угловой точке будем искать в виде
В каждой из подобластей $\varepsilon(\phi)$ является постоянной, поэтому после применения процедуры разделения переменных решение в секторах представляется следующим образом:
где $i=1,\dots,L{-}1$ и уравнения (12) соответствуют металлической границе, а (13)--(14) --- границе раздела двух диэлектриков.
Характеристическое уравнение относительно $\sigma$, полученное из условия равенства нулю определителя системы (12)--(14), определяет нетривиальные решения данной системы. В случае если существуют положительные корни характеристического уравнения ${\sigma<1},$ то в угловой точке градиент решения будет иметь особенность, что с физической точки зрения соответствует неограниченному росту напряженности электрического поля.
Отдельно следует рассмотреть значение $\sigma_0=0,$ для которого система (12)--(14) имеет несколько другой вид. В этом случае общее решение будет выражаться через линейную функцию угла:
однако из условий (7)--(9) и того, что ${\varepsilon_i>0},$ следует, что существует только тривиальное решение.
Пусть $\sigma_i$ --- положительные корни характеристического уравнения системы (12)--(14), в этом случае решение в точке $M_0$ будет ограничено. Обозначим соответствующие решения системы (12)--(14) с учетом возможной кратности корня как
Сингулярная часть решения задачи (6)--(9) представляется как линейная комбинация функций вида (15)
Рассмотрим неоднородную задачу (6)--(9). В силу линейности уравнения (6) получим (${P+1}$) подзадачу:
Представляя решение каждой из задач (17), (7)--(9) в виде $u_n=r^s\Phi_n(\phi)$, получим следующее уравнение:
Набор задач (18)--(20) решается численно.
Таким образом, приближенное частное решение неоднородной задачи имеет следующий вид:
Решение исходного уравнения (6)--(9) в окрестности угла $\Pi$ представляется как
где первое слагаемое описывает гладкую часть решения, а конечная сумма содержит слагаемое с особенностью в угловой точке.
Получение асимптотики решения для диэлектрического угла, лежащего внутри области $\Omega,$ в целом повторяет аналогичный вывод для металлического, за исключением некоторых деталей.
Ось полярной системы координат в данном случае вводится по направлению любой из линий разрыва диэлектрической проницаемости. Граничное условие (7) исключается, а условие сопряжения при $\phi=0.2\pi$ в соответствии с условиями периодичности принимают следующий вид:
При $\sigma_0=0$ однородная задача, в отличие от случая с металлическим углом, имеет нетривиальное решение --- константу, которое должно быть учтено в разложении (16).
Для произвольной конфигурации диэлектрического угла ФСР системы (12)--(14) при различных $\sigma_i$, как правило, состоит из одного столбца, но при углах, кратных $\pi,$ возможно появление кратных корней характеристического уравнения и ФСР, состоящих из нескольких столбцов.
Для диэлектрического угла с двумя диэлектриками поведение особенности вблизи вершины было получено в [bb9] с помощью метода Кондратьева. Показатели $\sigma_i$ степени $r$, полученные из характеристического уравнения системы (13), (14), соответствуют значениям, полученным в работе [bb9].
Вернемся к исходной задаче (1)--(4). Пусть в области $\Omega$ есть $S$углов, в которых имеет место особенность градиента, и для каждого из них получено представление решения в виде (22). Выделим окрестности углов $\Pi_s$ (рис. 3), в которых приближенное решение задачи (1)--(4) ищется в виде конечной линейной комбинации функций со степенями $r$ не выше $L$:
где $I_{\rm max}$ определяется из условия ${\sigma_i<L}$. На введенных границах $\partial\Pi_s$ поставим условия сшивания решения, заключающееся в его непрерывности и непрерывности его нормальной производной.
Умножим уравнение (1) на пробную функцию $v$ и проинтегрируем полученное выражение по области ${\Omega\backslash\Pi,}$ где $\Pi=\bigcup\limits_{s=1}^S\Pi_s$:
Так как граничное условие (2) является главным, то ${v\big|_{\partial\Omega}=0}$ и первый интеграл по границе $\partial\Omega$ обращается в нуль. Из условия сшивания решения на границах $\partial\Pi_s$ и его представления (23) следует
Таким образом, исходная задача, решаемая в области $\Omega$, сводится к решению задачи (26) в области ${\Omega\backslash\Pi.}$
Приближенное решение (26) ищется методом конечных элементов. Обозначим через $\{\omega_i\}_{i=1,\dots,N}$ узлы сетки в области ${\Omega\backslash\Pi,}$ полученной с помощью триангуляции (рис. 4). В качестве базиса берутся лагранжевы элементы первого порядка $\{\psi_i(x,y)\}_{i=1,\dots,N.}$ При этом к столбцу неизвестных значений коэффициентов в разложении решения по базису МКЭ, являющихся значениями решения в узлах сетки, добавляются неизвестные коэффициенты $C_i^{(s)}.$ Условие сшивания нормальной производной на границе $\partial\Pi$ уже учтено при выводе задачи (26), поэтому при ее решении необходимо удовлетворить только условию непрерывности, являющемуся главным.
Условие непрерывности решения на границах $\partial\Pi_s$ реализуется следующим образом. На границу $\partial\Pi_s$ помещается $K_s$ узлов $\{\widetilde\omega_i\}_{i=1,\dots,K_s}.$ Так как решение задачи в подобласти $\Pi_s$ определяется формулой (23), то вместо стандартных лагранжевых элементов $\{\psi_i(x,y)\}_{i=1,\dots,N}$ с вершинами в $\{\widetilde\omega_i\}_{i=1,\dots,K}$ на границе берутся их линейные комбинации $\widetilde\psi_i^{(s)}(x,y)$:
где последняя сумма содержит $N-K=N-\sum\limits_{s=1}^SK_s$ исходных базисных функций с вершинами строго внутри области ${\Omega\backslash\Pi}$. Поскольку на границах $\partial\Pi_s$ все ${\psi_j(x,y)=0,}$ то для построенного решения из конечных элементов условие непрерывности на $\partial\Pi_s$ выполняется автоматически.
Сохранив старые обозначения, объединим два семейства базисных функций $\bigl\{\widetilde{\psi}_i^{(s)}(x,y)\bigr\}_{i=0,\dots,I_{\rm max}}^{s=1,\dots,S}$ и $\{\psi_j(x,y)_{j=1,\dots,N{-}K}\}$ в $\{\psi_i(x,y)\}_{i=1,\dots,\widetilde{N}}$ и соответствующие им столбцы неизвестных коэффициентов $C_i^{(s)}$ и $A_j$ в $\{C_i\}_{i=0,\dots,\widetilde{N}}$, где $\widetilde{N}=N-K+I_{\rm max}+1$ --- общее количество базисных функций. Приближенное решение (29) запишется в виде
После подстановки (30) и (31) в (26) получается система уравнений относительно столбца неизвестных коэффициентов ${\bf C}=\{C_i\}_{i=0}^{\widetilde N},$ которая имеет вид
Отметим, что при решении полученной задачи (32)--(33) также находятся коэффициенты при функциях, описывающих сингулярную часть решения в каждой окрестности $\Pi_s$.
Для проверки предложенного метода задача (1)--(4) дополнительно решалась во всей области с помощью стандартного метода конечных элементов без учета особенности в угловых точках на более густой сетке.
На рис. 5 представлено решение задачи (1)--(4) для ${f(x,y)=1,}$ ${\varepsilon_1=1,}$ ${\varepsilon_2=2,}$ ${\varepsilon_3=3}$ стандартным методом конечных элементов на сетке из 11 992 узлов с сильным сгущением в окрестности рассматриваемых угловых точек.
На рис. 6 представлена относительная погрешность решения на мелкой сетке без выделения особенности к решению, полученному предложенным методом на грубой сетке с 1264 узлами. Расхождение в окрестности углов не превышает 0.4 %. При использовании такой же грубой сетки в обычном методе конечных элементов без выделения сингулярной части решения погрешность в окрестности углов достигает 20 %, что при вычислении градиента приводит к существенным ошибкам в значении напряженности электрического поля.
В настоящей работе представлен алгоритм расчета электростатической задачи для уравнения Пуассона в области, содержащей металло углы. Предложенный алгоритм обладает высокой эффективностью и позволяет с хорошей точностью получить распределение электрического поля вблизи угловых точек. Также данный метод обладает хорошей гибкостью, позволяя единообразно обрабатывать металлические, диэлектрические и металло углы с различными конфигурациями без существенных изменений в способе решения.
Метод может быть использован для решения широкого класса задач электростатики и электродинамики, таких как задачи анализа и синтеза в областях со сложной геометрией и заполнением, содержащих угловые точки различных типов.
Работа М. И Светкина выполнена при финансовой поддержке РФФИ (грант 16-01-00690). Работа А. И. Ерохина выполнена при финансовой поддержке РФФИ (грант 16-31-60084-мол_а_дк). Работа А. Н. Боголюбова и И. Е. Могилевского выполнена при финансовой поддержке РФФИ (15-01-03524).
A hybrid method for the solution of the Poisson equation in the domain of metal--dielectric corners
A. N. Bogolyubov, A. I. Erokhin, I. E. Mogilevsky, M. I. Svetkin
Department of Mathematics, Faculty of Physics, Lomonosov Moscow State University.
Moscow 119991, Russia.
E-mail: $^a$mihail-svetkin@mail.ru.
The Poisson’s equation in a two-dimensional domain that contains star metal--dielectric corners is considered. The corners can be only metal or only dielectric. In their neighborhood, the solution is bounded, but the gradient has power-law singularities. A numerical algorithm for solving this problem is proposed. The algorithm is based on the finite-element method and takes the asymptotic behavior of the solution in the neighborhood of metal-dielectric corners into account.
Keywords: Poisson equation, metal--dielectric corner, finite-element method.
PACS: 02.70.Dh.
Received 27 June 2016.
English version: Moscow University Physics Bulletin. 72, No. Pp.
footnotesize
Сведения об авторах
