РАСПРЕДЕЛЕНИЕ ПОТЕНЦИАЛА В ОКРЕСНОСТИ ОПОРЫ ЛИНИИ
ЭЛЕКТРОПЕРЕДАЧИ Пацюк В.И.
Институт Энергетики АНМ
Аннотация. В статье рассматривается новый подход применения метода конечных объемов для расчета электрического поля в пространственно-неоднородной трехмерной среде. Сформулирована задача Дирихле с построением расчетной сетки на основе разбиения пространства, известного как триангуляция Делоне с применением в расчетной схеме ячеек Вороного. Предлагается численный алгоритм для расчета потенциала и напряженности электрического поля в пространстве в окрестности опоры линии электропередачи. Разработан алгоритм и программа численного счета в среде приложения Matlab. Приведены результаты расчета распределения в пространстве потенциала и напряженности электрического поля.
Ключевые слова: Электрическое поле, неоднородная среда, потенциал,
напряженность электрического поля, численный метод.
REPARTITIE TRIDIMENSIONALA A CAMPULUI ELECTRIC iN VECINATATEA
PILONULUI ELECTRIC Patiuc V.I.
Institutul de Energetica al A§M Rezumat. in lucrare se examineaza o noua abordare a utilizarii metodei volumelor finite pentru a calcula campul electric in mediul tridimensional neomogen spatial. S-a formulat problema Dirichlet cu construirea retelei de divizare denumita trianghiularea Delaunay §i celulele Voronoi. S-a propus algoritmul de calcul numeric a potentialului §i intensitatii campului electric in spatiul in vecinatatea pilonului electric. A fost elaborat algoritmul §i programa de calcul in mediul aplicatiei Matlab. Sunt prezentate rezultatele calculelor repartitiei potentialului §i a intensitatii campului electric in spatiu.
Cuvinte-cheie: Camp electric, mediu neomogen, potential, intensitatea campului electric, metoda numerica.
POTENTIAL DISTRIBUTION IN VICINITY OF TRANSMISSION LINE SUPPORT
Patsiuk V.I.
Institute of Power Engineering of the Academy of Sciences of Moldova
Abstract. The paper examines a new approach to finite volume method which is used to calculate the electric field of spatially inhomogeneous three-dimensional environment. The Dirichlet problem is formulated and the computational grid is constructed on the base of space partition, which is known as Delaunay triangulation with the use of Voronoi cells. To calculate the potential and the electric field strength in the vicinity of the electric line support the numerical algorithm is proposed. We have developed the algorithm and the corresponding software using Matlab application. There are presented some numerical results of calculations of the distribution of potential and of the electric field strength.
Keyword: Electric field, inhomogeneous medium, potential, electrical field strength, numerical method.
1. Введение
Целью работы является расчет потенциала и напряженности электростатического поля в окрестности опоры высоковольтной линии электропередачи напряжения 500 кВ.
2. Постановка краевой задачи
Рассмотрим задачу определения трехмерного распределения потенциала u (x, у, z) электростатического поля в многосвязной области Q, в которой относительная диэлектрическая проницаемость s a (x, у) принимает кусочно-постоянное значение: sa =
1 в воздухе и sa = 6 внутри изолятора. На рис. 1 изображена область Q, которая представляет собой параллелепипед со сторонами длиной 2Lx, Ly, Lz, окружающий опору с подвешенными проводами. Нижние три провода находятся под напряжением
500 кВ, а верхние два имеют нулевой потенциал. Плоскость г = 0 является плоскостью поверхности земли с нулевым потенциалом. Потенциал опоры также равен нулю. Предполагается, что плоскости х = Ьх, х = -Ьх и г = Ьг расположены на достаточно большом расстоянии от опоры и поэтому на них также потенциал равен нулю. На плоскостях у = 0 и у = Ьу задаются условия симметрии.
-Ьх О Ьх х
Рис. 1. ВЛ 500 кВ стандартного исполнения с горизонтальной подвеской проводов Функция и(х, у, 2) внутри области О удовлетворяет уравнению Пуассона
ё1у(в а §гаёи) = -а( х, у, г)
(1)
где а(х, у, 2) - плотность распределения свободных зарядов. Если в области О таковые отсутствуют, то уравнение (1) превращается в уравнение Лапласа ё1у(в а§гаёи) = 0. На границе Г = дО области О считаются известными либо значения и (х, у, 2 )
и( х, у, г)| Г = ^( х, у, г),
(2)
либо производная по направлению внешней нормали равна нулю (условия симметрии)
- = 0. (3)
дп г
Вектор напряженности поля E определяется через потенциал и по формуле E =-gradи, а вектор электрического смещения D = saE. На границах раздела
разнородных сред выполняются условие непрерывности: [и] = 0 и [(D, n)] = 0, где квадратными скобками обозначена разность предельных значений слева и справа от границы раздела, а n - вектор нормали к ней.
3. Построение дискретной модели
Численное решение сформулированной задачи (1)-(3) в объемной области Q = Q + Г состоит из нескольких этапов:
A. Построение дискретной сетки в трехмерной области, окружающей опору.
B. Дискретизация уравнения Лапласа (1) и получение системы линейных
алгебраических уравнений.
C. Решение системы и получение распределения поля потенциала.
D. Расчет напряженности электрического поля на базе поля потенциала.
Из перечисленных этапов самым сложным является первый. Процесс получения сетки в трехмерном пространстве состоит из нескольких шагов:
404- 'О............О........о •••<);•
. о . . - о- • О - • О • О- О' О- О • О ■ О • О- О- О О
о- о -о- с о о о • с • о о- о- о • о • о- о- о- о
о . о- о - О О О’ О- О • О • О • О- • О • О • О • • О- • О • • • • О.................................................................................................о-
30}
о°°со о
о о о о о о оо оооороооо ооооо;с
°0<Р
о о о о о о
О о о о о о
а)
6)
Рис. 2. Сетка в окрестности опоры а) и триангуляция области б)
a) В двумерной области, перпендикулярной проводам, строится сетка с использованием фронтального алгоритма (advancing front method) [1, 2]. Рассматриваются две двумерные области: содержащая опору и без опоры (только с проводами). В случае области с опорой задача разбивается еще на две: прямоугольная область в непосредственной близости от опоры и область вдали от опоры без прямоугольника вокруг опоры. После построения сетки в этих двух областях производится объединение двух сеток с исключением общих узлов. На рис. 2 приведено изображение сетки на границе области расчета и триангуляция области, построенная путем применения фронтального алгоритма.
b) Финальная трехмерная сетка собирается из набора плоских двумерных сеток, которые вблизи опоры располагаются на достаточно близком расстоянии друг от друга, а с удалением от опоры расстояние между плоскостями увеличивается
c) Затем средствами Матлаба производится трехмерная триангуляция Делоне и строятся ячейки Вороного, которые используются в дальнейшем для дискретизации уравнения Лапласа.
Трехмерная триангуляция состоит из набора тетраэдров (пирамидок), множество которых обозначим как , где h - максимальная длина всех сторон пирамидок.
*
Введем также дуальную сетку Th , которая состоит из так называемых ячеек Вороного. Каждая ячейка Вороного окружает один из внутренних узлов разностной сетки. На рис. 3 приведены примеры трехмерной ячейки Вороного для базовых узлов с номерами 13 и 2012. Видно, что базовый узел связан с другими узлами сетки, которые называются соседними узлами, а ячейка Вороного представляет собой многогранник с зелеными гранями. Каждая грань ячейки ортогональна отрезку, соединяющему базовый узел с соседним узлом, а точка пересечения грани и отрезка расположена посередине отрезка.
б)
Рис. 3. Ячейки Вороного для базовых узлов с номерами 13 а) и 2012 б), имеющие
форму многогранника
*
Обозначим базовый узел через Ро, а ячейку Вороного через Кро. Вершины
*
ячейки Вороного Кр обозначим буквами Qi, которые являются центрами сфер,
описанных около тетраэдров, имеющих в качестве вершины точку Ро .
В качестве приближенного решения задачи (1)-(3) будем рассматривать кусочнолинейную функцию ик (х, у, 2), которая должна быть непрерывной в области О и линейной на каждом тетраэдре К е Т^. Такие функции называются линейными сплайнами. На множестве тетраэдров Тк функцию ик (х, у, г) можно задать следующим образом.
Пусть тетраэдр К = Р1Р2Р3Р4 (рис. 4) является некоторым элементом множества Т^ и Р(х, у, z) - произвольная точка этого элемента. В этом треугольнике для каждой
вершины вводятся функции формы N 1 (х, у), 1 = 1,4, которые удовлетворяют условиям: функции формы являются линейными и в вершинах тетраэдра принимают значения 1
Г1, 1 = к
или 0: Ni (Рк) = \ . Функции формы можно представить и в явном виде через
[0,1 к
координаты вершин треугольника:
Ыг (х, y, г) = ™,гх + ^2,гУ + Щз,гг + ™4, •
(4)
Здесь , щ2г., щ3і и щ4і являются компонентами векторов Щ,і = 1,4. Для определения векторов Щ нужно решить 4 системы уравнений АЩ і = /, і = 1,4. Коэффициенты матрицы А формируются из координат вершин р = Рі (х і, уі, г і ),і = 1,4 тетраэдра
А =
' х1 У1 г1 1І
Х2 У 2 г 2 1
х3 Уз гз 1
V х4 У 4 г 4 1
и /г = С/1,г , У3,г , /Лг1 /к,г ='
1, і = к 0, і Ф к
Используя функции формы, для каждого узла (внутреннего и граничного) сетки введем базисную функцию ф; (х, у, z), 1 = 1, 2, ..., п, п+1, ..., п1, где через пи п1
обозначены соответственно число внутренних и всех узлов сетки. Функция ф; (х, у, z) является кусочно-линейной, т.е. непрерывной и линейной на каждом тетраэдре со значением, равным единице в узле Р и равным нулю во всех остальных узлах. Тогда приближенное решение ик (х, у, z) может быть представлено в виде линейной комбинации базисных функций
'Ч
иЬ (^ ^ г) = 2иг% (^ У, г).
(5)
г=1
Нетрудно проверить, что коэффициенты и/ в (4) равны искомым значениям потенциала в узле р (х/, у{, zi): иь (х/, у/, zi) = и/.
Здесь следует отметить, что при решении задачи (1)-(3) методом конечных элементов применяется метод Галеркина, который состоит в следующем. Подставим (5) в уравнение (1) и запишем условие ортогональности полученного выражения к базисным функциям ф к (х, у, z) для внутренних узлов сетки
I^у(8а§гаёий )фкёУ = -|°фкЛУ, к = 1, п; (6)
О
| ^у(в а вгаё]Г и,.ф,. )ф кйУ = £ а ыи, = ~к; (7)
О /=1 /=1
~к =-\офкЛУ; аы = |&у(вавгаёф, )фkdV = -[ва^ьёф,.§гаёфкЛу.
О О О
Так как в граничных узлах значения решения известны, система (7) принимает
вид
2 ак(и1 = рк , рк = рк - 2 ак/Мч, к = 1, п. (8)
/=1 /=п+1
В отличие от метода конечных элементов в методе конечных объемов применяется обобщенный подход Галеркина, который заключается в том, что в условиях ортогональности (6) используются базисные функции у к (х, у, z)
пространства ^2 (О) = ¿2 (О), определяемого следующим образом. Введем новые базисные функции ук(х,у,z) для дуальной сетки Т* по правилу: функция ук(х,у,z) принимает постоянные значения равные единице в ячейке Вороного для внутреннего узла Рк и равна нулю в остальной области, т.е. функции имеют вид кусочнопостоянной функции. Тогда условие ортогональности (6) с функциями ук (х, у, z) принимает вид
I^у(8а§,Га&ик)укЛУ = -|^укЛУ, к = 1 п (9)
ОО
или в силу того, что функция ук (х, у, Z) отлична от нуля только в ячейке Вороного Крк
|ё1у(ва§гаёий)ЛУ = - |а ЛУ, (10)
где Кр*к - ячейка Вороного для узла Рк .
Поэтому для получения системы линейных алгебраических уравнений для неизвестных значений функции и* в узлах сетки по методу конечных объемов можно
К
К
поступить следующим образом. Рассмотрим в трехмерном пространстве с декартовыми координатами Охух уравнение Пуассона ё1у(в а §гаёи) =—о(х, у, I) и проинтегрируем
это уравнение по объему ячейки Кр .
Тогда получим формулу, совпадающую с (10)
|ё1у(ва§гаёи)аУ = - |о(х, у, г)ёУ. (11)
Кк Къ
К левой части (11) применим теорему о дивергенции
ди
|ё1у(ва§гаёи)ОУ = |ва§гаёи ОБ = |ва(§гаёи,п)ёБ = | ва —аБ, (12)
к 1 дк' дк' дк' д-
где дК* - полная поверхность многогранника К* ; п - внешняя нормаль к
поверхности дК*, а ди / дп - производная от функции и по этой нормали. В этом случае уравнение (11) принимает вид
ди ---1
дп
I Ва д-аБ = -Го(х,у, 1)йУ . (13)
* гп •>
Таким образом, решение задачи (1)-(3) по методу конечных объемов сводится к аппроксимации соотношения (13) для ячеек Вороного внутренних узлов разностной сетки. Аналогичная процедура присуща и методу конечных разностей для сетки с ячейками в форме параллелепипедов, поэтому метод конечных объемов может рассматриваться как обобщение метода конечных разностей для поблочной дискретизации с ячейками произвольной формы. По этой причине метод конечных объемов сохраняет все преимущества метода конечных разностей, а по сравнению с методом конечных элементов более простым оказывается сам алгоритм построения конечно-разностных соотношений и отпадает необходимость в построении локальной и глобальной матриц жесткости для формирования разрешающей системы уравнений типа (8).
Пусть в ячейке Вороного К* (рис. 3) для базового узла Р0 буквами р, I = 0,8
обозначены узлы сетки; Б1,1 = 1,12 - площади граней, ортогональных отрезкам Р0р ; М,1 = 1,12 - точки пересечения отрезка Р0Pi с гранью Б1. Тогда интеграл по
*
поверхности дКр0 из формулы (13) аппроксимируем следующим образом: г ди х2 Г ди ч и (р) — и (р0) 0
| В а ~аБ =Х]В а д-аБ = ^В а (М, ) --- Б
дК Р0 дп <=1 Б дп ¡=1
р р 1 01 i
где
р р 1 01 i
- длина отрезка р0 р1 .
Интеграл в правой части (12) аппроксимируем по формуле
Ja( x, y )dV = a( Po)Vo,
где У0 - объем ячейки Вороного Кр . Тогда аппроксимацию уравнения (13) можно
Po
представить в виде
J s „ (M,)"(P)—(Po) S, =-о( Po)Vo.
, =1
В окончательном виде уравнение для узла Po выглядит так
12
a o" (Po) + Za,"(P) = -^(Po)Vo; (14)
i=1
= ea(M,)ф=,, = 1,12; ao =-Ja,
" " i =1
P P
1 o1 i
Для каждого внутреннего узла сетки запишем уравнение вида (14), а для граничных узлов, в которых задано значение потенциала, используется условие (2). Граничные узлы, в которых задано условие симметрии, не включаются в систему уравнений (14) и соответствующие значения коэффициентов а, исключаются из
матрицы системы уравнений (14). В результате получаем систему линейных алгебраических уравнений с симметричной матрицей. Следует отметить, что при решении прикладных задач число уравнений системы исчисляется тысячами или десятками, сотнями тысяч. Однако поскольку каждое уравнение типа (13) содержит всего несколько ненулевых элементов (обычно от 9 до 25), то итоговая матрица получается достаточно разряженной. В реализованном алгоритме в памяти компьютера хранятся только ненулевые значения элементов матрицы. Для решения построенной системы используется итерационный метод сопряженных градиентов, который для задач рассматриваемого типа сходится очень быстро.
Полученное в области Q решение "h (x, y, z) позволяет построить поле потока
вектора напряженности E = (Ex, Ey, Ez) =-grad" . Обозначим через v поток вектора
E, проходящий через площадку единичной площади, ортогональной вектору E. Тогда изолинии "(x, y, z) = const и v(x, y, z) = const образуют взаимно ортогональные семейства. Функция v(x, y, z) может быть получена путем вычисления следующего контурного интеграла
( x,y, z)
v( x y) = j(Exdx + Eydy + Ezdz), (15)
( ^ Уo, zo)
где х0, у0, г0 - координаты произвольной фиксированной точки из области О, а путь интегрирования расположен внутри нее. В случае многосвязной области путь
K
интегрирования также не должен пересекать разрезы области, сводящей ее к односвязной структуре.
Емкость С между двумя проводящими телами вычисляется по формуле
С =
ч
и - и 2
(16)
где (и1 - и2) - разность потенциалов этих тел. Заряд ч тела, расположенного внутри некоторой трехмерной области У, вычисляется в соответствии с теоремой Гаусса как интеграл по поверхности Б = дУ от вектора напряженности Е
Ч = в| Е • = -в| §гаёи •
-в[(§гаёи • п)ё5 = -в [—.
•’ дп
(17)
Здесь Б - произвольная поверхность, содержащая в себе заряженное тело; п - вектор внешней нормали к ней; в - диэлектрическая проницаемость.
4. Практическая реализация
Описанный выше алгоритм был реализован в виде комплекса программ в среде приложения
МаИаЪ.
По разработанным алгоритму и программе была решена задача о трехмерном распределении потенциала и (х, у, г) в параллелепипеде Р = {(х, у, г): х е [-Ьх, Ьх ], у е [0, Ьу ], 2 е [0, Ьг ]] со
сторонами длиной Ьх = Ьг = 160 м, Ьу = 200 м. Размеры опоры представлены на рис. 5. Провода 3АС 400/51, ^0= 1,13 см, расстояние между составляющими - 40 см, радиус расщепления - 17,3 см. Эквивалентный радиус расщепленной фазы - 10,05 см.
Для оптимизации размера сетки шаги сетки варьировались в больших пределах. На большом расстоянии от опоры и провода шаг равнялся 4 м, с приближением к проводу и опоре шаг уменьшался и в окрестности провода равнялся 0,14 м. Общее количество узлов сетки равнялось 33024. Далее объем области был разбит на 1909507 тетраэдров.
Решение задачи представлено на рисунках 6. На рисунке 6 изображены линии равного потенциала (сплошные линии) и вектор напряженности (стрелки) для различных сечений вертикальными плоскостями при у = 36,23; 1,37; 0,675 и 0 м. Следует отметить, что результаты для у = 36,23 м уже практически не отличаются от результатов при у = 200 м, т.е. на расстоянии 36,23 м от опоры ее влияние на распределение поля потенциала не ощущается. Максимальное значение модуля вектора напряженности электрического поля достигается при у = 0 м в окрестности изоляторов и равняется 264,9 кВ/м. При удалении от опоры значение уменьшается и в сечении у = 36.23 м достигает величины 193,9 кВ/м.
По полученному распределению потенциала по формулам (16), (17) были посчитаны значения погонной емкости провода для различных расстояний от опоры.
у =152,4 м : С =
(7.27 1.80 0.61^1
1.80 7.15 1.80
0.61 1.80 7.27
пФ/м
у =0,675 м : С =
(10.32
0.88
0.25
0.88
10.39
0.88
0.25 ^
0.88
10.32
пФ/м.
г*
40
V.
•20й Ч*
%
200
.. . .
* ¿40 " '
•200 • •
•160' • * •КО ’ •
• «о
_*________*_____
*
АО
4
*
&
Во
*
о
»0
-30 -25 -20 -15 -10 -5 0
у = 36.23 м
5 10 15 20 25 30 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30
у = 1.37 м
Рис. 6. Линии равного потенциала (черные линии) и вектора напряженности (стрелки) для различных
сечений при у = 36,23; 1,37; 0,675 и 0 м
Выводы
1. Решена задача расчета численным методом трехмерного распределения потенциала u (х, у, z) и напряженности электростатического поля в многосвязной области Q, в
которой абсолютная диэлектрическая проницаемость s a (х, у) принимает кусочнопостоянное значение.
2. Предложен алгоритм решения задачи на основе метода конечных объемов. С целью уменьшения объема запоминаемой информации в реализованном алгоритме в памяти компьютера хранятся только ненулевые значения элементов матрицы.
3. Метод конечных объемов сохраняет все преимущества метода конечных разностей, а по сравнению с методом конечных элементов более простым оказывается сам алгоритм построения конечно-разностных соотношений. При этом отпадает необходимость в построении локальной и глобальной матриц жесткости для формирования разрешающей системы уравнений для емкостных коэффициентов.
4. Для решения построенной системы используется итерационный метод сопряженных градиентов, который для задач рассматриваемого типа сходится очень быстро.
Литература
[1] Немировский Ю.В., Пятаев С.Ф. Автоматизированная триангуляция многосвязных областей со сгущением и разрежением узлов. Вычислительные технологии, т. 5, N 2, 2000, с. 82-91
[2] Пашков С.В. Численное моделирование разрушения в динамических задачах
механики деформируемого твердого тела (МДТТ).
http://ps300.narod.ru/fr3d/mesh.htm
Сведения об авторе
ОПацюк В.И. Доктор физико-математических наук, ведущий научный сотрудник Института энергетики АНМ и доцент Государственого Университета Молдовы. Область научных интересов: численные методы расчета, математическая физика, механика твердого тела, теоретическая электротехника и энергетика. Автор более 100 научных публикаций, в том числе 1 патента и 10 монографий. E-mail: [email protected]