Гибридный метод решения уравнения Пуассона в области с металло-диэлектрическими углами

Во многих задачах электростатики и электродинамики требуется найти электромагнитное поле в сложных областях, содержащих металлические или диэлектрические входящие углы. В качестве примера подобных задач можно привести расчет распределения поля внутри многозазорной волноведущей системы, заполненной диэлектриками с различными коэффициентами диэлектрической проницаемости.

Наличие особенности электромагнитного поля, заключающееся в неограниченном росте некоторых его компонент вблизи металлического входящего угла, известно как теоретически [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)}$ слагаемого:

где $C_0=\mathop{\rm const}\nolimits$. В окрестности $\Pi$ заменим правую часть уравнения (1) на ее аппроксимацию (5) и рассмотрим вспомогательную задачу: найти частные ограниченные в окрестности нуля решения уравнения На лучах $\phi = \omega_1,\omega_L$ ставятся однородные граничные условия Дирихле На лучах, соответствующих линиям раздела сред с различными диэлектрическими проницаемостями, ставятся условия сопряжения

На внешней границе окрестности $\Pi$ в данный момент граничные условия не ставятся. В дальнейшем на ней будут задаваться условия непрерывности решения и его нормальной производной.

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

В областях строго внутри диэлектриков нетривиальные решения однородной задачи (6)--(9) в полярных координатах с центром в угловой точке будем искать в виде

В каждой из подобластей $\varepsilon(\phi)$ является постоянной, поэтому после применения процедуры разделения переменных решение в секторах представляется следующим образом:

где $A_i$ и $B_i$ --- некоторые константы. После подстановки решения (11) в условия сопряжения (8), (9) получается однородная СЛАУ относительно этих констант

где $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}$) подзадачу:

с дополнительными условиями (7)--(9). Решение исходной задачи (6) есть сумма решений задач (17), (7)--(9)

Представляя решение каждой из задач (17), (7)--(9) в виде $u_n=r^s\Phi_n(\phi)$, получим следующее уравнение:

откуда следует, что $s=n+2,$ и получается задача для $n$-й угловой части $\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)$:

где матрица с постоянными коэффициентами $\bigl\{A_{ij}^{(s)}\bigr\}$ определяется из условия равенства на границе $\partial\Pi_s$ функций $\chi_i^{(s)}(r,\phi)$ и $\widetilde{\psi}_i^{(s)}(x,y)$. Неизвестные коэффициенты перед введенными функциями выбираются равными $C_i^{(s)}.$ Дополнительно из стандартных конечных элементов на границе строится функция $\widetilde{\mu}^{(s)}(x,y)$ где коэффициенты $B_j^{(s)}$ определяются из условия ее равенства гладкой части решения $\mu^{(s)}(r,\phi)$ на границе $\partial\Pi_s$. Таким образом, в методе конечных элементов решение в области ${\Omega\backslash\Pi}$ ищется в виде разложения

где последняя сумма содержит $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).