The accuracy of numerical methods for the analysis of electrostatic fields

Full Text

Введение

Существуют многочисленные методы решения задач оценки. Эти задачи имеют первостепенное значение при решении вопросов защиты, построении экранов, оценки безопасностей уровней электрического поля и других [1–4].

Большинство из этих методов опирается на численные расчеты, поскольку аналитические методы имеют узкую область решений, так как граничные условия для таких задач должны совпадать с координатными плоскостями и уровнями. Как правило, аналитические методы решаются методом разделения переменных. Разумеется, точность таких методов абсолютна. Однако для решения прикладных задач они плохо приспособлены, поскольку реальные задачи имеют сложную геометрию, кусочно-линейное представление параметров среды, в частности электрических и магнитных.

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

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

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

Вариация. Основная идея метода конечных элементов была предложена Р. Курантом в 1943 г. и основана она на использовании функциональной трактовки законов электростатики. Связывая функционал с дифференциальным уравнением Лапласа можно получить эффективный способ решения задачи о распределении электростатического поля.

Задача определения стационарного электрического поля в общем случае сводится к минимизации функционала вида [5]:

$Ф\left[\text{φ}\left(x,y\right)\right]=\underset{\text{Ω}}{\overset{}{}}F\left[x,y,\text{φ},{\text{φ}}_x,{\text{φ}}_y\right]dxdy,$ (1)

где $\text{φ}\left(x,y\right)$ – функция потенциала электрического поля;

${\text{φ}}_x$, ${\text{ φ}}_y$ – частные производные по координатам, соответственно, x и y;

Подынтегральная функция $F\left[\dots \right]$ имеет непрерывные производные до второго порядка включительно.

Искомая функция $\text{φ}$ на границе Г области $\text{Ω}$ должна удовлетворять значениям граничной функции $f$.

Экстремум такого функционала находится на основе определения функциональной производной в стационарной точке, то есть в такой точке, где функционал достигает своего минимума или максимума. Здесь под точкой понимается такое ее обобщение, при котором для определения этой точки требуется бесконечное количество координат [6], иными словами это просто функция двух переменных в двумерном гильбертовом пространстве $H^2$. При этом функция $\text{φ}$ определена в области $\text{Ω}$. Графическое представление $\text{φ}\left(x,y\right)$ дано на Рис. 1.

 

Рис. 1 - График функции потенциала электрического поля

 

Для вычисления функциональной производной воспользуемся ее определением [7, 8]:

$\underset{\begin{array}{c}diam\text{ ω}\to 0\\ \begin{array}{c}\max \\ \text{ξ},\text{ ζ}\in \text{ω }\end{array}\begin{array}{c}\text{ δφ}\to 0\\ \end{array}\end{array}}{\lim }\frac{\text{Ф}\left[\text{φ}+\text{δφ}\right]-\text{Ф}\left[\text{φ}\right]}{{\displaystyle {\iint }_{\text{Ω}}\text{δφ}dxdy}}=\frac{\text{δФ}\left[\text{φ}\right]}{\text{δφ}}.$ (2)

Здесь$\text{ diam ω}$ означает диаметр многообразия $\omega$, определяемый как

$diam\text{ω}=su{p}_{a,b\in \text{ω}}d\left(a,b\right),$ (3)

где $d$ – квадратичная метрика.

Вариация $\text{δφ}$ имеет игольчатый характер, представленный на Рис. 1 элементом $\left(\text{φ}+\text{δφ}\right)$.

Применение операции (2) для вычисления экстремального значения функционала (1) составляет суть теоремы Остроградского-Эйлера.

Воспользовавшись результатом этой теоремы, можем записать соответствующее уравнение Остроградского-Эйлера:

$\frac{\partial F}{\partial \text{φ}}-\frac{\partial }{\partial \text{x}}\left(\frac{\partial F}{\partial {\text{φ}}_x}\right)-\frac{\partial }{\partial \text{y}}\left(\frac{\partial F}{\partial {\text{φ}}_y}\right)=0.$ (4)

В частном случае, когда подынтегральная функция имеет следующий вид:

$F=\frac{1}{2}\left[{\left(\frac{\partial \text{φ}}{\partial x}\right)}^2+{\left(\frac{\partial \text{φ}}{\partial y}\right)}^2\right]$ (5)

уравнение Остроградского-Эйлера становится эквивалентным уравнению Лапласа:

$\frac{{\partial }^2\text{φ}}{\partial x^2}+\frac{{\partial }^2\text{φ}}{\partial y^2}=\mathrm{0,}$ (6)

с граничным условием

${\left.\text{φ}\right|}_{\text{Г}}=f.$ (7)

Конечные элементы. Предположим, что носитель функции $\text{φ}$ имеет вид треугольного симплекса в двумерном пространстве (плоскость x0y на Рис. 2).

 

Рис. 2 - Линейная потенциальная функция $\phi \left(x,y\right)$

 

В вершинах симплекса потенциал будем считать на первых порах заданным и равным ${\phi }_1, {\text{φ}}_2{\text{, φ}}_3$. Тогда текущее значение потенциала в точках, принадлежащих носителю $\left(x, y\right)\in ∆\left(x,y\right)\in ∆$, определим линейной формой

$\text{φ}=a+bx+cy.$ (8)

Для определения коэффициентов $a, b и c$ линейной формы (8) необходимо составить систему уравнений для угловых точек симплекса и разрешить ее относительно неизвестных коэффициентов [5].

Затем произведя покрытие всей области $\text{Ω}$ треугольными элементами (Рис. 1), с учетом равенства потенциалов на границах смежных элементов получаем два множества точек. В точках, которые лежат на границе Г, потенциалы заданы и определяются функцией (7), а потенциалы внутри области  неизвестны и находятся исходя из минимума функционала (1).

Отметим, что подстановка линейных функций  в функционал (1) превращает его в функцию N=n+m переменных, где N – количество узловых точек конечных элементов, n – неизвестные потенциалы,
m – известные. В частном случае вида (5) эта функция является положительно определенной квадратичной формой, относительно потенциалов . Матрица, определяющая эту форму, получила название матрицы Дирихле. Поэтому в этом случае операция определения экстремума функционала превращается в обычное дифференцирование по n переменным. Подчеркнем, что дифференцирование осуществляется только по неизвестным потенциалам. Таким образом, схема метода конечных элементов заключается в последовательном выполнении следующих операций:

• разбиение области на треугольные элементы;
• определение матрицы Дирихле для каждого элемента и объединение ее на всей области;
• разделение потенциалов на известные и искомые;
• определение производной функционала и приравнивание ее к нулю;
• решение системы линейных уравнений и определение неизвестных потенциалов в узлах симплексов.

Оценка точности

Абсолютную точность, как указывалось ранее, дают аналитические решения, однако, в большинстве практических случаев получить эти решения не представляется возможным. Тем ни менее во многих практических приложениях используются конструкции, позволяющие найти эти аналитические решения. Например, коаксиальные кабели, волноводы, установки по созданию сверхмощного импульса тока РС-20, электростатические фильтры для отчистки технических жидкостей и газов [9, 10] и другие устройства.

Одновременно такого рода задачи, имеющие аналитические решения, можно использовать в качестве эталонных для определения точности численных алгоритмов [11].

Выберем следующую задачу [12, 13]: задан коаксиал с радиусом внутреннего проводника r₁ = 1 мм, внутренним радиусом внешнего трубчатого проводника r₂ = 4 мм. Потенциал внутреннего проводника ${\left.\text{φ}\right|}_{r_1}=400 В$, а трубчатый проводник заземлен ( ${\left.\text{φ}\right|}_{r_2}=0 В$ ).

В цилиндрической системе координат уравнение Лапласа для данной задачи имеет известный вид [13]. Решением данного уравнения является распределение потенциала $\text{φ}$ внутри коаксиала. Аналитическое представление потенциала выглядит:

$\text{φ}=C_1\ln r+C_2,$ (9)

где $C_1$ и $C_2$ – константы, определяемые из граничных условий.

Для простоты рассмотрения и учитывая осесимметричность ограничимся четвертью коаксиальной области.

Для оценки точности численного решения будем разбивать область между проводниками на различное количество симплексов. Для примера рассмотрим два варианта, представленные на Рис. 3.

 

Рис. 3 - Варианты разбиения на конечные элементы: а - 16 элементов, б - 200 элементов

 

Матрица Дирихле каждого элемента легко вычисляется [5], поскольку производные в (5) от линейной функции (8) являются константами, хотя выкладки носят громоздкий характер. Для примера приведен элемент этой матрицы:

$S_{12}^{\left(e\right)}=\frac{1}{4A}\left\{\left(y_2-y_3\right)\left(y_3-y_1\right)+\left(x_3-x_2\right)\left(x_1-x_3\right)\right\},$ (10)

где A – площадь треугольного элемента.

Остальные восемь получаются путем циклической перестановки индексов координат. Энергия одного элемента определяется на основании формул (1) и (5):

$Ф\left[\text{φ}\right]=\frac{1}{2}{\text{φ}}^T{S}^{\left(e\right)}\text{φ},$ (11)

где T – символ транспонирования вектора потенциалов ${\text{φ}}^T=\left[{\text{φ}}_1,{\text{φ}}_2,{\text{φ}}_3\right]$.

На примере (11) видно, что функционал для данного элемента имеет квадратичную форму трех переменных. Поэтому определение экстремума функционала по всей области $\text{Ω}$ сводится к дифференцированию энергии по свободным потенциалам и приравниванию результата к нулю. Полученная система линейных уравнений имеет симметричную ненулевую матрицу и потому имеет решение.

Проведем исследования погрешности вычислений в зависимости от количества конечных элементов, покрывающих область $\text{Ω}$. Будем использовать энергетическую метрику для оценки погрешности расчетов потенциалов

$\text{ε}=\sqrt{\underset{\text{Ω}}{\overset{}{}}{\left|{\text{φ}}_a\left(x,y\right)-{\text{φ}}_d\left(x,y\right)\right|}^{2}dxdy},$ (12)

где ${\text{φ}}_a$ – аналитическое решение (9);

${\text{φ}}_d$ – численное решение.

При этом погрешность вычисления энергии составляет

$\text{ξ}=\left|\text{Φ}\left({\text{φ}}_a\right)-\text{Φ}\left({\text{φ}}_d\right)\right|.$ (13)

Для удобства интерпретации погрешности выражены в относительных единицах в %. В качестве базовых значений для вычисления относительных погрешностей взята норма аналитического потенциала ${\phi }_a$ и соответствующая энергия.

Данные расчёты представлены на Рис. 4.

 

Рис. 4 - Результат вычисления относительных погрешностей

 

Произведем анализ результатов расчетов. Представим, что неточность вычислений потенциалов есть функция $h\left(x,y\right)$, норма которой $h$ меньше некоторой положительной малой величины. Причем функция ${\left.h\right|}_{\text{Γ}}$ на границе области $\text{Ω}$ имеет значение равное нулю в силу того, что в этих точках потенциалы заданы и не вычисляются.

Рассмотрим изменение энергии, вызванное неточностью вычислений

$\text{Φ}\left(\text{φ}+\text{θ}h\right)=\text{Φ}\left(\text{φ}\right)+\text{θ}\underset{\text{Ω}}{\overset{}{}}\left(\nabla \text{φ}\right)\left(\nabla h\right)dxdy+\frac{1}{2}{\text{θ}}^2\underset{\text{Ω}}{\overset{}{}}{\left(\nabla h\right)}^{2}dxdy,$ (14)

где $\text{θ}$ – малая положительная величина;

$\nabla \nabla$ – градиент.

Используя теорему Грина данное выражение можно привести к виду [5]:

$\text{Φ}\left(\text{φ}+\text{θ}h\right)=\text{Φ}\left(\text{φ}\right)+{\text{θ}}^2\text{Φ}\left(h\right).$ (15)

Анализ формулы (15) показывает, что теоретически ошибка в вычислении потенциала пропорциональна $\theta$. При этом ошибка в определении энергии пропорциональна ${\text{θ}}^2$. Однако, проведенный вычислительный эксперимент, показывает, что для данной задачи это выполняется в диапазоне от 72 до 200 элементов. При дальнейшем увеличении теория предсказывает, что ошибка по энергии должна уменьшаться существенно быстрее, чем ошибка по потенциалам, однако в действительности происходит все ровно наоборот. Причина такого поведения погрешностей может быть объяснена наличием вычислительных погрешностей, связанных с ограничением разрядности. Эти ограничения имеют случайный характер и поэтому проявляется неустойчивость ввиду уменьшения размеров конечных элементов. Отсюда следует вывод о том, что для дальнейшего уменьшения ошибки необходимо использовать на конечных элементах нелинейные поверхности более высоких порядков.

Величины погрешностей, при незначительном увеличении количества симплексов, лежат в пределах допустимых величин, характерных для метода конечных элементов [14].

Практическая значимость

Рассмотренный метод может быть использован для анализа внешних воздействий при решении задач электромагнитной совместимости. Множество подобных проблем существует в электротехническом комплексе железнодорожного транспорта [15]. Например, известна проблема электротермической деградации волоконно-оптических кабелей, подвешенных на опорах контактной сети. Одной из причин этого явления является образование на оболочке кабеля электропроводящей пленки из частиц пыли. В середине пролета между опорами потенциал поля максимальный, а крепления кабеля на опорах заземлены и имеют минимальный потенциал, это притягивает заряженные частицы к креплению и в итоге они оседают на оболочке кабеля вблизи опоры. Далее наводимый в этих пленках ток разогревает оболочку и способствует ее ускоренной деградации. Знание распределения электрического поля в данной задаче позволит усовершенствовать методы защиты волоконно-оптического кабеля, тем самым снизив число аварий в системах обеспечения безопасности движения поездов.

Заключение

Применение численных методов позволяет решать сложные задачи даже в случае отсутствия аналитического решения [1, 14, 16]. Наряду с этим применение метода конечных элементов дает дополнительные преимущества. К ним относится возможность расчета конструкций, либо задач с граничными условиями произвольной формы, например, как при расчете электромагнитных полей реальных объектов. Кроме того, простота перехода к неравномерным покрытиям расчетной области и высокая технологичность, масштабируемость и т.д., особенно проявляющиеся при алгоритмизации и программировании задач.

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

Авторы заявляют, что:

1. У них нет конфликта интересов;
2. Настоящая статья не содержит каких-либо исследований с участием людей в качестве объектов исследований.

About the authors

Vladimir N. Taran

Rostov State Transport University; Don State Technical University

Email: vladitaran@rambler.ru
ORCID iD: 0000-0002-1308-3883
SPIN-code: 7763-2292
Scopus Author ID: 57205870949

Professor, Doctor of Physics and Mathematics Sciences, Professor of the Department "Communication on the Railway Transport", Professor of the Department of Radioelectronics

Russian Federation, 344038, Rostov-on-Don, Rostovskogo Strelkovogo Polka Narodnogo Opolchenya sq., 2.; 344000, Rostov-on-Don, Gagarin Square 1

Maxim V. Shevlyugin

Russian University of Transport

Email: mx_sh@mail.ru
ORCID iD: 0000-0003-3300-5193
SPIN-code: 4123-6584
Scopus Author ID: 24922001300

Associate Professor, Doctor of Technical Sciences, Head of the Department of Electric Power Engineering of Transport, Professor of the Department of Electric Power Engineering of Transport

Russian Federation, 127994, Moscow, st. Obraztsova, 9, bldg. 9.

Aleksey V. Shandybin

Rostov State Transport University

Author for correspondence.
Email: shav850@mail.ru
ORCID iD: 0000-0003-2052-5480
SPIN-code: 4444-3213
Scopus Author ID: 57205099726

Head of the Laboratory of the Department «Communication on the Railway Transport»

Russian Federation, 344038, Rostov-on-Don, Rostovskogo Strelkovogo Polka Narodnogo Opolchenya sq., 2.

References

Supplementary files