Метод конечных объемов для расчета электростатического поля трехмерной сферы с кубической полостью Текст научной статьи по специальности «Математика»

FINITE VOLUME METHOD FOR CALCULATION OF ELECTROSTATIC FIELD IN THREE-DIMENSIONAL SPHERE WITH CUBIC CAVITY

Patsiuk V.

Institute of Power Engineering of the Academy of Sciences of Moldova

Abstract. The paper presents the numerical approach based on finite volume method, which is used to calculate the electric field in three-dimensional environment. For this reason for formulated Dirichlet problem it is proposed the construction of computational grid based on space partition, which is known as Delaunay triangulation with the use of Voronoi cells. The numerical algorithm for calculating the potential, electric field strength and capacitance in the 3-d sphere with cubic cavity is proposed. Numerical results for potential distribution, electric field strength and capacitance in three-dimensional sphere with cubic cavity are represented. Keywords: Electric field, potential, capacitance, numerical method.

METODA VOLUMELOR FINITE PENTRU CALCULAREA CAMPULUI ELECTROSTATIC

PENTRU SFERA CU GOLUL CUBIC Patiuc V.

Institutul de Energetica al A§M Rezumat. in lucrare se examineaza utilizarea metodei volumelor finite pentru a calcula campul electric in mediul tridimensional. S-a formulat problema Dirichlet cu construirea retelei de divizare a spatiului, denumita trianghiularea Delaunay §i utilizarea in schema de calcul numeric a celulei Voronoi. S-a propus algoritmul de calcul numeric a potentialului, intensitatii campului electric §i capacitatii sferei cu cavitatie de forma cubica. Sunt prezentate rezultatele calculelor repartitiei potentialului, intensitatii campului electric §i capacitatii sferei cu cavitate de forma cubica.

Cuvinte-cheie: Camp electric, potential, intensitatea campului electric, capacitatea, metoda numerica.

МЕТОД КОНЕЧНЫХ ОБЪЕМОВ ДЛЯ РАСЧЕТА ЭЛЕКТРОСТАТИЧЕСКОГО ПОЛЯ ТРЕХМЕРНОЙ СФЕРЫ С КУБИЧЕСКОЙ ПОЛОСТЬЮ Пацюк В.И.

Институт Энергетики АНМ Аннотация. В статье рассматривается применения метода конечных объемов для расчета электрического поля в трехмерной среде. Сформулирована задача Дирихле с построением расчетной сетки на основе разбиения пространства, известного как триангуляция Делоне с применением в расчетной схеме ячеек Вороного. Предлагается численный алгоритм для расчета потенциала, потока вектора напряженности электрического поля и емкости в трехмерной сфере с кубической полостью. Приведены результаты расчета распределения в пространстве потенциала, потока вектора напряженности электрического поля и емкости сферы с кубической полостью.

Ключевые слова: Электрическое поле, потенциал, напряженность электрического поля, емкость, численный метод.

1. Введение

Многие явления и процессы окружающего нас мира могут быть описаны в виде начально-краевых задач для уравнений в частных производных. Разработанные математические модели должны адекватно соответствовать объекту и соответствовать требованиям решаемых инженерных задач. Электротехника и электроэнергетика являются относительно новыми научными дисциплинами, которые ассимилируют знания и из других областей науки с более длительной историей развития, например, из механики сплошной среды [1-3]. Математическое моделирование широко используется при изучении процессов в физике плазмы, синергетике и других разделах современного естествознания, где численное интегрирование многомерных эволюционных уравнений с учетом различного рода нелинейностей [4-6] является широко используемым подходом.

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

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

Использование уравнений Максвелла для изучения явлений электромагнитного взаимодействия в задачах электродинамики является достаточно сложной математической задачей из-за необходимости учесть взаимодействие электромагнитного поля с веществом в неоднородной среде. Сложность проблемы заключается в том, что вещество состоит из громадного количества частиц, движение которой в отдельности невозможно описать. Для того, чтобы обойти эту трудность используются некоторые модели среды. Отметим, что разработка и обоснование таких моделей не является тривиальной задачей, особенно для неоднородных сред.

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

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

2. Постановка краевой задачи

Рассмотрим задачу определения трехмерного распределения потенциала и(х, у, £) электростатического поля в многосвязной области О, в которой абсолютная диэлектрическая проницаемость еа (х, у) принимает кусочно-постоянное значение. Функция и(х, у, г) внутри области О удовлетворяет уравнению Пуассона

div(s а gradu) = -ст(х, у, г), (1)

где ст(х, у, г) - плотность распределения свободных зарядов. Если в области О таковые отсутствуют, то уравнение (1) превращается в уравнение Лапласа div(s а gradu) = 0. Граница Г = дО области О разбивается на 2 части: Ц и Г2, Ц пЦ = 0. На Ц значения и(х, у, г) считаются известными

и(х, у, ц = х, у, г) , ф)

а на Ц задаются условия симметрии в виде

ди дн

= 0, (2Ь)

г,

где н - вектор внешней нормали к Г2.

Вектор напряженности поля Е определяется через потенциал и по формуле Е = -grad и , а вектор электрического смещения Ъ = ваЕ. На границах раздела разнородных сред выполняются условие непрерывности: [и] = 0 и [(Ъ, п)] = 0, где квадратными скобками обозначена разность предельных значений слева и справа от границы раздела, а п - вектор нормали к ней.

3. Метод конечных объемов для уравнения Пуассона

Для численного решения сформулированной задачи Дирихле объемная область 0 = 0 + Г разбивается на конечное множество малых объемных элементов, имеющих форму тетраэдров (пирамидок). Вершины пирамидок называются узлами разностной сетки. На рис. 1 представлен пример такого разбиения объемного тела в форме сферы с полостью в виде куба. На рисунке представлена 1/8 часть сферы. При фиксированном положении узлов разностной сетки может быть построено большое количество различных разбиений трехмерной области на пирамидки. Наилучшим считается разбиение, которое называется триангуляцией Делоне. Это такое разбиение, при котором в сферу, описанную вокруг конкретной пирамиды, не попадают другие узлы разностной сетки.

Рис. 1. Сфера с кубической полостью, разбитая на 4501 тетраэдров на трехмерной сетке

с 955 узлами.

Множество пирамидок обозначим как Т *, где Н - максимальная длина всех сторон

*

пирамидок. Введем также дуальную сетку Т* , которая состоит из так называемых ячеек Вороного. Каждая ячейка Вороного окружает один из внутренних узлов разностной

сетки. На рис. 2 изображен пример трехмерной ячейки Вороного для базового узла с номером 74. Видно, что базовый узел связан с узлами с номерами 2, 73, 75, 80, 176, 182, 188 и 332, которые называются соседними узлами, а ячейка Вороного представляет собой многогранник с зелеными гранями. Каждая грань ячейки ортогональна отрезку, соединяющему базовый узел с соседним узлом, а точка пересечения грани и отрезка расположена посередине отрезка.

Обозначим базовый узел через р , а ячейку Вороного через Кр . Вершины ячейки Вороного К* обозначим буквами Qi, которые являются центрами сфер, описанных около тетраэдров, имеющих в качестве вершины точку Р0 .

Рис. 2. Ячейка Вороного для базового узла с номером 74, имеющая форму многогранника

В качестве приближенного решения задачи (1), (2) будем рассматривать кусочно-линейную функцию ий (х, у, 2), которая должна быть непрерывной в области О и линейной на каждом тетраэдре К е Тк. На множестве тетраэдров ТА функцию ий (х, у, 2) можно задать следующим образом.

Пусть тетраэдр К = рР2Р3Р4 (рис. 3) является некоторым элементом множества Т* и Р( х, у, г) - произвольная точка этого элемента. В этом тетраэдре для каждой вершины вводятся функции формы Ni (х,у, г), г = 1,4, которые удовлетворяют условиям: функции являются линейными и в вершинах тетраэдра принимают значения 1 или 0:

Г1,г = к

Ni (Р ) = ^ . Функции формы можно представить и в явном виде через координа-

[0,г * к

ты вершин тетраэдра:

N (х, у, г) = wг(1) х + у + wг'-3) г + w.

(2),

,(3),

,(4)

(3)

Здесь wг(1), wг(2).

w ,(3) и w (4)

являются компонентами векторов , г = 1,4 . Для определения векторов щ нужно решить 4 системы уравнений = ^, г = 1,4. Коэффициенты матрицы А формируются из координат вершин р = р (xi, уг, zi),г = 1,4 тетраэдра

А =

(х УХ 7х

х2 у 2 7 2 1

хз уз 73 1

V х4 у4 7 4 Х,

и I =(/,(1),(2),/,(3\(4)) ^)

(1, г = к

0, г * к.

Используя функции формы, для каждого узла (внутреннего и граничного) сетки введем базисную функцию ф г (х, у, г), г = 1, 2, ..., п, п+1, ..., щ, где через п обозначено число внутренних и граничных из Г узлов сетки, а через щ - число всех узлов сетки. Функция ф (х,у, г) является кусочно-линейной, т.е. непрерывной и линейной на каждом тетраэдре со значением, равным единице в узле р и равным нулю во всех остальных узлах. Тогда приближенное решение ий (х, у, г) может быть представлено в виде линейной комбинации базисных функций

"1

(X, У, 2) = Х и' ф1 (Х' У' 2)- (4)

1 =1

Нетрудно проверить, что коэффициенты и, в (4) равны искомым значениям потенциала в узле р (X' у,' и * (X' У1' гг) = и.

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

|&у(ва%гади*)фкйУ = -|стфкйУ, к = 1,п; (5)

□ □

п1 п1 ~ | в иг ф, Уфкйу = ^а к1иг = ~к; (6)

□ 1=1 1=1

~к =-\ъфкЛу; ак, = |вБгаёф,)фк^у = -|бвБгаёф,Бгаёфк^у.

□ □ □

Так как в граничных узлах границы ^ значения решения известны, система (6) принимает вид

п ^ " _

Еаши, = Рк, Рк = Рк- Еак,^> к=п. (7)

,=1 ,=п+1

В отличие от метода конечных элементов в методе конечных объемов применяется обобщенный подход Галеркина, который заключается в том, что в условиях ортогональности (5) используются базисные функции (х, у, 2) пространства

Ж2°(П) = ^(О), определяемого следующим образом. Введем новые базисные функции

*

(х, у, г) для дуальной сетки Т* по правилу: функция (х, у, 2) принимает постоянные значения равные единице в ячейке Вороного для внутреннего узла Рк и равна нулю в остальной области. Тогда условие ортогональности (5) с функциями (х, у, г) принимает вид

|ё1у(8аgгaduhкйУ = -|стукйУ, к = 1,п (8)

□ □

или в силу того, что функция (х, у, г) отлична от нуля только в ячейке Вороного Кр

|<у(8аБгаёи)йУ = - |айУ, (9)

где К* - ячейка Вороного для узла Рк.

Поэтому для получения системы линейных алгебраических уравнений для неизвестных значений функции и в узлах сетки по методу конечных объемов можно по-

и

К Р

К Р

ступить так. Рассмотрим в трехмерном пространстве с декартовыми координатами Охуг уравнение Пуассона ёгу(еа§гаём) = -ст(х, у, г) и проинтегрируем это уравнение по объему ячейки К* .

Тогда получим формулу, совпадающую с (9)

|йу(еаgгadu)ау = - /ст(х, у, г)ау. (10)

К р Ки

К левой части (10) применим теорему о дивергенции

|div(sagгadu)ау = /^§гаёи аь = |еа^га<3и,п)аь = / ва—аь, (11)

■с'л дК * "'* дп

ди

gгadu, п) аь = | ва

дКк

где дК* - полная поверхность многогранника К* ; п - внешняя нормаль к поверхно-

Рк

сти дК* , а ди / дп - производная от функции и по этой нормали. В этом случае урав

нение (10) принимает вид

| га—аь = - /а(х,у,г)а¥ . (12)

дК'к дп

Таким образом, решение задачи (1), (2) по методу конечных объемов сводится к аппроксимации соотношения (12) для ячеек Вороного внутренних узлов разностной сетки. Аналогичная процедура присуща и методу конечных разностей для сетки с прямоугольными ячейками, поэтому метод конечных объемов может рассматриваться как обобщение метода конечных разностей для поблочной дискретизации с ячейками произвольной формы. По этой причине метод конечных объемов сохраняет все преимущества метода конечных разностей, а по сравнению с методом конечных элементов более простым оказывается сам алгоритм построения конечно-разностных соотношений и отпадает необходимость в построении локальной и глобальной матриц жесткости для формирования разрешающей системы уравнений типа (7).

Пусть в ячейке Вороного К * (см. рис. 2) для базового узла Р0 буквами р, , = 0,8

обозначены узлы сетки; , , = 1,8 - площади граней, ортогональных отрезкам Р0р ; М{, , = 1,8 - точки пересечения отрезка рр с гранью £. Тогда интеграл по поверхности дК* из формулы (12) аппроксимируем следующим образом:

I ^а ^ = ± /в. ^ ,£ва (М, ) £

дп ' дп

Р Р р 0р

где

р р 1 0 р

- длина отрезка р р . Интеграл в правой части (12) аппроксимируем по формуле

/о( х, у )ау = о( РЖ,

дК *

к :

где У0 - объем ячейки Вороного Кр . Тогда аппроксимацию уравнения (12) можно представить в виде

£ в. (и,) иЩт з, =-« РоЖо.

,=1

Р Р р 0р

В окончательном виде уравнение для узла Ро выглядит так

а о и(Ро) + ги(Рг ) = "°(Р0)Ж0 ;

1 =1

(13)

3 __8

а, =8 в (М, , I = 1,8; а о =-£а,

Р Р

1 о Р

1=1

Для каждого внутреннего узла сетки запишем уравнение вида (13), а для граничных узлов используется условие (2а). В результате получаем систему линейных алгебраических уравнений с симметричной матрицей. Следует отметить, что при решении прикладных задач число уравнений системы исчисляется тысячами или десятками тысяч. Однако поскольку каждое уравнение типа (13) содержит всего несколько ненулевых элементов (обычно от 9 до 25), то итоговая матрица получается достаточно разряженной. В реализованном алгоритме в памяти компьютера хранятся только ненулевые значения элементов матрицы. Для решения построенной системы используется итерационный метод сопряженных градиентов, который для задач рассматриваемого типа сходится очень быстро.

4. Определение потока вектора напряженности

Вектор напряженности в двумерном случае вычисляется по формуле

(ди 1

Е =

( Е' 1 I Еу J

= -§гаё и = -

дх ди

ду

(14)

Поток вектора напряженности v(xy) определяется как интеграл по линии равного потенциала

(15)

у( х, у) = | (Е, пуН,

где п - вектор нормали к линии I.

При решении двумерной задачи потенциал и и поток вектора напряженности V связаны условиями Коши-Римана: ди / дх = дv / ду; ди / ду = -дv / дх [1о]. Из этого следует, что обе функции являются гармоническими и сопряженными, линии уровня этих функций взаимно ортогональны и функция V может быть получена путем вычисления следующего контурного интеграла от известной функции и

(х, у)

Ч х у) = |

( x0, уо)

д^ . дv , ^ — ах н--ау

^дх ду

(х,у)

' ( x0, уо)

ди , ди , ^

--ах н--ау

ду дх

(16)

где х0, уо - координаты произвольной фиксированной точки из области □ е Я2, а путь интегрирования расположен внутри нее. Так как выполняются условиями Коши-Римана, то значение интеграла не зависит от пути интегрирования. В случае многосвязной области путь интегрирования также не должен пересекать разрезы области, сводящей ее к односвязной структуре.

Полученное в трехмерной области Ое Я3 решение ик (х, у, г) позволяет вычислить

вектор напряженности Е = (Ех, Е , Бг) = ^гаёи. Однако задача построения потока

вектора напряженности в трехмерном пространстве не имеет единственного решения. Поэтому поступим следующим образом: выберем некоторую плоскость Р и найдем пересечение П2 = □ о Р . Область О2 является двумерным подмножеством области □ . По функции и (х, у, 2) методом интерполяции построим двумерную функцию и* л), л) е □ 2, а затем определим функцию -й л), градиент которой ортогонален градиенту и* (£, л) : (вгаё и* (£, л), вгаё V* (£, л)) = 0.

В трехмерном случае формула (16) будет давать неверные результаты, так как фун-ции и и V уже не удовлетворяют условиям Коши-Римана. Поэтому для построения на плоскости Е,Ол функции V, двумерный градиент §гаё2V = (д- / дv / дл) которой ортогонален вектору §гаё2 и = (ди / ди / дл), поступим следующим образом. Запишем условие ортогональности в виде равенства нулю скалярного произведения градиентов функций и и V:

дv ^С^и д- ^С^и

у,Вга<2и )= — — + — — = °. (17)

дс, дс, -л -л

Так как функция и является известной, то соотношение (1 7) является гиперболическим уравнением типа уравнения переноса для определения неизвестной функции -.

5. Вычисление емкости

Емкость С между двумя проводящими телами вычисляется по формуле

С = ® , (18)

и и 2

где (и - и2) - разность потенциалов этих тел. Заряд Q тела, расположенного внутри некоторой трехмерной области У, вычисляется в соответствии с теоремой Гаусса как интеграл по поверхности £ = СУ от вектора напряженности Е

Q = 8^ Е • й£ = -8^ §га<м • й£

£ £

-8| (§гас1и • Я)й£ = -8|■дUd£. (19)

Здесь £ - произвольная поверхность, содержащая в себе заряженное тело; п - вектор внешней нормали к ней; 8 - диэлектрическая проницаемость.

6. Практическая реализация

Описанный выше алгоритм был реализован в виде комплекса программ в среде приложения МайаЬ. Комплекс программ включает в себя блок построения трехмерной

£

£

сетки, блок решения уравнения Лапласа для определения потенциала, блок построения потока вектора напряженности и блок вычисления емкости.

В первом блоке формируется трехмерная сетка. Вначале формируются двумерные поверхности, форма которых плавно меняется от сферической вблизи внешней границы шара до кубической поверхности вблизи внутренней полости. Затем на каждой из этих поверхностей формируется треугольная сетка по алгоритму, описанному в [11] (рис. 4). Трехмерные координаты всех узлов запоминаются в едином массиве X и затем для построения трехмерной триангуляции Делоне используется Matlab-функция delaunayn(X).

Во втором блоке для определения потенциала решается уравнение Лапласа с граничными условиями Дирихле на внутренней и внешней поверхностях и с условиями симметрии на остальных границах. Задача решается методом конечных объемов, описанным в пункте 3. В результате требуется решить систему уравнений вида (13), порядок которой может составлять сотни тысяч. Матрица этой системы является сильно разряженной: с учетом симметрии максимальное количество ненулевых элементов составляет обычно 8-10. Поэтому не имеет смысла хранить в памяти компьютера всю матрицу системы. Для хранения сильно разряженной матрицы используется CSIR-схема [12], а система уравнений решается методом сопряженных градиентов.

Рис. 4. Триангуляция двумерных поверхностей

В третьем блоке вычисляется поток вектора напряженности. Для этого, как было описано в пункте 4, выбирается плоскость сечения Р вида х tgф — у = 0 при заданном

угле ф. Получаемая при этом двумерная область 02 = Р показана на рис. 5 и 6. Алгоритм определения функции V из гиперболического уравнения (17) состоит в следующем. На внешней границе области, имеющей форму четверти окружности, строим сетку и задаем начальные значения функции V равномерно возрастающие от 0 до 1. Затем строим второй слой узлов как вершины треугольников (рис. 5а). В каждом из треугольников значения функции и известны во всех трех узлах, а значения функции V в двух узлах известны из начальных данных, а для определения значения V в вершине треугольника поступаем следующим образом. Обозначим номера вершин треугольника через 1, 2 и 3. Будем предполагать, что функции и и V в треугольнике являются линейными. Тогда они могут быть представлены в виде: и(£, л) = А£ + Бц + С и v(£, л) = а£ + Ьц + с. Коэффициенты определяем из систем

^ ц 1Л( А Л (щ Л (£ ц 1У а Л

£ 2

^ 3

Л1 Л 2 Лз

в

/V С /

V из у

£ 2

^3

Л1 Л 2 Лз

Ь

/V с у

V ^ /

(20)

Тогда система (17) в каждом треугольнике может быть записана в виде

ду ди ду ди . , „ --+--= аА + ЬВ = 0

д£ д£ ^ дЛ

или

А(Лз —Л 2) + В(£ 2 —£ 3) А(Л1 —Л 2) + В(£ 2 —£1)

а

= 1 — а.

(21)

(22) (23)

Уравнение (22) обеспечивает устойчивость вычислений при условии 0 < а - 1. Ес-

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

(£з. Лз)

т.е.

положение вершины треугольника.

и

V

2

2

v3 = а V +а v2

а =

Рис. 5. Сетка и триангуляция для определения потока V

После того, как по формулам (22) определены значения V во всех внутренних узлах второго слоя, в граничных узлах значения доопределяются соответственно единицей и нулем. Затем строится следующий слой узлов и треугольников (рис. 5Ь) и по формулам, аналогичным (22), вычисляются значения V. Процесс повторяется до достижения внутренней границы (рис. 5^ d, e). Алгоритм формирования следующего слоя узлов предусматривает плавное изменение формы получаемой кривой от круговой до прямоугольной.

Результаты расчета полей потенциала и вектора напряжений в двух сечениях представлены на рис. 6.

Рис. 6. Линии равного потенциала (сплошные линии) и потока вектора напряженности (штриховые линии) для сечения плоскостью х = 0 и х - у = 0

В четвертом блоке комплекса программ по формулам (18), (19) вычисляется емкость сферы с кубической полостью. Для сферы диаметра 20 м и кубической полостью со стороной 10 м емкость равняется С = 2,20 нФ. Для сравнения приведем значения емкостей сферы диаметра 20 м со сферическими полостями диаметрами, вычисленные по точной формуле:

к, — г '

где - радиус сферы, г3 - радиус полости, в - абсолютная диэлектрическая проницаемость.

Диаметр сферы, м Диаметр сферической полости, м Емкость, нФ

20 10 1,113

20 14,14 2,686

20 17,32 7,192

Выводы

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

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

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

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

Литература

[1] Пацюк В.И., Римский В.К. Волновые явления в неоднородных средах. Т.1. Теория распространения упругих и неупругих волн. - Кишинев: Издательско-полиграфический центр МолдГУ, 2005. - 254 с.

[2] Пацюк В.И., Римский В.К. Волновые явления в неоднородных средах. Т.2. Нестационарное деформирование многосвязных тел. - Кишинев: Издательско-полиграфический центр МолдГУ, 2005. - 239 с.

[3] Самарский А.А., Михайлов А.П. Математическое моделирование: идеи, методы, примеры. - М.: Наука, 1997. - 320 с.

[4] Тихонов А.Н., Самарский А.А. Уравнения математической физики. - М.: Наука, 1977. -736 с.

[5] Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. - М: Наука, 1978. - 592 с.

[6] Самарский А.А., Попов Ю.П. Разностные схемы газовой динамики. - М.: Наука, 1975. -352 с.

[7] Li R., Chen Z. and Wu W. Generalized difference methods for differential equations. Numerical analysis of finite volume methods. New York-Basel: Marcel Dekker, Inc., 2000. - 459 p.

[8] Римский В.К., Берзан В.П., Пацюк В.И. и др. Волновые явления в неоднородных структурах. Т.5. Теория и методы расчета электрических цепей, электромагнитных полей и защитных оболочек АЭС. - Кишинев: Типография АНМ, 2008. - 664 с.

[9] Patsiuk V. Mathematical methods for electrical circuits and fields calculation. Chi§inau, Central Editorial-Poligrafic al USM, 2009. 442 p.

[10] Демирчян К.С., Нейман Л.Р., Коровкин Н.В., Чечурин В.Л. Теоретические основы электротехники. Том 3: Учебник для вузов. 4-е издание. С.-Пб: Питер, 2004, 384 с.

[11] Пашков С.В. Фронтальные алгоритмы построения треугольной сетки. http://ps300.narod.ru/fr3d/mesh.htm

[12] Баландин М.Ю., Шурина Н.П. Методы решения СЛАУ большой размерности. -Новосибирск: Изд-во НГТУ. - 70 стр.

Сведения об авторе.

Пацюк Владимир Ильич. Доктор физ.-мат.наук, доцент Государственного Университета Молдовы, ведущий научный сотрудник Института Энергетики АН Молдовы. Области научных интересов: математи-ческая физика, численный анализ, механика сплошных сред, теоретическая электротехника, силовая электроника. Автор более 140 научных публикаций, в том числе 9 монографий, одного изобретения.