CC BY
14ISSN 0868-5886
НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2006, том 16, № 1, c. 85-93
ОРИГИНАЛЬНЫЕ СТАТЬИ
УДК 517.925.4: 62 - 408.64
© С. И. Шевченко
О РЕШЕНИИ ДВУМЕРНОГО УРАВНЕНИЯ ПУАССОНА, КОГДА ОБЛАКО ОБЪЕМНОГО ЗАРЯДА ИМЕЕТ ГЛАДКИЕ КРИВОЛИНЕЙНЫЕ ГРАНИЦЫ. I. ПЛОСКАЯ ГЕОМЕТРИЯ
Приведен алгоритм решения двумерного уравнения Пуассона в плоской геометрии методом граничных интегральных уравнений на случай, когда облако объемного заряда имеет гладкие криволинейные границы.
ВВЕДЕНИЕ
В представленной работе, являющейся продолжением работ [1, 2], рассматривается развитие метода граничных интегральных уравнений при решении двумерного уравнения Пуассона (вычислении электростатического поля) на случай, когда облако объемного заряда имеет гладкие криволинейные границы, не сводимые к ступенчатым функциям.
Для вычисления потенциала электростатического поля (ЭСП) будем пользоваться интегральным уравнением, эквивалентным уравнению Пуассона [3]:
0(го) =\ао йь + -—|до ds, (1)
(Ь) 4Л£о (S)
где —0 = 8.8542 -10-12 Ф/м — электрическая постоянная; д — функция плотности объемного заряда (ПОЗ); <г(г) — функция плотности поверхностного заряда (ППЗ); 0(г, г0) — ядро интегрального уравнения (ядро для задачи Дирихле для уравнения Лапласа); г0— точка, в которой ищется потенциал (точка наблюдения (ТН)); первый интеграл берется по длине контура всех электродов (границы), а второй — по площади, которую занимает объемный заряд (ОЗ); коэффициент 1/ 4п—0 , который должен быть перед первым членом в правой части (1), включен в функцию а.
Ядро для задачи Дирихле для уравнения Лапласа в плоской геометрии имеет вид [3]
0(г, го) = -1п |г -г„| =
= -0.51п|(х -о)2 + (у - у о)2|, (2)
где г0 = (г0, г0), г = (г, г)— координаты точки интегрирования (ТИ).
Подробности решения уравнения (1) можно найти в [1, 2, 4, 5]. Отличие решения уравнения Пуассона от хорошо разработанного решения уравнения Лапласа [5, 6] заключается в необходимости вычисления как на первом этапе решения, так и на втором этапе решения второго интеграла в правой части (1), представляющего собой вклад в потенциал или в компоненту напряженности (КН) ЭСП от всего облака ОЗ. Этот интеграл, в котором д и О считаются известными, можно находить как численным методом (например, повторным интегрированием методом Гаусса), так и аналитически.
Для удобства вычислений область облака ОЗ покрывают некоторой прямоугольной сеткой, причем возможна как однородная так и неоднородная сетка. Все ячейки сетки можно отнести к одной из трех групп: внешние ячейки, полностью находящиеся вне облака ОЗ, внутренние ячейки, полностью находящиеся внутри облака ОЗ, и ячейки, которые пересекаются границей облака ОЗ и которые поэтому находятся частично вне а частично внутри облака ОЗ. В последнем случае часть площади полной ячейки, заполненную ОЗ, будем называть граничной ячейкой.
Методы вычисления двумерного интеграла имеют наиболее простой вид, если область интегрирования является прямоугольником (см. [1, 2]) или когда границу облака ОЗ можно представить ступенчатой функцией (линией), с тем чтобы можно было далее разбить область, занимаемую ОЗ, на совокупность прямоугольных областей.
Несколько другое дело, когда облако пространственного заряда имеет гладкую криволинейную границу. В этом случае границу облака ОЗ уже нельзя заменить ступенчатой функцией, особенно если облако ОЗ имеет резкую границу (четкую, не диффузную размытую границу, при которой плотность ОЗ плавно спадает к границе до нуля). Наличие криволинейной границы облака ОЗ застав-
ляет вводить в алгоритм решения уравнения Пуассона некоторые особенности, которые состоят в вычислении суммы вкладов в потенциал (или в КН ЭСП) от всех граничных ячеек и которые пока не были отражены в научной литературе.
Переход в представленном в [1, 2] алгоритме решения уравнения Пуассона от прямоугольных ячеек сетки к треугольным ячейкам может быть применен для повышения как точности вообще, так и точности вычисления в некоторой области, в которой, например, плотность ОЗ имеет некоторые особенности.*
Представленная работа посвящена разработке алгоритма вычисления вклада в потенциал (и/или КН ЭСП) от граничных ячеек облака ОЗ. При этом возникает ряд вопросов:
по каким критериям следует относить ячейки к одной из трех упомянутых выше групп;
как определять, принадлежит ли ТН к рассматриваемой граничной ячейке;
как находить вклад от ОЗ граничной ячейки?
1. АНАЛИЗ И ПОЗИЦИОНИРОВАНИЕ ЯЧЕЕК СЕТКИ, ИНТЕРПОЛЯЦИЯ
Будем полагать, что каждую ячейку сетки гладкая граница ОЗ может пересекать только один раз. (В противном случае следует сделать сетку более мелкой). Если граница облака ОЗ в рамках некоторой ячейки не является гладкой (имеет излом), то следует в данном месте так сдвинуть сетку, чтобы угол границы ОЗ пришелся на линию сетки, разделяющую две ячейки. В результате этого можно считать, что граница ОЗ в рамках каждой ячейки является гладкой функцией. В данной работе в качестве базового случая эта гладкая линия представлена прямой линией. Рассмотрено также обобщение на случай более сложной линии границы ОЗ в пределах одной ячейки.
1.1. Анализ
Вопросы геометрии пересечения граничных ячеек сетки с границей, отнесения каждой граничной ячейки к одному из установленных типов пересечения и относительного положения ТН и рассматриваемой граничной ячейки подробно рассмотрены в [7]. Показано, что всего возможны 12 типов пересечения ячеек сетки границей расчетной области (рис.). Выписан алгоритм анализа полученных вариантов пересечений ячеек сетки границей расчетной области, позволяющий для каждой ячейки установить реализуемый тип пересечения. Отличие данной работы состоит в том, что
* Частные беседы и частная переписка с Ивановым В.Я., 1978-2000 г.
¡IIIIIIIIIIIIIIIIIIIIJIII! Я я
1 2 3 4
111 f
А lililí 1 К
9 10 11 12
Типы пересечения ячейки сетки прямой линией. Заштрихованная область — граничная ячейка типа (1-12). Воспроизводится из [7]
вместо границы расчетной области используется граница облака ОЗ.
1.2. Позиционирование
При рассмотрении второго интеграла в правой части (1) важно знать, где расположена ТН (как она позиционируется в пределах граничной ячейки). Возможны следующие варианты: ТН лежит вне, внутри, точно на границе ячейки сетки. Возможность реализации этих вариантов зависит от того, на каком этапе решения уравнения Пуассона мы работаем:
На первом этапе решения уравнения (формирование матрицы, нахождение функции плотности поверхностного заряда) все ТН представляют собой точки коллокации, расположенные на границе расчетной области (электродах), а облако ОЗ расположено чаще всего вдали от электродов или более редко — касаясь некоторых частей электродов. Ясно, что ТН находятся или вне и далеко от ячеек, или в редких случаях расположены в пределах граничных ячеек (на той стороне граничной ячейки, которая является как границей облака ОЗ, так и поверхностью одного из электродов).
На втором этапе решения уравнения Пуассона (расчет потенциала или КН ЭСП) возможны самые различные варианты того, что считать и где считать. Однако одними из самых часто встречающихся вариантов являются:
1) расчет потенциала электростатического поля на сетке для последующего построения эквипо-тенциалей;
2) расчет КН ЭСП на сетке для последующего вычисления траекторий;
3) расчет вклада в КН ЭСП во внутренних точках граничных ячеек в процессе прямого расчета значений КН ЭСП для непосредственного их использования в программе расчета траекторий.
В ряде этих случаев значения КН ЭСП рассчитываются в узлах сетки и в точках пересечения линий сетки с граничной линией расчетной области (как это описано, например, в [7]).
1.3. Интерполяция по площади граничной ячейки
Для проведения интегрирования по площади граничной ячейки необходимо знать функцию плотности ОЗ. Вместо точной функции, которая практически никогда не известна, в формулу (1) обычно подставляется ее интерполяция. Мы ограничимся в качестве базового случая простейшей интерполяцией билинейной функцией, о которой можно найти информацию в справочниках [8, 9].
При проведении интерполяции функции плотности ОЗ по площади граничной ячейки большую роль играет форма этой ячейки (треугольная, четырехугольная, пятиугольная). Интерполяция для каждой из этих форм рассмотрена в [7].
Для треугольной формы:
Q(х, У) = co,o + Ci,ox + Co,i У
(3)
ячейки (рассмотрим граничную ячейку с конфигурацией N = 3) определяется выражением (второй член в правой части (1))
Фр =
Ук(х )
4п—
i K 7К\Л>
— Jdх J dy • Q(x ;y )• G(x,y;xa,y0), (6)
Ум(x )
Для четырехугольной формы:
Q(x, y) = Co,o + Ci o x + Co,i y + Ci,i xy . (4)
В случае пятиугольной формы граничной ячейки функция, интерполирующая ПОЗ, имеет пять точек (шаблона), в которых ПОЗ определена. Поэтому в данной работе пятиугольная граничная ячейка разбивается на треугольную и четырехугольную так, чтобы обе они опирались своими вершинами на точки, в которых определены значения плотности ОЗ. Для каждой их этих частей можно интерполяцию функции плотности ОЗ проводить независимо (т. е. по формуле (3) и формуле (4)).
Для повышения точности вычислений для всех форм ячеек (треугольной, четырехугольной и пятиугольной) может быть рассмотрена биквадра-тичная интерполяционная функция
Q(x, y) = Co o + Ci,o x + Coi у + cUi xy +
+ C2,ox2 + c2Jx2y + Co,2 y2 + Cu xy2 + c2 2x2y2. (5)
Для нахождения коэффициентов этой интерполяционной функции необходимо использовать не только базовые точки граничной ячейки (узлы сетки плюс точки пересечения границы облака ОЗ с линиями сетки), но и некоторые дополнительные точки, лежащие вне рассматриваемой ячейки. В качестве таких дополнительных точек использовались узлы сетки, ближайшие к рассматриваемой ячейке и лежащие внутри облака ОЗ.
2. АНАЛИТИЧЕСКОЕ ИНТЕГРИРОВАНИЕ ПО ПЛОЩАДИ ГРАНИЧНОЙ ЯЧЕЙКИ
Вклад в потенциал от произвольной граничной
где пределы внутреннего интегрирования: yN (х)— нижний предел и ук (х)— верхний предел, и пределы "внешнего" интегрирования: xN — нижний предел и хк — верхний предел определяются конфигурацией граничной ячейки, по площади которой проводится интегрирование.
Вычисление вклада в потенциал от граничной ячейки, осуществляемое аналитическим интегрированием, проводимым в формуле (6), по переменным (х, у), весьма схоже с его вычислением, проведенным в [1, 2]. Этот вклад имеет вид
где
. =
Ф =--!-р 8—
Yk (x)
(
É C,j ^,
i, j=o
(7)
= / хт а х / у" а у • 1п [(х - х0)2 + (у - Уо)2 ]. (8)
X« Г« (х)
Отличие заключается в том, что в интегралах 1тп данной работы один из пределов "внутреннего" интегрирования (по переменной у) является зависящим от "внешней" переменной (х). При этом "внешнее" интегрирование становится более громоздким, однако непреодолимых трудностей нет: все двукратные интегралы сводятся к совокупности табличных интегралов. Выражения для всех интегралов 1т п приведены в Приложении для случая, когда переменным является верхний предел "внутреннего" интегрирования (ук = ук (х) = = «• х + в ).
Порядок (последовательность переменных) интегрирования и то, какой предел "внутреннего" интегрирования (нижний или верхний) является переменным, определяются конфигурацией граничной ячейки. Для каждого из интегралов (8) нет необходимости рассматривать каждую из конфигураций граничной ячейки. Рассмотрим ячейку конфигурации номер 3, для которой интеграл (8) записывается (в системе координат, привязанной к ТН) в виде
ax+в
I,(3) =
J x'd x J yj d y • ln [ x2 + y2 ] .
x
Обозначим
этот
интеграл
через
3 (', 3, Уы,а, в) • Если обратиться к ячейке с конфигурацией N = 1, для которой этот интеграл имеет вид
хк Ук
/(3 = | х' а х | у3 а У • 1п [х2 + У2 ],
хы ах+в
и провести замену переменной интегрирования у —^ —у, то в результате получается
/£ = (—1)3 • 3(3)(',3, — Ук, —а, —в).
Т. е. удается свести к интегрированию по площади граничной ячейки, имеющей 3-ю конфигурацию.
Аналогично:
/® = (—1)' • 3(3)(3,— хк, —а, —в), /^ = 3 (3)( 3,хы ,а, в).
Для треугольных граничных областей все вполне аналогично, т. к. каждая из граничных треугольных областей может быть рассмотрена как предельный случай соответствующей четырехугольной граничной области. Пятиугольные граничные области мы просто разбиваем на треугольную и четырехугольную области, интегрирование по каждой из которых сводится к интегрированию по области с N = 3.
Таким образом, при вычислении интегралов по площади граничной ячейки достаточно рассмотреть всего одну конфигурацию граничной ячейки.
з. ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ МЕТОДОМ
ГАУССА
Метод одномерного численного интегрирования Гаусса характеризуется тем, что в нем точки интегрирования никогда не приближаются вплотную к пределам (нижнему и верхнему) интегрирования [10]. Поэтому в случае, когда ТН находится в одном из узлов сетки, в одной из точек пересечения линии сетки и граничной линии или лежит на граничной линии, ТИ не может совпасть с ТН,
и, следовательно, в этих случаях вопросы сингулярности ядра можно не рассматривать.
Другое дело, когда ТН лежит внутри граничной ячейки, где возможно прямое (или чрезвычайно близкое) совпадение ТН и ТИ, что может дать сингулярный (или неоправданно большой) вклад в интегральную сумму.
В качестве выхода можно предложить следующую простую процедуру: делим рассматриваемую граничную ячейку с помощью прямой линии, проведенной через ТН и параллельной оси Х (или оси
У), на две части. В каждой из полученных частей ТН находится точно на границе и поэтому при численном интегрировании по площади каждой из этих частей по методу Гаусса не может ни совпасть, ни очень сильно приблизиться к ТИ, что избавляет алгоритм интегрирования от проблемы сингулярности ядра.
Преимущество численного интегрирования в формулах для вклада в потенциал от облака пространственного заряда состоит в том, что имеется возможность использовать точные (а не интерполированные) значения функции плотности ОЗ в точках-узлах численного интегрирования. Кроме того, для формулы границы облака ОЗ в пределах ячейки, по площади которой проводится численное интегрирование, можно использовать как аналитические выражения, так и интерполяцию высокого порядка.
Конкретная реализация численного интегрирования зависит от формы занимаемой ОЗ части ячейки (треугольная, четырехугольная и пятиугольная).
3.1. Треугольная или четырехугольная граничная ячейка
Рассмотрим граничную ячейку типа 9 (или 3). Уравнение граничной линии, пересекающей ячейку, рассмотрено в работе [7] для всех двенадцати типов пересечения в виде Ук (х) = ах + в (для граничной ячейки типа 9 или 3). Оно определяет переменный верхний предел "внутреннего" интегрирования.
Выражение для вклада в потенциал от треугольной или четырехугольной формы граничной ячейки, реализующее повторное численное интегрирование методом Гаусса [10]:
Ф =——
Р 4П£
Ах х
ы1х
хХ 4
'=1
ёУ
АУ ( х )Х АУт б( х', Ут ) ° ( х, Ут \ Xo, Уо)
(9)
где Ах = Рх['х +1] — Р8хК['У ] — для треугольной граничной ячейки или Ах = Рх['х + 1] — Рх['х ] — для четырехугольной граничной ячейки; Рх['х ] — линия сетки, параллельная оси X и имеющая номер 1х; Р8хК['У] — х-координата точки пересечения линии сетки Ру['У ], параллельной оси У, и линии границы облака пространственного заряда;
АУ (х) = У к (х)— ру['У ]; ы8х, ыгУ— число узлов интегрирования вдоль направлений X и У; Ах1, АУт — коэффициенты квадратурной формулы
Гаусса; х1 = Рх[/х ] + Ах • хЦ) — координаты узлов интегрирования вдоль направления X; ут = ту ] + Ау(х1) • х«® (т) — координаты узлов
интегрирования вдоль направления У; х«х Ц) и
х(т) — корни каждого узла интегрирования;
xi, у^) — значения плотности объемного заряда в узлах интегрирования (х1, у.).
Полный вклад в потенциал (или КН ЭСП) от всего облака ОЗ получается суммированием вкладов от всех ячеек (внутренних и граничных), содержащих ОЗ.
3.2. Пятиугольная граничная ячейка
Особенностью пятиугольных граничных ячеек является то, что в каком бы порядке переменных не проводилось повторной интегрирование, по одной из переменных в качестве одного из пределов получается немонотонная линия, имеющая излом. Чтобы избежать этой трудности, пятиугольная граничная ячейка разбивается на две: треугольник и трапецию, как описано в разделе интерполяции. После чего и для трапециевидной, и для треугольной частей граничной ячейки интегрирование проводится согласно п. 3.1.
Если в формулы численного интегрирования подставлять вместо функции плотности ОЗ ее интерполяцию (4) (или(3)), то для вклада в потенциал приходим к выражениям (7, 8), в которых интегралы (8) вычисляются повторным численным интегрированием. Это дает возможность сравнивать (проверять) аналитическое и численное вычисления интегралов.
4. СМЕШАННОЕ (АНАЛИТИЧЕСКИ-ЧИСЛЕННОЕ) ИНТЕГРИРОВАНИЕ
Этот тип интегрирования состоит в том, что первая фаза интегрирования ("внутреннее" интегрирование) проводится аналитически (формулы записаны в системе координат, привязанной к ТН):
ук(х)
Ко(х) = / у0ау 1п[х2 + у2],
у«(х)
ук(х)
К1(х) = / у а у 1п [х2 + у2 ] ,
уN ( х )
а вторая фаза интегрирования ("внешнее" интегрирование) осуществляется численным интегрированием по методу Гаусса.
Понятно, что на второй фазе интегрирования возможна ситуация близкого расположения ТН и ТИ, в результате чего может получиться неоп-
равданно большой вклад в потенциал. Чтобы избежать этого для случая, когда ТН принадлежит ячейке, по площади которой проводится интегрирование, используется описанный выше метод с разбиением отрезка интегрирования на два и отдельным интегрированием по каждому из этих отрезков.
Приведенные ниже результаты тестов показывают, что смешанное интегрирование дает как для базовых интегралов (8), так и для вклада в потенциал от всего облака ОЗ результаты, практически совпадающие по точности с результатами аналитического интегрирования.
Еще одним положительным качеством этого типа интегрирования (или преимуществом по сравнению с аналитическим интегрированием) является то, что форма границы облака ОЗ во "внутреннем" интегрировании участвует, только когда подставляются пределы в найденную первообразную. Это позволяет легко использовать для описания границы облака ОЗ как аналитические выражения, так и интерполяцию высокого порядка.
К недостаткам смешанного интегрирования можно отнести то, что для его реализации обязательно следует провести интерполяцию функции плотности ОЗ в пределах каждой ячейки сетки, в которой присутствует ОЗ.
5. ОБСУЖДЕНИЕ
Для демонстрации свойств описанных выше алгоритмов, приведем некоторые результаты, появляющиеся при решении уравнения Пуассона для сферического конденсатора. Рассмотрим первый квадрант сечения сферического конденсатора плоскостью, по краям которого (на двух радиусах) поставим граничные условия, представляющие собой теоретические значения потенциала. Выберем треугольную граничную ячейку типа 9 (по классификации работы [7]).
При использовании равномерной сетки 10^10 расчет вклада в потенциал от выбранной граничной ячейки как численным методом двукратного интегрирования по методу Гаусса, так и аналитическим интегрированием дают результаты, совпадающие с точностью более чем шести знаков (см. табл. 1), что для решения уравнения Пуассона является более чем достаточным.
Формулы расчета вклада от граничных ячеек каждого из двенадцати упомянутых в [7] типов, с одной стороны, более громоздки, чем формулы расчета вклада от прямоугольной ячейки. Поэтому время расчета вклада в потенциал (или КН ЭСП) от одной граничной ячейки каждого из типов несколько больше (до нескольких десятков процентов) времени такого же расчета для прямоугольной ячейки. Однако полное число граничных
Табл. 1. Сравнение результатов расчета интегралов ¡(^ тремя методами.
Здесь Nх X Nу — параметры разбиения сетки, М — метод интегрирования: аналитический
метод (А), численный (Б), аналитически-численный (АО)
N х Ny М 10,0 I 1,0 10,1 I1,1
10 х10 A D AD -7.69371380e- 4 -7.70201629e- 4 -7.69371380e- 4 -9.39313818e- 7 -9.39311282e- 7 -9.39313818e- 7 -1.00278051e- 6 -1.00278395e- 6 -1.00278051e- 6 -5.4771779e-10 -5.4770984e-10 -5.4771779e-10
100 х 100 A D AD -1.41143729e- 5 -1.41078450e- 5 -1.41143729e- 5 -2.16196574e- 9 -2.16196775e- 9 -2.16196574e- 9 -3.33638627e- 9 -3.33646770e- 9 -3.33638627e- 9 -5.0178797e-13 -5.0178810e-13 -5.0178797e-13
Табл. 2. Примеры расчета вклада в потенциал от всего облака пространственного заряда в точке, расположенной посередине между обкладками. Расчет выполнен тремя методами для разного разбиения сетки и для различных распределений объемного заряда, расположенного между обкладками сферического конденсатора: равномерного (const), пропорционального r , r2 и r3
Nx х Ny M const r r2 r3
10 х10 A D AD 0.0655802521 0.0655803745 0.0655802521 0.0547396069 0.0547397292 0.0547396069 0.00181795542 0.00181795884 0.00181795542 0.00146042412 0.00146042754 0.00146042412
100 х100 A D AD 0.0664095014 0.0664094910 0.0664095014 0.0662341763 0.0662341658 0.0662341763 0.00188729439 0.00188729408 0.00188729439 0.00188192443 0.00188192413 0.00188192443
300 х 300 A D AD 0.0664544829 0.0664544827 0.0664544829 0.0664349726 0.0664349725 0.0664349726 0.00189311098 0.00189311097 0.00189311098 0.00189253603 0.00189253603 0.00189253603
ячеек (меньшее 2 • (Ых + Ny)), как правило, намного меньше числа прямоугольных ячеек (меньшее Nx • Ny). Поэтому для достаточно большой сетки
(Nx = Ny > 102) время расчета вклада в потенциал
(или КН ЭСП) от всех граничных ячеек может составлять единицы процента и менее времени расчета вклада в потенциал (или КН ЭСП) от всех полных (прямоугольных) ячеек.
Время интегрирования по площади граничной ячейки численным методом оказывается примерно на порядок большим времени реализации формул аналитического интегрирования.
В большинстве случаев (реально встречающихся задач) точность решения уравнения Пуассона не может быть высокой, т. к., как правило, весьма низка точность определения правой части этого уравнения. Поэтому и высокая точность интерполяции и частая сетка, в большинстве случаев, не имеют смысла.
Точность вклада в потенциал (или КН ЭСП) в сравнении аналитического и численного методов интегрирования в зависимости от разбиения сетки для полных (прямоугольных) ячеек сетки для случая плоской геометрии была рассмотрена в [1]. Для граничных ячеек точность остается практически такой же, однако полная величина вклада от всех граничных ячеек должна быть намного (вплоть до порядков) меньше полного вклада от всех полных ячеек.
В каждом из трех описанных в данной работе методов нахождения вклада в потенциал (или КН ЭСП) использована предварительная интерполяция функции плотности объемного заряда, поэтому точность во всех трех методах в значительной степени определяется точностью этой интерполяции.
Все три метода нахождения вклада в потенциал (или КН ЭСП) в совокупности ценны тем, что исключают (резко уменьшают) вероятность наличия ошибок в алгоритме. Это особенно полезно для основного рассматриваемого в данной работе
алгоритма аналитического вычисления интегралов, который реализует весьма громоздкие формулы, вывод которых и программная реализация весьма трудоемки, и маловероятно, чтобы был выполнен без ошибок. Сравнение результатов всех трех методов позволило выявить и исправить ошибки вывода и реализации.
Для демонстрации точности интегрирования тремя описанными методами интересно вычислить интегралы этими методами для разных
разбиений пространства на ячейки. Когда ТН не принадлежит ячейке сетки, по которой проводится интегрирование, то результаты интегрирования тремя методами совпадают не менее чем в девяти значащих цифрах. Поэтому сравнение таких результатов не приводится. Мы приводим результаты для случая, когда ТН принадлежит ячейке, по площади которой проводится интегрирование (табл. 1).
Для демонстрации точности вычисления вклада в потенциал от всего облака пространственного заряда приведем результаты расчета вклада в потенциал в точке, расположенной посередине между обкладками сферического конденсатора. Расчет выполнен тремя методами для разного разбиения сетки и для различных распределений объемного заряда, расположенного между обкладками: равномерного (const), пропорционального r, r2 и r3 (табл. 2).
В таблице мы не указываем размеры сферического конденсатора, потенциалы на обкладках и параметры (плотность) объемного заряда, а также координаты точки, в которой проводятся вычисления вклада в потенциал, т. к. нам важны лишь соотношения между результатами, полученными для различных параметров задачи, чтобы показать возможности получения точности.
Сравнение результатов, представленных в табл. 2, показывает, что при переходе к более частой сетке ошибки в вычислении вклада в потенциал довольно быстро уменьшаются, что объясняется эффектом улучшения интерполяции. Поэтому вполне очевидный рецепт, как уменьшить ошибки вычисления вклада в потенциал в области с ОЗ с большой неоднородностью функции распределения, — перейти в этой области к более частой сетке.
Несколько слов о возможностях и методах повышения точности решения уравнения Пуассона. Очевидно, что при любом из описанных выше методов нахождения вклада в потенциал или КН ЭСП можно выделить два метода повышения точности решения уравнения Пуассона:
• повышение точности интерполяции функции плотности ОЗ (если используется интерполяция функции ПОЗ);
• повышение точности интерполяции границы облака ОЗ.
Повышение точности интерполяции функции плотности ОЗ означает переход от билинейной интерполяции функции ПОЗ к, например, би-квадратичной интерполяции (5). Это ведет к повышению максимальных степеней переменных х и у в интегралах (8). Выражения для этих интегралов, когда одна из степеней х или у выше первой, будут приведены в следующей статье, посвященной решению уравнения Пуассона в аксиальной геометрии, для которой рассматриваемые интегралы являются базовыми. Отметим, что эти интегралы по мере повышения степеней х и у имеют все более громоздкий характер. Это ведет к увеличению времени интегрирования, которое становится сравнимым со временем аналогичного интегрирования численным методом. Т. е. теряется главное преимущество метода аналитического интегрирования.
Повышение точности интерполяции границы облака ОЗ означает переход от кусочно-линейного представления границы к более точному представлению (кусочно-квадратичному или сплайн-представлению).
При аналитическом интегрировании в случае квадратичного полиномиального выражения (или сплайн-выражения) для формы границы "внутреннее" интегрирование в интегралах (8) проводится без изменения (по сравнению с кусочно-линейным представлением границы). Во "внешнем" интегрировании возникают интегралы, являющиеся практически неинтегрируемыми. Все попытки преодолеть эту трудность приводили к такому увеличению времени счета, которое делало бессмысленным применение численного интегрирования.
При смешанном аналитически-численном интегрировании и двукратном численном интегрировании не возникало никаких дополнительных сложностей по сравнению со случаем линейного представления границы, и время вычисления вклада в потенциал или КН ЭСП практически не увеличивалось.
ЗАКЛЮЧЕНИЕ
В представленной работе описаны и реализованы три алгоритма решения плоского уравнения Пуассона для случая, когда облако объемного заряда имеет гладкие криволинейные границы. Эти алгоритмы отличаются от ранее рассмотренных в [1] тем, что в них учитываются вклады от ячеек, которые пересекаются границей облака ОЗ. Для нахождения вклада в потенциал (или КН ЭСП) проводится двумерное интегрирование по треугольным, четырехугольным или более сложным (но сводимым к треугольным или четырехугольным) областям.
Для граничных ячеек проведено сравнение аналитического, численного и смешанного (аналитически-численного) методов интегрирования при нахождении вклада в потенциал в случае, когда ТН находится вне пределов рассматриваемых граничных ячеек и когда ТН находится внутри граничной ячейки. Показано, что три рассмотренных метода интегрирования дают вполне удовлетворительные по точности результаты (особенно метод аналитического интегрирования и смешанный метод аналитически-численного интегрирования).
Метод аналитического интегрирования при нахождении вклада в потенциал при билинейной интерполяции функции плотности ОЗ и кусочно-линейном представлении границы дает значительный (до порядка) выигрыш во времени счета. Однако более точная интерполяция плотности ОЗ настолько увеличивает время счета, что делает неце-
Приложение. ВИД ИНТЕГРАЛОВ I,,. (,, у = 0,1)
лесообразным применение данного метода.
Метод смешанного интегрирования выполняется до трех раз более быстро, чем численный метод, и дает вполне удовлетворительную точность. Однако требует интерполяции как функции плотности, так и границы облака ОЗ. Переход к более точной интерполяции любого из этих объектов не привносит в алгоритм трудностей и не ведет к заметному увеличению времени счета.
Численный метод (двукратного численного интегрирования по формуле Гаусса) является самым медленным из рассмотренных. И по точности он уступает вышерассмотренным методам. Однако в рамках этого метода возможно не только использование более точной интерполяции для функции плотности ОЗ и для границы облака ОЗ (чем в аналитическом методе), но и применение для них аналитических функций.
Ун )_'
Обозначим dx _ хк _ хн, Тогда
__ __ =2 _в/ —к _ Хк х1, _ Х1 , Х1 _а Х2 , Х2 _ /(а2 + 1)
(К _ ХН )+в (Хк _ Хн ) ■ [1п(«2 + 1) _ 2] _ Ун ■ [ Хк 1п(хК + У2) _ Хн 1п(хН + Ун) ] +
+ 2 ■ dx _ 2 у
_ 1п( —2н + х1)
+ 2(в_а^ Х1)х,
аг^ Х— _ аг^
Ун Ун
+ 1п( —2К + Х22) ■
—(¿"к + х\ ) + (£_«■ Х1) Zк
а
—(¿н + Х22) + (в _а ■ Х1) 2
— 2 аг^ —— _ аг^—
_—(г2к _ _ 2(в_а■ х1)(2к _ —н ) +
2 а ■ хк + в 2 а ■ хн + в Хк аг^-^- _ х2 аг^-
в
dx _ х2 аг^ — + х2 аг^ —
[1п(—к + х2) _ 1п(+х2)] +
а +1-1
х2 ■ 1
а2 +1 х2
1 — — аг^ —— _ агй§ —— + ( х2к аге1в ^ _ х2 аге1в ^^ + Ун ■ dx _ Ун ■
Х2 _ Х2 Х2 ] ^ хк хн
аг^—— _ аг^-
Ун
11,0 _
а (к _ хЪн ) + в (( _ х2 ) +1) _ 2] + Ун■(( _ х2н )_
_ Ун [(( + Ун )1п(Хк + У2н) _ (н + У2н )1п(4 + У2н) _ (( _ хн )] +
—к 1п( —к + х\) _ —^Щ—н + х22) __ х2) + 2х2 ■ dx _ 2х2
( — 2 ^ аг^ —— _ аг^ ——
Х2 Х2 yJ
+ —
3
Хк аг^
а ■ хк + в
х
_ Хн аг^
а ■ хн + в
Л
к
в
2в
3(а2 +1)
{ у? \ {
12 Х2
2
аг^ —— _ аг^—
3(а +1) ^ 2в
(_ х2+зх2) [1п(—к + х2) _ 1п(—н+х2)] +
3(а2 +1)
^ г2 _ г2 ^^ _ 3х, ■ dx
х1 аге1в^ _ хън аге1в^ _ У-[1п(хк + у2) _ 1п(х2 + У2)] + У-(х2 _ х2);
Хг
■У
ХАТ
I0 1 = yAdx -1 0Д 2 2
а
3
- + ав(x2K - xN) + вdx
ln(a2 +1)
(а2 +1) Хк з Xn + ав(x2K - x2N) + в2dx
(а2 +1)
i-V
' 3 2 Л
XK yNXK
f z3 ^ + X 2 z
3 -гл.^ K
\
/
..2 i ,,2
ln( zK + X22) -
f x-3
fz 3
+x2 z 3 -г N
V
ln( zN + x2) - zK - zN)
(а2 +1)
X 2 dX + X3 3X2 x 3 X2
arctg —— - arctg —1N
in( xK + y2)+
2 Л
XN yNXN
\
ln( xN+yN) -
2 yN f
Л f 3 - 3 2
XK : XN XK XN 1 2y N
arctg—— - arctg
2
2 Л
dx
11,1 ^ (xk XN) а
+ (а2 + 1)ln( zK + x2)
--ав-
6 3
4 - 4 3
zK X2 X1z к
-в'
- + 0.51п(а2 +1)
(а2 +1)-
■2ав:
в2
- (а2 + 1)ln(zN + x2)
- 3 x1 х2(а +1)
arctg —— - arctg —N
Л
+ (а2 +1)
^2" K 2
4 - 4 .. - -
zк zN J2 zк zN
ZM х х Z л.
16
х zk zn - ххг ^х + х х2dx
х 9 Х1 Х2 3 Х1Х2 х
/
СПИСОК ЛИТЕРАТУРЫ
1. Иванов В.Я., Шевченко С.И. О расчете плоских электростатических полей в приборах, имеющих области, заполненные объемным зарядом // Научное приборостроение. 1999. Т. 9, № 4. С. 88-94.
2. Шевченко С.И. О расчете аксиально-симметричных электростатических полей в областях, заполненных объемным зарядом // Научное приборостроение. 2002. Т.12, № 2. С. 23-29.
3. Соболев С.Л. Уравнения математической физики. М., 1966. 400 с.
4. Ильин В.П. Численные методы решения задач электрофизики. М.: Наука, 1985. 334 с.
5. Иванов В. Я. Методы автоматизированного проектирования приборов электроники. Новосибирск: Изд-во СО АН СССР, 1986. 193 с.
6. Шевченко С.И. Алгоритм получения предельной точности в электростатических расчетах элементов электронно- и ионно-оптических
приборов, имеющих плоскую симметрию // Научное приборостроение. 1997. Т. 7, № 1-2. С. 45-53.
7. Шевченко С. И. О расчете траекторий заряженных частиц вблизи криволинейных поверхностей // Научное приборостроение. 2005. Т. 15, № 1. С. 70-79.
8. Корн Г. и Корн Т. Справочник по математике. М., 1970. 720 с.
9. Абрамовиц М., Стиган И. Справочник по специальным функциям. 1979. 830 с.
10. Крылов В.И., Шульгина Л.Т. Справочная книга по численному интегрированию. М.: Наука, 1976. 370 с.
Институт аналитического приборостроения РАН, Санкт-Петербург
Материал поступил в редакцию 13.12.2005.
ON THE SOLUTION OF TWO-DIMENSIONAL POISSON'S
EQUATION FOR THE CASE OF A SPACE-CHARGE CLOUD WITH SMOOTH CURVILINEAR BOUNDARIES.
I. PLANE GEOMETRY
S. I. Shevchenko
Institute for Analytical Instrumentation RAS, Saint-Petersburg
The paper presents an algorithm for solving two-dimensional Poisson's equation in plane geometry by the boundary integral equation method for the case of space charge cloud with smooth curvilinear boundaries.
