Спросить
Войти
УДК 519.6:51-73:004.942
А.В. Павельчук
ИСПОЛЬЗОВАНИЕ МЕТОДА УСТАНОВЛЕНИЯ ДЛЯ ЧИСЛЕННОЙ РЕАЛИЗАЦИИ МАТЕМАТИЧЕСКОЙ МОДЕЛИ ЭЛЕКТРОСТАТИЧЕСКОГО ПОЛЯ
В работе представлена конечно-разностная реализация модели расчета потенциала электрического поля, созданного заданным распределением электрических зарядов. Вычислительная схема решения уравнения Пуассона построена на основе метода установления и метода прогонки для решения полученной системы сеточных уравнений с учетом цилиндрической симметрии многомерной задачи. Проведен расчет параметров вычислительного процесса с использованием практической оценки сходимости схемы.
ESTABLISHING METHOD APPLICATION FOR NUMERICAL IMPLEMENTATION OF ELECTROSTATIC FIELD MATHEMATICAL MODEL
The paper deals with a finite-difference implementation of calculation model for electric field potential generated by continuous charge distribution. A computational scheme of the Poisson&s equation was based on the establishing method together with the sweep method for solving the obtained system of finite-difference equations subject to cylindrical symmetry of multidimensional problem. The calculation of control parameters was performed with the use ofpracti-cal estimation of scheme convergence.
Введение
Уравнение Пуассона лежит в основе линейной модели математической физики, описывающей пространственное распределение основных характеристик электростатического поля - напряженности и потенциала, создаваемых системой электрических зарядов [1]. Как известно, в случае произвольной геометрии среды, сложных граничных условий или специально заданной функции пространственной плотности распределения зарядов аналитические методы работают только для ограниченного класса задач. В практике инженерных и научных вычислений широкое распространение получили методы численного анализа [2-4]. Существует множество подходов для построения численного решения многомерных уравнений эллиптического типа, в числе которых можно выделить методы установления [3-4]. Основная идея методов этой группы заключается в преобразовании стационарной задачи в нестационарную, причем рассматривается релаксация системы к стационарному состоянию, с заменой шага по времени номером итерации. Конечно-разностная аппроксимация при этом проводится на базе уже разработанных для решения нестационарных задач методов расщепления [5-6].
Для осуществления итерационных процессов используются следующие разностные схемы: явная разностная схема, схема расщепления, схема переменных направлений, схема предикат-корректор и др. [2]. Проводя сравнительный анализ методов, можно отметить, что при использовании явной разностной схемы имеем второй порядок аппроксимации по пространственным переменным и первый - по времени. При этом схема является условно устойчивой и требуемое количество итераций
составляет величину порядка -1 . Схема расщепления будет абсолютно устойчивой, количество итеh
раций составит — , схема имеет также второй порядок аппроксимации по пространственным переh
менным и первый - по времени. При использовании схемы переменных направлений и схемы предиктор-корректор достигается второй порядок аппроксимации не только по пространственным переменным, но и по времени; схемы абсолютно устойчивы и количество необходимых итераций - есть 1
величина порядка — .
Задача решения уравнения Пуассона является важной составляющей более общей проблемы моделирования динамического процесса зарядки диэлектрических материалов при воздействии электронного облучения. В работах [7-8] представлены основные концепции этой модели и варианты численных схем ее реализации, в частности в работе [9] построен алгоритм решения нестационарной задачи для определения пространственного распределения объемной плотности зарядов на основе метода расщепления. В связи с чем цель данной работы - конструирование эффективной конечно-разностной вычислительной схемы, основанной на методе переменных направлений, для решения модельной многомерной задачи о расчете потенциала при заданном пространственном распределении объемной плотности системы зарядов.
Постановка задачи
Сформулируем математическую постановку стационарной задачи численного решения нахождения потенциала и напряженности электрического поля, создаваемых заданным пространственным распределением электрических зарядов. Математическая модель формулируется в виде уравнения Пуассона и выражения для расчета напряженности поля:
ДФ = --Р-,
ее0 (1)
Е = - grad ф,
где р - объемная плотность распределения зарядов, Кл/м 3; ф - распределение потенциала, В; Е - напряженность поля, В/м; в - удельная диэлектрическая проницаемость среды;в0 - диэлектрическая постоянная, Кл/(В-м).
Примем во внимание цилиндрическую симметрию пространственной задачи с соответствующей записью оператора Лапласа:
д2ф 1 дф д 2ф Дф = —2 +—- + —-.
дг г дг дz
Пусть для (1) (гг )е G , Г - граница области G, G = G + Г = {0 < г < R, 0< z < Z} - прямоугольник. Будем рассматривать смешанную краевую задачу, дополнив уравнение (1) граничными условиями на Г: дф
= 0, ф г=*, = 0. (2)
г=0, г= 2
Конструирование вычислительной схемы для уравнения Пуассона
Построение численных решений стационарных задач математической физики, описывающих равновесные состояния, основано на рассмотрении релаксации динамического процесса к установившемуся режиму. Применение такого подхода для целого класса прикладных задач дает эффективный инструмент, который позволяет проводить вычисления проще, чем прямой расчет равновесного состояния. Использование схемы переменных направлений в концепции метода установления для решений уравнений с частными производными эллиптического типа характеризуется быстрой сходимостью и абсолютной устойчивостью решения.
Метод установления заключается во введении в исходное уравнение фиктивной производной по времени:
дф Л Р
— = Дф +-1— . дt 880
при этом функция ф( г г ) становится функцией трех аргументов: ф( г г ( ) .
Уравнение вида (3) относится к уравнениям с частными производными параболического типа, в силу независимости от времени граничных условий с течением времени производная по времени будет стремиться к нулю, а решение нестационарной задачи ф(гг ( ) - к решению исходной стационарной задачи ф(гр).
Введем в рассмотрение сетку узлов с шагами К,, К2, т соответственно по переменным гр t :
= {г = &К, ; = О, I, г. = jh2, j = 0,3 }к = к т, к = 0,1, 2,...} и на этой сетке будем аппроксимировать
уравнение Пуассона задачи (1) на верхнем временном слое t + = (к + 1)т:
к+12 к+ 12 1 ф+1,2 - ф ,, /2 + Ф
„, к к , к +1 ф. - ф ,,, + ф 2
к +1
к + X
г,-1 +1 ф,+62
к + К
к +К
+1 ^ к+1 , +1 ! +X — фк+У2
л к+1 /■* к+1 , к+1 л /1 /, +1 ф;,;+1 - ф и + ф;,+1 ф;+1,, - Ч>,-1. + Р;,„
Рис. 1. Вычислительный шаблон для разностной схемы, построенной на основе метода переменных направлений.
к +12 ф;+и
к+>2 ;-и
к + К
+ <Р,„
к , т Г к п. к , к
= ф& + гК7 ,+1 - ф ; + ^ +
Р;,;т 2880
На втором полушаге разностная задача представляется в виде (вторая подсхема):
к+1 ф++1
т + фк+1 1 + 4" К22 + фк+1 т
- 2к2 + ф&,; -1 2К22
к + V
&к+12 к+12
к + К
р ■ . т
к + к
Для согласования порядка аппроксимации производной по координатным направлениям в исходном уравнении и граничных условиях по границам г=0, z=0, введем вдоль этих границ фиктивные
узлы ф0 . для 1 = 1,М и ф.. 0 для & = 1,Ы и воспользуемся подходом, описанным еще в раннем обзоре
[6]. Аппроксимируем граничные условия со вторым порядком точности по координатам:
Ф2, , " ф
= 0 для 1 = 1М и
ф, 2 - Ф,,: 2К
= 0 для & = 1,N .
Далее вычислительный процесс организуем по полученным формулам (6)-(7), исключая с помощью (8), значения потенциала в фиктивных узлах для первых уравнений полученных систем. Согласно [6], в такой концепции заложен глубокий смысл: если при решении задач математической физики разностными методами граничные условия исключены из рассмотрения, то для формирования вычислительного алгоритма уже не приходится специально заботиться об удовлетворении граничных условий, поскольку они автоматически учтены в модифицированных разностных уравнениях.
Кроме того, учтем граничные условия вдоль границ r=R, z=Z в последних уравнениях сформированных систем: ф = 0 для 1 = 1,М и ф. м = 0 для & = 1,N .
Каждая из подсхем схемы переменных направлений решается с помощью метода матричной факторизации (прогонки).
Коэффициенты, соответствующие уравнению вида а1 фк
X X , , X XX
а =--г--, Ь. = 1- —, с. =--- +-,
& 2% 4г И & И & 2% 4%
X ф!/+1 - ф I + фк1 -1 Р., 1 _
к+1 "¡У&+1, /
+ Ь ф*+ + с.ф^1, = dk., имеют вид:
, 1 -1, 1 , 1
а. = с. =--1 1 2И2
X к + 1
Ь=1+ИТ "к&& =
к + К
к +12 ^ к +12 к Xф¡+^,/2 - ф &,/2 + ф1,1
На каждом дробном шаге достигается частичная аппроксимация, полная аппроксимация достигается на последнем дробном шаге. Для обеих подсхем выполняется достаточное условие сходимости прогонки [3]:
а \\- + \\с\\ = 2И22 +
=И2"
Расчет итераций (переход од одного фиктивного временного слоя к другому) происходит по схеме (6)-(7), пока итерационный процесс не сойдется с заданной точностью 5:
ф - ф < О .
Исследование устойчивости схемы
Исследование устойчивости схемы переменных направлений, аппроксимирующей уравнение
Пуассона (1), проведем с помощью метода гармонического анализа. Для этого отбросим член ,
который не оказывает влияния на устойчивость разностной схемы, и представим решение в виде гармоники [2]: ф^- = ак ехр пг1 + А,^ ^ . Тогда для
первой подсхемы, получим:
ak+12 а*= а k+12
— ( exp(/l П - 2 + exp(-zl Л)) + —^ ( exp(/X Л) + exp(-zl Л))
T77 (eXP(&^ mh2) - 2 + eXP(-?&^ mh2) ) 2h2
Откуда, учитывая, что
ехр(/А, ) - ехр(-/А, ) = ^(А, ) - ^(-А, ) = 0 и ехрО&А тк^) - ехр(-7&А ,^2) =
= cos(Amh2) - cos(-Amh2) = 0 , получим
В последнем выражении: exp(/^ nh1) - 2 + exp(-/^ nh1) = -4 sin21 —21
Аналогично P2 = 4sin | m 2 |. Параметры: G =—-, G == - Pi.
Проводя аналогичные рассуждения для второй подсхемы, получим:
а». i а i / — а i /
k+1 k +12 k +12
— (exp(/l „h1) - 2 + exp(-/l „h1)) + —— (exp(/l„h1) + exp(-/l „h1)) 2h1 4r h
+а.
—y (exp(/lmh2)- 2 + exp(-/lmh2))
откуда
аk+1 1 - G1P
ак+X 1 + g2p2
. Составляем отношение амплитуд гармоник:
а+1 а+1 а+к 1 - G1 P 1 - g2p2
ак а k+у2 а k 1+g2p2 1+G1P
Последнее неравенство означает, что разностная схема переменных направлений (6)-(7) абсолютно устойчива, что обеспечивает возможность произвольного выбора значения шага итерации. При этом порядок точности аппроксимации по т более высокий (второй), что позволяет использовать при решении более грубый шаг итерации, ускоряя сходимость итерационного процесса и уменьшение количества итераций.
Численный пример
Продемонстрируем решение задачи (1)-(2) на модельном примере, для которого инициализируем следующий набор параметров: R = 5 усл.ед., 2 = 5 усл.ед., \\= ^ = 0.1; т = 0.5 ; р = -1 усл.ед., 880 = 1 усл.ед. Ввиду того, что независимо от используемого критерия окончания вычислений, во избежание зацикливания следует поставить ограничение на количество выполняемых итераций. Выберем максимальное число итераций ктах = 100 , с учетом скорости сходимости метода.
Представим практическую оценку сходимости вычислительной схемы (6)-(7). Используем
указания по выбору оптимального значения параметра т , предложенные в [4]: т
sin ()h )
. Однако
рекомендации являются ориентировочными и требуется проверка скорости сходимости при варьирования шага итерации. Приведем значение погрешности 3 для последнего приближения, зафиксировав максимальное значение числа итераций на значении кт=100. Примем также, что для каждой подсхемы шаги по координатам совпадают /? И2. параметры (за исключением координатных шагов и шага по времени установим такими же, как и в примере). Будем варьировать значения шага по времени т в окрестности рекомендуемого оптимального значения т (при
/7а=/72=ОЛ тор1 « 0.06 ).
Результаты оценки погрешности (с учетом различных способов вычисления нормы) представлены в таблице.
Оценка погрешности решения на последней итерации
h=h2 I Ы+1 - ф1 (норма-максимум) Ы+1 - ф1 (Евклидова норма)
Анализируя полученные данные, можно заключить, что наиболее оптимальным представляется выбор шага итерации т ( 0.6 .
Заключение
Таким образом, построена вычислительная схема для решения смешанной краевой задачи для уравнения Пуассона с учетом цилиндрической симметрии многомерной задачи на основе метода переменных направлений. Для сконструированной схемы (6)-(7) были исследованы устойчивость и сходимость. Показано, что схема является абсолютно устойчивой, имеет второй порядок точности аппроксимации по т, позволяющий использовать при решении более грубый шаг итерации, что способствует ускорению сходимости итерационного процесса и уменьшению количества итераций. Введенную в рассмотрение вычислительную схему можно использовать для реализации математических моделей сложных процессов, описывающих пространственное распределение основных характеристик электростатического поля - напряженности и потенциала. В частности, для решения прикладной задачи о моделирования динамической зарядки диэлектриков, индуцированной электронным облучением.
Z, усл.ед.
Рис. 2. Пространственное распределение потенциала для модельной задачи.
