Расчет трехмерных электростатических полей методом граничных элементов с выделением сингулярностей ядра около поверхностей электродов Текст научной статьи по специальности «Физика»

ISSNÜ868-5886

НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2ÜÜ4, том 14, № 4, c. 2Ü-38

ОРИГИНАЛЬНЫЕ СТАТЬИ

УДК 519.635.4: 621.319.7 © А. С. Бердников

РАСЧЕТ ТРЕХМЕРНЫХ ЭЛЕКТРОСТАТИЧЕСКИХ ПОЛЕЙ МЕТОДОМ ГРАНИЧНЫХ ЭЛЕМЕНТОВ С ВЫДЕЛЕНИЕМ СИНГУЛЯРНОСТЕЙ ЯДРА ОКОЛО ПОВЕРХНОСТЕЙ ЭЛЕКТРОДОВ

В работе рассмотрена модификация метода граничных элементов, предложенная для вычислений с высокой точностью (относительная ошибка < 10-5 ) трехмерных электростатических полей при использовании электродов практически произвольной формы. Продемонстрированы и проанализированы эксперименты с аналитическими тестовыми задачами.

ВВЕДЕНИЕ

Хотя на рынке программного обеспечения уже имеются программы, предназначенные для расчета трехмерных электростатических полей [1-8], особенно точные вычисления в оптике заряженных частиц по-прежнему требуют более высокой точности, чем точность, достигаемая этими программами. Более того, поскольку для всех сторонних пользователей эти программы недоступны на уровне исходных кодов, это затрудняет построение гибких и быстродействующих специализированных приложений на их основе. Чтобы справиться с указанными трудностями, была разработана специализированная библиотека модулей, предназначенная для точного вычисления электростатических полей методом граничных элементов, которая может быть интегрирована с программами пользователя в единый выполняемый модуль. Предварительное тестирование показывает, что предложенные численные алгоритмы могут вычислять электростатическое поле с точностью 10-5 вплоть до поверхности электродов с разумными расходами компьютерной памяти и компьютерного времени1-1. Сравнение с другими программами показало, что на сегодняшний день эта точность является рекордной (конечно же, при условии, что разбиение границ на конечные элементы производится с примерно таким же шагом дискретизации). Численные эксперименты показали, что сходимость разработанного алгоритма равна 4), где И — шаг разбиения поверхностей электродов на отдельные граничные элементы. В данной статье описываются некоторые характерные особенности разработанного алгоритма, которые позволяют достичь таких результатов.

1 На самом деле в силу специфики метода граничных элементов точность 10-5 на поверхности электродов означает, что вдали от электродов точность расчета много лучше.

1. ОСНОВНЫЕ ПРИНЦИПЫ МЕТОДА ГРАНИЧНЫХ ЭЛЕМЕНТОВ

Среди общеприменимых известных методов, предназначенных для решения уравнения Лапласа и расчета электростатических полей, а именно метода конечных разностей, метода конечных элементов и метода граничных элементов, только метод граничных элементов (он же — метод интегральных уравнений, или метод зарядовой плотности) обеспечивает необходимую высокую точность [9]. Хотя с точки зрения классического математика метод граничных элементов позволяет решить уравнение Лапласа со сколь угодно высокой точностью (то есть с нулевой погрешностью в случае идеальных плотностей заряда), компьютерные реализации метода неизбежно привносят численные ошибки. Вот список основных источников погрешности расчета поля в случае применения метода граничных элементов.

• Приближенное представление точных границ с помощью дискретизованных поверхностей.

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

• Отклонения дискретно вычисленных интегральных выражений, используемых в методе граничных элементов, от математически точных аналитических значений интегралов.

• Резкое увеличение погрешности численного интегрирования вблизи граничных поверхностей и на граничных поверхностях в силу сингулярного (нерегулярного) поведения интегрального ядра около граничной поверхности.

• Приближенное представление континуального граничного условия набором дискретных коллокационных точек или каким-либо еще дискретизованным аналогом истинного граничного условия.

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

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

• Увеличение влияния неточности вычисления поверхностного заряда для точек, расположенных вблизи граничных поверхностей, в силу увеличения весовой функции (интегрального ядра) в этих точках.

• Увеличение погрешности в аппроксимации "истинной" плотности заряда вблизи краев электродов, которое является следствием сингулярного поведения плотности заряда в этих точках (поэтому плотность нельзя приблизить непрерывными дискретизованными поверхностными функциями без потери устойчивости метода, если не принимать специальных мер).

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

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

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

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

П( р, д) =

х = X (р, д),

у = у (р, д),

г = 2 (р, д),

то есть для каждой точки (р, д) в пространстве параметров (будем считать, что (р, д) принадлежит некоторой прямоугольной области ръ < р < ръ, дъ < д < дъ , специфичной для данной поверхности П) мы можем вычислить трехмерную точку (X, У, 2), которая принадлежит поверхности П, а также вычислить касательные вектора, нормальный вектор и другие дифференциально-геометрические элементы поверхности, если они будут нам необходимы.

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

а(р, д) аП(р, д)

х - X(р, д))2 + (у - У(р, д))2 + (г - 2(р, д))2

а(р, д) (X(р, д) - х) аП(р, д)

^(х - X(р, д))2 + (у - У(р, д))2 + (г - 2(р, д))2 ^

^(р, д) (У(р, д) - у) ¿П(р, д)

^(х - X(р, д))2 + (у - У(р, д))2 + (г - 2(р, д))2 ^

^ (р, д) (2(р, д) - г) ап(р, д)

л/(х - X(р,д))2 + (у - У(р,д))2 + (г - 2(р, д))2

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

Хотя в начале нашего вычислительного процесса мы и не знаем поверхностную плотность зарядов, но нам известны граничные условия (например, значения потенциала), заданные на поверхностях. Эта информация в принципе позволяет восстановить поверхностные плотности зарядов (процедура, как это можно сделать, описывается, например, в [10]). А как только плотность зарядов становится известной, для произвольной точки (х, у, г) трехмерного пространства с помощью приведенных выше интегральных формул можно вычислить потенциал и компоненты напряженности электростатического поля (предполагая, что трудности с вычислением интегралов отсутствуют).

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

2 , 2,2. \ конечно удаленной точке х + у + г ), эти

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

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

1.1. Поверхности

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

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

3

3

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

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

Чтобы достигнуть по возможности наивысшей точности расчета, будем рассматривать граничные поверхности "как есть" — без приближения более грубыми, но и с большей легкостью интегрируемыми элементами. Этот подход имеет свои недостатки (см., например, [5], где вполне убедительно аргументируются преимущества фасеточных поверхностей при описании сложных объектов), но зато позволяет избежать очевидных ошибок, связанных с дискретизацией исходной поверхности. На всех уровнях алгоритма, кроме самого нижнего, граничные поверхности рассматриваются как аналитически точные "черные ящики": граничные поверхности описываются с помощью абстрактного параметрического представления, реализованного как набор подпрограмм, и формулы, "зашитые" внутри подпрограмм, извне — основной части алгоритма, вообще говоря, неизвестны. То есть, если мы можем вычислить значения параметризующих функций X(р, д), У(р, д), 2(р, д), а также всех их производных в произвольной точке (р, д), мы знаем поверхность.

С точки зрения организации программного кода библиотека модулей включает в себя замкнутую подсистему — библиотеку трехмерных геометрических объектов и библиотеку трехмерных геометрических преобразований. Библиотека трехмерных геометрических объектов реализует подпрограммы функций X(р, д), У(р, д), 2(р, д) для элементарных поверхностей типа плоского прямоугольника, плоского диска, кругового коль-

ца, цилиндрического, сферического или тороидального сегмента, цилиндрического или гиперболического цилиндра и т. д. (в эту библиотеку можно добавлять любые новые элементы, для которых существует параметрическое представление через аналитически заданные функции, обладающие достаточной степенью гладкости). Библиотека трехмерных преобразований позволяет перемещать элементарные объекты в нужную точку пространства и с нужной ориентацией относительно базовой координатной системы, комбинируя элементарные перемещения (параллельный перенос, повороты и симметрии). Необходимо подчеркнуть, что функции X(р, д), У (р, д), 2(р, д) должны быть достаточно гладкими, чтобы допускать вычисление производных высокого порядка (вплоть до четвертого включительно) — эти значения необходимы, чтобы выделять сингулярности интегрального ядра для пробных точек, расположенных около поверхности и на поверхности (см. далее разделы 1.3, 2.2, 2.3.).

1.2. Поверхностная плотность зарядов

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

Стандартный подход — использовать дискретную аппроксимацию, при которой функция ст(р, д) заменяется конечной суммой с неизвестными коэффициентами от известного и фиксированного набора базовых функций

a(p, д) д).

г. ]

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

Использование для поверхностной плотности заряда дискретной аппроксимационной схемы фиксированного порядка было бы неизбежным,

если бы не имелось более гибкой и более успешной идеи. На самом деле нам не требуется знать поверхностную плотность заряда — нам надо знать интеграл от поверхностной плотности заряда. При предположении, что поверхностная плотность заряда это произвольная (хотя и достаточно гладкая) функция, интегрирование может быть только численным. Численное интегрирование (одномерное и двумерное) базируется на приближенных формулах типа конечных сумм:

ь

I = Г /(х)ах = ^ /(х ^

а 1 =1,^

а ь

^ = | &I&х f(x, У) “ ?(х и) =

с а 1=1,N

1 =1М

_ ^Л1Л]f (хг, у}).

1=1, N

1 =1,М

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

)_ГГ, =

->/(х - X (Р,Я))2 + -о аП(Ру , Я и )

= ^а1, Пч I-----------------2------

и 1 х -X(р1,^)) + ...

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

2)

На самом деле это не совсем так. Порядок сходимости интегралов ограничивается, помимо всего прочего, еще и тем фактом, что подынтегральное выражение не является полиномом и, следовательно, сходимость про-

Достоинства и недостатки "квадратурной " аппроксимации заряда

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

• Мы ограничили себя фиксированной квадратурной формулой и потеряли возможность гибко выбирать алгоритм интегрирования в зависимости от положения пробной точки. Это не очень страшно, поскольку (в случае необходимости) всегда можно построить дискретную аппроксимацию нужного порядка по известным узловым значениям поверхностной плотности и пересчитать зарядовую плотность для промежуточных узлов. Но в таком случае нужно снова строить явную аппрок-симационную схему по заданным узлам — и это означает, что если исходной квадратурной формулы вдруг оказывается недостаточно, то в концептуальном смысле мы ничего не выиграли.

• Мы однозначно связали число неизвестных, характеризующих распределение заряда на заряженной поверхности, с числом узлов, необходимых для точного вычисления интегралов. Это не самый правильный выбор, поскольку сплошь и рядом поверхностную плотность заряда можно с достаточной точностью аппроксимировать с гораздо меньшим количеством узловых значений, чем требуется квадратурных узлов для вычисления поверхностных интегралов (с той же, естественно, точностью, что и точность аппроксимации поверхностного заряда). Поскольку число неизвестных, связанных с поверхностью, — узкое место метода граничных элементов (размер полной

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

• На самом деле мы никуда не ушли от дискретной аппроксимации с использованием базисных функций — просто теперь она спрятана внутри квадратурной формулы (коэффициенты квадратурной формулы как раз и выводятся в предположении, что интегрируемая функция заменяется полиномом, построенным по значениям функ-

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

ции в узлах интегрирования (р 1}-,д 1}-)). Большинство проблем, связанных с построением для поверхностной плотности заряда хорошей аппрок-симационной схемы, сохраняются и при таком выборе неизвестных (например, размещение фиксированного числа узлов так, чтобы оптимальным образом приблизить данную плотность заряда, или расходимость аппроксимации при сильной разно-масштабности разбиения).

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

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

1.3. Сингулярные и нерегулярные

компоненты интегрального ядра

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

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

Стандартный способ, как это сделать: представить поверхностную плотность в виде усеченного ряда Тейлора, выделить сингулярные составляющие и проинтегрировать их отдельно — аналитически или с помощью специально выбранной квадратурной схемы (подобный вариант компенсации сингулярных компонентов используется, например, в [11, 12]):

и _ ГГ_Д(р,£)4р£^_ о ^Дх2 + Ду2 + Дг2 (сг(p,д) - Стр -Стррр -афд - ...)ф с!д +

д/Дх2 + Ду2 + Дг2

+ст„ II ар * +

О д/Дх2 + Ду2 + Дг2

, _ гг р Ф ад ,

+ °р0ГГ I 2 2 2 +

а а/Дх +ДУ +Дг

, _ гг д Ф ад ,

+ °д0.!■! / 2 2 2 + •••

а У Дх + Ду + Дг

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

Этот способ, который мы назвали библиотекой эталонных сингулярностей, состоит в следующем. Пусть имеется набор аналитически заданных функций фк (р, д, х, у, г), которые, с одной стороны, могут быть проинтегрированы по Р и д в аналитической форме, а с другой стороны, обладают требуемым сингулярным поведением в окрестности точки (р0, д0) (в разделе 2.3 показано, как эти функции могут быть сконструированы в явном виде). После этого мы можем вычесть из исходного подынтегрального выражения все сингулярности, которые мешают численному интегрированию, подбирая должным образом коэффициенты Хк:

и =

И

р, ч)

л

д/Дх2 + Ду2 + Д

Xя* Я^*(P, Ч’Х’ У’г) Ф аЧ •

-ХЯ*^*(Р’ Ч’Х У’г)

ар ач+

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

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

1.4. Сингулярные "всплески" плотности заряда вблизи краев

Важной особенностью метода граничных элементов применительно к задачам электростатики является то, что поверхностная плотность зарядов, как правило, имеет сингулярные особенности на краях электродов [13]. Эта сингулярность имеет другую природу, нежели сингулярность интегрального ядра при приближении пробной точки к заряженной поверхности, и должна приниматься во внимание при аккуратном расчете поля. Практические эксперименты показывают, что для таких случаев попытки найти хорошее решение для поверхностной плотности заряда в классе полиномиальных или сплайновых функций приводят к чис-

ленной неустойчивости. Поэтому хороший алгоритм должен специальным образом обрабатывать края электродов (см., например, [5]), однако полноценное рассмотрение этого аспекта метода граничных элементов требует отдельной работы.

2. ДЕТАЛИ АЛГОРИТМА

2.1. Базовые сингулярности интегрального ядра

Рассмотрим интегральное ядро, соответствующее потенциалу поверхностного заряда:

К (, Я,х, У,2 ) =

=______________________1_____________________

л/(х (а Я )- х )2 + ( ( Я)- У)2 + ( (Р, Я)- 2 )2

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

X (р, я )= Х о + Хро (Р - Ро)+ Хяо (Я - Яо)+ -,

У(Р, ЯЬ Уо + ¥ро(Р - Ро)+ ¥яо (Я - Яо)+ -,

% (Р, Я) = % о + %ро (Р - Ро)+ ^яо (Я - Яо)+ -

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

+

к (р, Ч У

УІА(р + рн )2 + в{ч + чн )2 + 2С(р + рн )(ч + чн )+ Е __________________________________1__________________________________

у/а2 •(р + рн )2 + Ь2 -р + Чн )2 + 2-с-а- (р + рн )-Ъ • ( + Чн )+е2

1

где а, Ь , с, рн , Ян , £ — это вещественные значения, выражающиеся через значения функций X (р, я ), У (р, я), % (р, Я) и их производных в точке

(р0 ’ Ч0)’ а также через координаты пробной точки (Х’ У’ г). При этом легко показать’ что для невырожденной точки поверхности справедливы

условия

А = а2 > 0, В = Ь2 > 0, Е = £ > 0,

т. е. запись соответствующих коэффициентов как квадратов некоторых вещественных чисел вполне

обоснована, а коэффициент с = С/^АВ удовлетворяет условию -1 < с < +1.

Это выражение эквивалентно записи К(р, д) ~

1 _ _

:, где переменные р ид — это

р2 + д2 + £2

нормализованные отклонения от точки (р0, д0), а £ является нормализованным расстоянием между пробной точкой (х, у, г) и ближайшей к ней точкой поверхности (р0, д0). Разложив плотность поверхностного заряда в ряд Тейлора по переменным р ид , получаем следующий набор компонентов интегрального ядра, которые являются сингулярными и нерегулярными (точнее, могут являться сингулярными и нерегулярными) при приближении пробной точки к граничной поверхности:

1р

Гр 2 + д2 д + £ 2 '

у1р 2 д рд + + £2

у1р 2 + д2 р3 + £2

р2

ч2

л/]

р2 д

ны при малых значениях £ отличаются от полиномиальных функций значительно (и следовательно, создают трудности при численном интегрировании), в то время как кубические компоненты выглядят вполне регулярными.

2.2. Ошибка интегрирования и сходимость процесса численного интегрирования для сингулярных компонентов интегрального ядра

В этом разделе мы проверим, как работает процедура численного интегрирования для сингулярных компонентов интегрального ядра, рассмотренных в предыдущем разделе. В Приложении на рис. П2-П5 показано отклонение результата численного интегрирования от аналитически точного значения для различных значений £ в случае, когда £ ^ 0. Для рассмотренных в этом примере вариантов максимальная относительная ошибка между численным результатом и аналитическим значением приведена ниже в виде пары (вид функции, ошибка):

•у/р2 + д2 +£2 д/р2 + д2 + £2

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

3) Вообще говоря, чтобы выполнить корректно вычитание сингулярных членов вплоть до членов второго порядка, необходимость которого обоснована в разделе

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

X (р, д), У (р, д), 7 (р, д). В этом случае указанные условия будут выполнены, если пробная точка (х, у, г) находится достаточно близко от точки (р0, д0), но они могут нарушаться, если точка (х, у, г) расположена слишком далеко от точки (р0, д0). Проверка, что квадратичная форма в знаменателе интегрального ядра все еще является положительно определенной, является задачей программного блока, ответственного за вычитание сингулярных и нерегулярных компонентов.

1/д/ р2 + д2 + £2 — 1.2 %,

рЦр2 + ч2+£2 — 0.2 %,

дЦр2 + ч2 +£2 — 0.2 %,

р2Цр2 + ч2+£ — 0.015 %,

д2Ц р2 + д2 + £2 — 0.015 %,

рчЦр2 + ч2 +£2 — 0.023 %,

р 3/л1 р2 + ч 2 +£ — 0.003 %,

р2 + д2 +£2 — 0.003 %,

р2 ч/^р 2 + ч 2 +£2 — 0.006 %,

рч 2Ц р2 + ч 2 +£2 — 0.006 %.

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

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

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

для функции і/д/ Р2 + Ч2 + £2 — « о(л),

—р/-у/р 2 + Ч2 +Є2 — = о(й2

—”— ч/4р 2 + ч2 + £І — = о(й2

—''— Р 2Л/р2 + ч2 +£2 — = о(и3

—”— РчЦр2 + ч2 +е2 — = о(й3

—”— ч ЦР2 + ч2 +е2 — = о(й3

—”— Р 3Л/Р2 + ч2 +е2 — = о(л4

—”— Р2 ^Р2 + ч2 +е2 — = о(л4

—”— рч 2/4р 2+ч 2+£І — = о(л4

—”— ч ЪЦ Р2 + ч2 +Є2 — = о(л4

Это означает, что если из исходного интегрального ядра аккуратно выделить нерегулярные компоненты нулевого, первого и второго порядков, то сходимость процедуры численного интегрирования для остатка будет равна ~ о(к4) независимо от порядка гауссовской квадратурной формулы, использованной для интегрирования (при условии, что порядок гауссовской квадратурной формулы не ниже 2). То есть, если ошибку численного интегрирования определить как = Скт , то увеличе-

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

Этот вывод является важным, поскольку показывает, что для аппроксимации поверхностной плотности заряда, если в ней возникает необходимость, достаточно использовать схему 4-го порядка — ведь все равно преимущества аппроксима-ционной схемы более высокого порядка будут испорчены более низким порядком сходимости процедуры численного интегрирования. Более детально этот эффект можно проследить на процедуре выделения "плохо интегрируемых" компонентов:

и =

я

о( Р, ч)

~°09(0) (Р, ч, X, У, 2)~

д/Дх2 + Ду2 + Д2

ч

- °роР(Р)(Р, ч, X У,2) - ^чо^(ч)(Р, ч, X У,2) -

- °рроР(рр) (Р, q, X ^ 2) - °рчо9(РЧ) ч, x, У, 2) -

Л

+ о,

+ о

+ о

~^ддо9(дд)(Р, Ч, X У, г) - ••• аР ¿д +

У

Ц^(0) (Р, Ч, X У, г)аР аЧ +

а

ар00 ИР(р) (р’ Ч, X У, г)аР аЧ +

а

!К(Р’ Ч, X У, г)аР аЧ +

а

р0 И Р(РР)(Р, Ч,X У, г)аР аЧ +

а

+ аРЧ0 ИР(РЧ) ^ Ч, x, У, г)аР аЧ +

а

+ ^ЧЧ0 II Р(чч)(Р, Ч, x, У, г)аР аЧ + •

а

Пусть функция о(р, Ч) нам известна с точно -стью о(ит), а погрешность численного интегрирования после выделения сингулярных членов оценивается как О {к" ). Кроме того, нам известно, что ошибка численного интегрирования по той же квадратурной формуле для функции Ф(0)(Р, Ч, х, У,2) составляет О(к), для функций

Ф(р)(Р,Ч,х,У,г) и Ф(Ч)(р,Ч,х,У,г) — о{к2), для функций ф( рр)( р, ч, х, У, z), Ф( рЧ)( Р, Ч, x, У, г) и Ф(чч)( Р, Ч, x, У, г) — о{к 3) и т. д.

Если бы численное интегрирование для всех этих выражений было абсолютно точным, то при любых коэффициентах при сингулярных членах получалось бы одно и то же выражение — исходный интеграл с ошибкой О{кт), соответствующей неточности поверхностной плотности заряда. Если бы поверхностная плотность заряда и коэффициенты при сингулярных членах были бы известны абсолютно точно, то ошибка интегрирования была бы равна О {к" ) — т. е. неточности численного интегрирования, соответствующей первому сохраненному сингулярному члену. Однако при комбинировании этих ошибок получается более интересная картина.

Действительно, в силу неточности аппроксимации поверхностной плотности заряда коэффициент О0 известен с точностью о{ит ), коэффициенты О 0 и Оч0 — с точностью о(кт 1), коэффициенты Оpp0, Орд0 и ОЧЧ0 — с точностью о(кт-2) и т. д. Это означает, что в исходном выражении, которое интегрируется численно, остались сингулярные компоненты, вклад которых пропорционален ошибке соответствующего коэффициента. В результате для нулевого сингулярного члена появляется дополнительная ошибка численного интегрирования о{кт )• о{к) (произведе-ние неточности интегрирования на амплитуду сингулярного члена), для линейных членов —

ошибки о{кт-1 )• о (к2 ), для квадратичных членов — ошибки о{кт~2 )• о (к3 ) и т. д. Все эти ошибки оказываются меньше (в плане порядка сходимости), чем вклад от неточности значений

поверхностного заряда, равный о{кт ). Тем самым ошибка вычислений при комбинировании двух ошибок оценивается как тах{о{кт ) о{к")}, что и требовалось доказать. (Из этих рассуждений, в частности, следует, что ошибка в определении

коэффициентов О0 , Оp0, ОЧ0, Оpp0, ОРЧ0, ОЧЧ0 в меньшей степени влияет на неточность результата и что, вообще говоря, можно вычислять значения этих коэффициентов достаточно грубым образом, не добиваясь аналитически точного значения).

2.3. Библиотека эталонных сингулярностей

Чтобы двигаться дальше, нам требуются новые аналитические выражения, а именно функции рк (р, Ч, х, У, I), которые могут быть проинтегрированы в аналитическом виде и которые ведут себя при приближении к граничной поверхности как сингулярные компоненты интегрального ядра, т. е. как линейные и квадратичные члены, ведущие себя нерегулярным образом при выполнении численного интегрирования и рассмотренные ранее в разделах 2.1 и 2.2. (В разделе 1.3 было показано, как при наличии таких функций можно сделать процесс численного интегрирования интегрального ядра равномерно регулярным без искусственного увеличения числа узлов интегрирования вблизи граничной поверхности).

Базой для библиотеки эталонных сингулярностей служат следующие функции, комбинируя которые с нужными коэффициентами можно добиться строгого соответствия результата двукратного дифференцирования по переменным Р и Ч нужному сингулярному поведению:

Я( р, ч) = д/ а2 (р + ръ )2 + Ъ2 (ч + чк )2 + 2с • а(р + Ръ)' ъ( + чк )+£2 РР(Р, ч) = \°$а{р + Ръ) + с • Ъ(ч + чъ) + К{Рг ч)], Р<2(Р, ч) = 1°§[ъ(ч + чк)+ с •а(( + Ръ) + R(Р, ч)]

РР(Р, ч) = агс!ап[2 • а(р + [)• ъ(ч + чк)- с •є2•є • ^Р, ч))];

здесь а > О, Ъ > 0, є> 0, -1 < с <+1

= V1 - с2 — это константы, выражающиеся через характеристики параметризованной поверхности в окрестности точки (р0, ч0), ближайшей к пробной точке (х, у, г) (см. раздел 2.1). В результате достаточно трудоемких вычислений, выполненных

с применением программы МмИвтаИса версии 5, получены выражения для эталонных сингулярных функций. Они приведены в Приложении (можно проверить, что для приведенных функций действительно выполняются дифференциальные соот-

ношения

д2Фтп (Р, ч)дРдч = Фтп (Р, ч)).

2.4. Модифицированные эталонные сингулярности: поведение на бесконечности

Суммируем результаты, полученные в предыдущих разделах. Нами были построены аналитические функции (ртп (X, у, г, р, р), которые ведут себя вблизи граничной поверхности как

Ртп (X У, ^ P, Ф =

= ртрп • [а2 (р + рк )2 + Ь2 (р + рк )2 +

+ 2с •а • (р + рк )Ь • (р + рк)+ е2 ] V2

и для которых существуют выражаемые в аналитической форме интегралы Ф тп (х, у, г, р, р), удовлетворяющие соотношениям

¡¡Ртп (р- Р, X, У, 2) Ф ¿Р = Ф тп (X У, р, Р).

С помощью этих функций потенциал и (х, у, г) может быть записан как

и(х, у, г ) = Ц

ст( р, р)

д/дх2 + Ду2 + Дг2

~^^тпРтп (р, Р, X у, 2)

Ф ¿р + ^ Атп ¡¡Ртп (р5 Р, X у> 2) Ф dР,

где коэффициенты Хтп зависят от поведения зарядовой плотности в окрестности точки поверхности (р0, р0), ближайшей к пробной точке

(х, у, г ), а функции ртп (р, р, х, у, г) зависят и от точки поверхности (Po, Чо), и от пробной точки (х, у, г), и от дифференциальных характеристик поверхности в окрестности точки (Po, ро ), и, что самое существенное, от расстояния между пробной точкой (х, у, г) и ближайшей к ней точкой поверхности (р0, р0). Первый член этого выражения ведет себя при численном интегрировании как регулярное выражение, а оставшиеся члены могут быть проинтегрированы в аналитической форме, несмотря на содержащиеся в них сингулярности:

¡¡Ртп (А ^ X У, г)(^р ¿р =

= Ф тп ( (, ^ ртах , ртах ) - Ф тп ( У, ^ р,шп, ртах ) -

- Ф тп ( У, ^ ртах, р^п ) + Ф тп (x, (, ^ р^п, р,шп )

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

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

ожиданные проблемы. А именно функции Фтп (х, у, г, р, р), определенные в разделе 2.3, стремятся к бесконечности при удалении пробной точки (х, у, г) от граничной поверхности, что вызывает значительные ошибки в процессе вычисления аналитических интегралов от сингулярных компонентов интегрального ядра, поскольку приводит к вычитанию больших величин, отличающихся друг от друга на незначительную величину.

Чтобы избежать накопления ошибок округления, возникающих при вычитании близких по значению чисел с плавающей точкой, функции Фтп (х, у, г, р, р) из раздела 2.3 необходимо модифицировать так, чтобы они стремились к нулю или хотя бы к константе, когда пробная точка (х, у, г ) стремится к бесконечности:

ИтФтп ( ^ X ^ г)= Oр),

х——^

ИтФтп ( ^ X ^ г)= ^

у—~>

Ит Фтп( q, х y, г )=°(1).

Для того чтобы сохранить дифференциальные соотношения д2 Фтп (р, р )дрдр = Ртп (р, р) между

функциями Фтп (х, у, р, р) и Ртп ^ Ц, X, у, z), и при этом обеспечить нужную асимптотику на бесконечности, достаточно применить алгебраическое преобразование

Ф тп (X у> ^ P, р) =

= Ф тп (X у, ^ P, р) - рФ Ртп (X у, ^ р) -

- рФ <2тп (X, у, ^ р) + Ф1 (X, у, г)

при выбранных подходящим образом функциях

Ф1 (X У, Z, Р), ФQmn (Х, У, Z, q) и Ф0mn (Х, У, Z)•

(При этом последняя чисто техническая трудность состоит в том, что функции Фmn (x, y, z, p, q) необходимо привести к виду, когда при их вычислении не происходит вычитания друг из друга больших и близких по значению чисел). Эта задача легко может быть выполнена любым опытным математиком, поэтому результирующие выражения в данной работе не приводятся, чтобы не загромождать текст лишними формулами.

3. ЧИСЛЕННЫЕ РЕЗУЛЬТАТЫ

Рассмотрим несколько тестовых примеров, чтобы продемонстрировать эффективность предложенного алгоритма.

1. Однородно заряженная сфера. Сферическая поверхность описывается с помощью стандартных формул сферической системы координат как x = Rcosp-cose , y = Rsinp- cose , z = = R sin в . Потенциал на поверхности сферы постоянен и равен U0. Потенциал внутри сферы

а

б

Рис. 1. Тестовый пример — однородно заряженная сферическая поверхность. Результаты: а) радиальное распределение потенциала; б) отклонение между аналитической формулой и численно рассчитанным значением (максимальная ошибка меньше 10-7)

f

О 10 20 30 40 50 60 70 80 90

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

(т. е. при х2 + у2 + г2 < Я2) также равен ио, потенциал вне сферы (при х2 + у2 + г2 > Я2) убывает по формуле и (х, у, г )= и 0 яЦх2 + у

2 2 2

2 + y + z

а

б

в

г

Рис. 1 демонстрирует результаты численных вычислений. Поскольку однородный заряд описывается идеально точно выбранной схемой аппроксимации поверхностной плотности заряда, то ошибка, которая видна на графике — это ошибка численного интегрирования поверхностной плотности заряда в чистом виде. Как и ожидалось, наибольшая ошибка обнаруживается на границе (т. е. на поверхности сферы), но благодаря аккуратному вычитанию сингулярных компонентов интегрального ядра она достаточно мала (< 8 -10-8). (Стоит также заметить, что эта точность получена при сетке зарядов на поверхности сферы всего лишь 8х8 !).

2. Сферическая поверхность, заряженная как идеальный диполь. Сфера параметризована как x = R cosp- cos0 , y = R sinp- cos0 , z = R sin0 . Поверхностная плотность заряда ведет себя как ~ sin(0). Потенциал вне сферы ведет себя как потенциал идеального диполя и описывается форму-

лои

U (x, y, z )= U 0 zR y |y

Потен-

циал внутри сферы соответствует однородному полю и описывается формулой U(х, y, z) = U0 z¡R . Потенциал на поверхности сферы ведет себя как U (в)= U 0 sin(0). На рис. 2 показаны результаты сравнения численного решения и аналитических формул. Максимальная ошибка составляет < 3 -10-5 и достигается, как легко видеть из графика, на поверхности сферы. (Острый провал графика на некотором расстоянии снаружи сферы вызван особенностями переключения режимов

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

3. Квадруполь, образованный четырьмя круговыми цилиндрами. При достаточной длине цилиндров поле в середине квадруполя должно быть неотличимо от двумерного плоского поля. Рис. 3 показывает результаты сравнения трехмерного поля, рассчитанного с помощью описанного в данной работе алгоритма и двумерного поля, рассчитанного с высокой точностью (< 10-11) независимой двумерной программой. Отклонение между результатами на рис. 3, б нормализовано к потенциалу, приложенному к цилиндрам. Эквилинии на рис. 3, б соответствуют шагу 10-5, так что максимальная относительная ошибка между двумя полями, не превышающая 7 -10-5, достигается внутри "коридора" между цилиндрами и составляет величину < 10 5 для рабочей области квадруполя.

ЗАКЛЮЧЕНИЕ

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

а

б

Рис. 3. Распределение потенциала в средней плоскости квадруполя, образованного четырьмя круговыми стержнями: а) представление распределения потенциала с помощью уровней яркости; б) отклонение рассчитанного по изложенной трехмерной методике распределения потенциала от распределения потенциала, рассчитанного двумерной "эталонной" программой

3

потенциала в непосредственной близости от граничной поверхности (т. е. когда расстояние от пробной точки до поверхности много меньше, чем расстояние между соседними узлами интегрирования) для тестовых задач оказалась равной примерно 10-5. Благодаря гармоническим свойствам интегралов, используемых для вычисления потенциала, это означает, что ошибка вычислений будет гораздо меньше (~1015), когда пробная точка оказывается вдали от граничной поверхности. Хотя

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

ПРИЛОЖЕНИЕ

1. Выражения для эталонных сингулярных функций

Фоо( p, я) =

ф оо( P, Я) =

--

а,8

• № (р, я)+1 я • ЕР(р, я)+^р(р, я)+1 р • ^е(р, я)+ %ре(р, я);

Фю( p, я) = р;

К

ф о(р,я) = - !аЪр„Я„8 + + се .ЕР(р,я)- ар^ • я • ГГ(р>,я)-

2а Ъs а

, 2 2 2 2 ,

Ъс 2 \ - - а р„8 , \ 1 2

2а

■•я2 я)+-

22

2а Ъs

е- р„

я)+ — -р2 •^^я)+' ,

2Ъ аЪ8

Фо1( p, я) = я;

К

+а; р;»^с+с- ^ я)-^рс+я.-р.Ро(р, я у

ф о1(P, я) = -:

2аЪ 8

2 - Ъ 2 2 2 1

а° 2 ^о(р,я)+ „ ,2я„8 -рр(р,я)+ — •я2-рр(р,я)+-я„-^^(р,я);

2Ъ

2а

а,8

Ф2о( P, я) = 4“

К

Ъ2( - 3с2) 3 / \ ъ(- Ъя„ ( - 3с2)+ 2ар„с) 2 / ч

ф2о(р,я) = —VI—'-я3 -*р(р,я)+ у „у , 3——ля -^р(р,я)-

6а 2а

- 2а 2 р„2 + Ъ 2 я„2 ( - 3^2)- 4аЪр„я„с + -2 • ч,ВР(р я) +

2а3

+82,я„(6а2р„ + 6аЪр„я„с -Ъ2я1(- 3с2 )К 3(-Ъя„ + 2ар„с)2 ^ ^р(р я)_

6а ,82

+ 3--P3 -FQ(p,q)+ Ph33P'2!2Г^).FQ(p,q)+-C■q2 • K(p,q) +

3b 3a bs 2a

+ -i7-pq • K(p, q )+ -О- p ■Kip, q )-5aPh - ^ ^^(p, q)-

6a2 6a 6a

- 5abPhqhs2 + 3b2qh2s2c + іЄc q)+ є(є2 - 3a2Ph2s2 ). FFq);

6a 3a bs

Фп( p, q) = pq ;

З-^ q3 -Fp(p,q)-aPh2+ bqhC •q2 • Fp(p,q.+

3a 2a

b2qhs2(3aph + bqhc)+ Зє2(aph + bqhcF Fp(p,q) pз •FQ(p,q)-

bq

•P3

6a 2b2 s2 V"’J/ 3b2

*+ aphc ^2 •FQ(p,q)+ a2phs2(3bqh + aphc) + Зє2(bqh + aphc) FQ(p,qF 2b2 6a2b2s2

+ tT! • p2 •K^’ q F+tV q2 •к(p, q F+■:p■r■ p^Cp q .+-V q^^ q У

3b 3a 6b 6a

a2pls 1 + -1 q¡S2 - є • K(p, q. є( + 3abw2 F ff(p, q);

6a b s 3a b s

Voi( p, q) = ^ ;

K

fÜLl^.p3 • FOip, q )+ a(- aph(l - 3C32 К ^ ) . p1 ,fQ(p, q).

6b3 lb

2 2 2 2

2

+

-1b 2 qh + a 2 ph(l - 3^2 )- 4abPhqhc + є2 • pFQ(p q) +

s2aPh(6b2qh + 6abPhqhc-a2Ph2(l-3c2F+ з(-aPh + 1-qhc)2 • FQ(p q.

6ab3s2

3--q3 •Fp(p, q.+^^brr^-.P’ q K^p2 ^(p* q F+

3a 3ab s lb

wK(p,q)+ -b-q)-5-q + -W •p^Lp,q)-

6b 6b 6b

5abPhqhs2 +3a2p2s2c +1є2c • K(p,q)+ є(є2 -3b2q¡s2 F. FF(p,q).

+

2. Графические материалы

Рис. П1. Нерегулярные компоненты интегрального ядра при малых значениях - :

а — ]Д/р2 + я2 + -2 , б — рД/р2 + я2 +— , в — я/\1р2 + я2 + -2 , г — р2Д/р2 + я2 + -2 , д — ря/4р2 + я2 + -2 , е — я2/л]р2 + я2 +-2 , ж—рр2 + я2 + —2 , з — р2^^р2 + я2 +— , и—ря2Д/р2 + я2 +— ,

к — я V V р2 + я2 + —

Numerical value and theory Relative error of numerical integration

области [о,1]х [од]: а — численные результаты и аналитические результаты; б — относительная

ошибка численного интегрирования

Рис. П3. Численное интегрирование функции /1о — )=Ц"

2 2 2 p + q + £2

=&р &я для прямоугольной

области [-1.9, 2.1 ]х [- 2.1,1.9]: а — численные результаты и аналитические результаты; б — относительная ошибка численного интегрирования

б

а

Relative error of numerical integration

2

Рис. П4. Численное интегрирование функции / — )= ГГ___________р_ф ¿д для прямоугольной об-

20 пу1р2 + я2 +—2

ласти [-1.9,2.1]х [- 2.1,1.9]: а — численные результаты и аналитические результаты; б — относительная ошибка численного интегрирования

б

а

Relative error of numerical integration

Рис П5. Численное интегрирование функции /п (g) = JJ

pq

■\jp 2 + q 2 +e2

б

dp dq для прямоугольной об-

ласти [-1.9,2.1]х[- 2.1,19]: а — численные результаты и аналитические результаты; б — относительная ошибка численного интегрирования

а

1Л? (order=3)

X min=-2,1846914308176 X max=-0,477121254719662 V min=-3,18386755779055 Y max=-1,38804888499853

p/R (order=3)

К min=-2,1846914308176 X max=-0,477121254719662 Y min=-5,93085950315715 Y max=-2,48337046974848

рл2Л? (order=3)

X min=-2,1846914308176 X max=-0,477121254719662 Y min=-8,89384872912636 Y max=-3,S8013166982066

pq/R (order=3)

X min=-2,1846914308176 X max=-0,477121254719662 Y min=-8,75266558482411 Y max=-3,47682146165538

б

а

в

г

Рис. П6. Сходимость процесса численного интегрирования для эталонных нерегулярных функций в зависимости от шага интегрирования: а — ГГ| ¡Д/р2 + я2 + — 2 б —Г Гг _/ I _2 , ,.2 , „2

-JJ|7V p2+q2 +е2 dq, б — Я( p^p2+q2+£l dq ’

JJ^ p 2! jp 2 + q 2 +£2 ^jdp dq , г — JJ^ pq/^p 2 + q 2 +£2 jdp dq

в

Благодарности. Автор благодарит Shimadzu

Research Laboratory (Europe), Ltd. за поддержку этого исследования; М. Судакова за полезные и конструктивные критические обсуждения в процессе работы и за подготовленный им рис. 3, который он разрешил использовать в данной публикации; С.И. Шевченко за идеи, которые явились следствием наших совместных дискуссий, и за ссылки на работы Б. Фрейнкмана, В.Я. Иванова,

за многолетнее дружеское сотрудничество и за предоставленный им в свое время код программы "Топаз"; М.А. Монастырского за неизменное дружелюбие, с которым он хранит математические секреты своих алгоритмов (чем в значительной степени стимулирует разработку своих собственных), и за то влияние, которое на меня оказали его блестящие работы и нестандартные идеи.

СПИСОК ЛИТЕРАТУРЫ

1. Yamamoto K., Nagami K. Difference Method for Electrostatic Field Equation in general Orthogonal Curvilinear Coordinates System // Optik. 1977. V. 48. P. 69.

2. Иванов В.Я. Методы автоматического проектирования в физической электронике (в двух томах). Новосибирск: Труды Института математики СОАН, 198б.

3. Rouse J. Three Dimensional Computer Modeling of Electron Optical Systems // Advances in Optical and Electron Microscopy. Harcourt Publishers Ltd., 1994. V. 13. 400 p.

4.Murata H., Ohye T., Shimoyama H. Improved 3D Boundary Charge Method // Proceedings SPIE. 1996. V. 2858. P. 103-114.

5. Greenfield D., Monastyrski M. Three-Dimensional Electrostatic Field Calculation with Effective Algorithm of Surface Charge Singularities Treatment Based on the Fichera’s Theorem // Nuclear Instruments and Methods in Physics Research, A. 2004.V. 519. P. 82.

6. http://www.electronoptics.com Программа CPO-3D. Описание программы и демонстрационная версия.

7. http://www.simion.com Программа SIMION 3D, версия 7.0. Описание программы и демонстрационная версия.

8. http://www.mebs.co.uk Программы 3D и Go-3D. Краткое описание.

9. Cubric D., Lencova B., Read F.H., Zlamal J. Comparison of FDM, FEM and BEM for Electrostatic Charged Particle Optics // Nuclear Instruments and Methods in Physics Research, A. 1999. V. 427, N 1-2. P. 357-362.

10. Hawkes P.W., Kasper E. Principles of Electron Optics. Vol. 1. 1989, London: Academic Press. 635 p.

11. Фрейнкман Б.Г. Вычисление электростатического поля вблизи заряженной поверхности // ЖТФ. 1979. Т. 49, № 11. С. 246.

12. Фрейнкман Б.Г. Выделение особенности в интегральных уравнениях трехмерного электростатического поля // ЖтФ. 1980. Т. 50, № 2. С.425.

13. Fichera G. The Asymptotic Behavior of Electric Field and Charge Density in the Vicinity of Singular Points of a Conducting Surface // Adv. Math. Sci. 1975. V. 183, N 3. P. 105.

Институт аналитического приборостроения РАН, Санкт-Петербург

Материал поступил в редакцию 7.07.2004.

CALCULATION OF THREE-DIMENSIONAL ELECTROSTATIC FIELDS BY THE BOUNDARY ELEMENTS METHOD WITH REVEALING KERNEL SINGULARITIES NEAR THE ELECTRODE SURFACES

A. S. Berdnikov

Institute for Analytical Instrumentation RAS, Saint-Petersburg

The paper offers a modified boundary elements method to be used in high accuracy (relative error < 10-5) calculations of three-dimensional electrostatic fields for electrodes of almost arbitrary shape. Experiments with analytical test tasks are demonstrated and analyzed.