Федоренко Р.П. Введение в вычислительную физику (1185915), страница 80
Текст из файла (страница 80)
(10) В дальнейшем мы обсудим необходимые свойства функции Д. Вторая операция — решение уравнения Пуассона (7). Уравнение дополняется условиями периодичности — это важное для дальнейшего обстоятельство. Оно означает, что все двумерное пространство заполнено частицами, положения и скорости которых периодически (по х, у) повторяются. Поэтому никаких других сил, кроме сил взаимного электростатического отталкивания (притяжения) частиц, в системе нет. Учитывать можно только взаимодействия частиц, расположенных в одном прямоугольнике в пространстве х, у, сторонами которого являются периоды по х, у соответственно. Действие всех остальных частиц учитывается периодичностью потенциала. 13 — »ззз ем известного закона Ньютона: если частица а действует на частицу Ь с силой Р, то частица Ь действует на частицу а с силой — г".
Следствием этого является отсутствие «самодействия»ч если частица в системе только одна, на нее никакие силы не действуют. В первых схемах метода облаков (или макрочастиц), составленных так, как было описано выше, был обнаружен недостаток: при единственной частице образовывалось поле, действующее на нее, и частица начинала движение, являющееся чисто счетным эффектом.
Расчетная схема не обладала законом «действие равно противодействию», т.е. законом сохранения импульса, что порождало нефизические эффекты. Конечно, формально такие эффекты исчезающе малы при возрастании числа частиц, но реальные расчеты приходится проводить при не столь уж большом числе частиц. И если есть возможность избавиться от такого эффекта, это следует сделать.
Изложим прежде всего способ анализа схемы метода макрочастиц с точки зрения закона Ньютона «действие равно противодействию». Проследим и формально опишем процесс формирования снл в какой-то фиксированный момент времени, Движение мы не рассматриваем, поэтому речь идет о следующей цепочке вычислений. Рассмотрим совокупность частиц, имеющих заряды р» и координаты (Х, У ) (к = 1, 2, ..., К). Первая операция — расчет плотности р, в узлах эйлеровой сетки. Отвлечемся от конкретного способа, связанного с той или иной формулой распределения заряда в облаке.
Ведь влюбом случае заряд к-й частицы д» распределяется по узлам сетки в соответствии с некоторой функцией Д(1, л»; Х, У) таким образом, что вклад 1»-й частицы в заряд в точке (1, т) есть р)»~ = д» Д(1, »л; Х», У»), а полный заряд в точке (1, »и) есть, очевидно, лгивлиженные методы вычислительной визики 1ч. и ззб Если бы мы рассматривали уравнение Пуассона с какими-то другими условиями, например с заданнымн значениями р на границе, это означало бы наличие внешних сил.
В такой системе импульс не обязан сохраняться. Решение разностного уравнения Пуассона может быть выражено через разностную функцию Грина: 'рь = Х б(1 т' 1 у) ре у. (11) В у Это очевидный факт. Функция б(1, т; 1, у) — просто матрица, обратная к матрице разностного оператора Лапласа. Здесь и в дальнейшем мы не будем аккуратно выписывать пределов суммирования: оно ведется по одной ячейке периодичности (или, если угодно, на торе).
Мы обсудим свойства б ниже, пока же оставим запись преобразования р в р в самой общей форме. Очевидно, потенциал всей системы частиц есть сумма потенциалов, порожденных каждой частицей: (12) р~в> = ~б(1,т;1,у)р~; еу Третья операция — вычисление сил (градиента р) в узлах сетки. Ограничимся только одной компонентой силы, так как для другой все будет точно так же, Итак, пусть сила у, есть — рв(ху, у.). В разностной реализации любая аппроксимация может быть записана в общем виде: г(1, у; 1, т) ~реп. (13) См Четвертая операция — вычисление силы, действующей на частицу с зарядом о, расположенную-в точке (х, у). Это достигается процедурой интерполяции сеточной функции у, в точку (х, у), которую можно записать в виде 5(х, у) = д ~; 1(х, у; г', у) У~ ..
(14) су Если нам нужно подсчитать силу, действующую на г-ю частицу со стороны й-й, необходимо провести следующую цепочку преобразований: 5(Х„, У„; Х,, Ув) = гУ,ц ~ У(Х„, У„; 1, У) х су х ~~, г(1, у; 1, т) ~ б(1, т; р, з) Ц(р, з; Х, 1' ) = См л, ю = д„д у(Х„У„; ') Г( °; *) б(*; «) д(в; Х, У„). (15) й 241 ИРизлиженное ннтегРиРОВАние уРАВнення ВлАЕОВА 337 Здесь, », в — символы «немых» индексов, по которым произведена свертка.
Сила 5(Х„, У„Хы Уь) есть элемент (г, к) матрицы, являющейся произведением некоторых других матриц. Введем С и Р— пространства функций, определенных в точках (х, у) и'(У, ун) соответственно. Алгоритм интегрирования уравнения Власова определяется прн конкретизации следующих операторов: Д: С- УУ вЂ” оператор «раздачи заряда» (из точки (х, у) в узлы сетки); О: Ь- УУ вЂ” разностный оператор Грина, обратный к оператору Лапласа с периодическими условиями; Р: Р— УУ вЂ” оператор вычисления первой разностной производной; У: Р- С вЂ” оператор интерполяции силы с узлов сетки на произвольную точку.
Таким образом, 5 есть оператор типа С- С, а действие й-й частицы на г-ю есть (множнтель д,д опускаем) Ю(Х„, У,; Х, У ) = (УРОД), Тогда действие г-й частицы на к-ю есть Л(Х,, У,; Х„, У„) = (УРОД) „„= (УРОД)', = (Д О'Р'У')„,. Для того чтобы в системе частиц был справедлив закон «действне равно противодействию» и, как следствие, сохранялся импульс, нужно обеспечить соотношение (16) Д'О'РР = — УРОД. Используем следующие почти очевидные свойства операторов.
а) О'= О (следствне самосопряженности оператора Лапласа в классе периодических функций). Ббльшая часть разностных аппроксимаций оператора на равномерной сетке автоматически наследует это свойство, хотя при желании его можно и нарушить, не потеряв аппроксимации, б) Р' = — Р (следствие кососимметричности оператора д/дх). Аппроксимация типа (9), как нетрудно проверить, сохраняет это свойство; аппроксимация (д//дх),. = (/, — /,,)/УР его нарушает: (д//дх); = -(/,.„— /,.)Н.
в) РО = ОР (следствие перестановочности операторов д/дх и А). Это свойство наследуется использованными выше разностными аналогами операторов. Таким образом, нужно обеспечить равенство (17) Д'РОУ' = УРОД. Оно выполняется, если Д = У". В схеме с гарантированным сохранением импульса из операторов Д, У только один строится независи- 388 ПРИВЛИЖЕННЪ|Е МЕТОДЫ ВЫЧИСЛИТЕЛЬНОД ФИЗИКИ 1ч. н мо, например оператор интерполяции 1. Оператор раздачи заряда Д после выбора 1 определяется автоматически.
Поясним смысл сказанного. Пусть оператор 1(х, у; ю', 1) реализован как кусочно-билинейный. Это означает, что для вычисления 7(х, у).=,'Р1(х, у; 1,1) /,. (1й) |,/ нужно найти индексы 1, т, при которых х м [1Ь, 1Ь + 1), у |= [тЬ, тЬ + Ь), и вычислить соответствующие коэффициенты интерполяции н, 8 (а, р О, 1). Тогда Л(~ У) = иао/с + тгьв А~нс + иа ~д +|+ м|,| А|н, | (1о) где ив Ф= Ь з(1Ь+ Ь вЂ” хНтЬ+ Ь вЂ” у) и т.д.
Элеменг «матрицы» 1 (1, 1; х, у), очевидно, вычисляется так: по точке (х, у) нужно найти индексы 1, т и коэффициенты интерполяции, после чего определюотся элементы матрицы 1'. 1'(1+ а, т + Р; х, у) = и, 8, а, Р = О, 1. Все остальные элементы 1 равны нулю. Итак, операция раздачи заряда аы находящегося в точке (Х, г ) состоит в следующем.
По значениям Хы Уь находится ячейка сетки (1, т), в которой эта точка находится, вычисляются соответствующие коэффициенты интерполяции и 8, и доли заряда а«|г раздаются в узлы (1+ а, т + Р): Р|+а,а»а' Р|»а,е»В+ |4 ар Если бы использовалась более точная интерполяция, мы столкнулись бы с таким «нефизичным» явлением, кзк раздача отрицательной доли заряда в какой-то узел. Видимо, поэтому на практике ограничиваются простейшим кусочно-билинейным оператором интерполяции, когда все и 8 |в О и ~ч' |г 8— - 1. Таким образом конструируются схемы интегрирования, обеспечивающие непрерывность плотности заряда р, по положениям частиц (Хы г'„) и сохранение импульса.
Быстрое дискретное преобразование Фурье. Имеется некоторая сеточная функция 1(Ь) (Ь = О, 1, ..., Ф вЂ” 1). Она может быть представлена разложением в сумму Фурье: У-| 1(Ь)= ' 7 1(л) ., = '", Ь=О,1,...,Ы 1, (20) и 0 з 241 ПРИБЛИЖЕННОЕ ИНТЕГРИРОВАНИЕ РРАВНЕНИЯ ЕЛАСОБА ззо Коэффициенты Фурье А(п) вычисляются по аналогичным формулам; И-1 А(п) =д~~ У(/с) 1Р "", и =О, 1, ..., Дà — 1.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
