Fedorenko-RP-Vvedenie-v-vychislitelnuyu-fiziku (810773), страница 81
Текст из файла (страница 81)
пгвглижзнныв методы зычислитвльноИ «изнки (ч. и Зз« Имея И," ' в узлах сетки, можно определить в узлах значения +1 е+! „л«! «~ы (Е )«+! гы ~-ь (Е )п«1 и»ы ~ ~ ~ (9) «П, 2Ь ' ( г)~ «2Л Опыт показал, что силу, действующую на к-е облако, нужно вычислять достаточно аккуратно. Например, с помощью интерполяции можно определить функцию Е(х, у). При вычислении правых частей уравнений (8) интегрируется функция 9(х, у) Е(х, у) по облаку с учетом распределения заряда в нем. Это достаточно сложная процедура, но прн упрощенном вычислении правых частей (8) проявляются нежелательные счетные эффекты.
Анализ привел к пониманию вызывающих их причин и к упрощению техники конструирования расчетных схем, свободных от таких дефектов. Отметим еще один момент, связанный с назначением для расчетов масс и зарядов макрочастиц. Прн моделировании явлений в плазме каждая макрочастица представляет 1О'« -» 1О'~ реальных частиц, заряды и массы которых — хорошо известные физические константы. Какими же брать прн моделировании массы и заряды макрочастнц? Есгесгвенный ответ «их нужно брать в 10г««- 1Ом раз ббльшими реальных» оправдан еще и тем, что, как нетрудно проверить, в такой «макроплазме» значения характерных величин (периодов электронных и ионных колебаний) оказываются теми же, что и в реальной плазме. Расчетная схема с законом сохранения импульса. Выше была описана схема, в которой преодолен один самый грубый недостаток, — флуктуации плотности, появляющиеся при примитивной технике расчета плотности в счетном узле сетки.
Следующим шагом было построение схем с более тонким свойством — сохранением импульса. Поясним суть проблемы. Мы имеем дело как бы с набором твердых частиц, обладающих в некотором смысле пространственной структурой. Но эта структура сказывается только на формулах вычисления сил, действующих между частицами. Реализуется следующая цепочка зависимостей. Положения частиц определяют плотность заряда в узлах сетки; заряд определяет потенциал э в узлах; потенциал определяет силу, действующую на частицы, находящиеся в каких-то точках пространства; в данной схеме центр частицы движется по тем же законам, что и «точка».
Итак, речь идет о системе точек, между которыми действуют силы, вызывающие их движение. А плотность заряда, потенциал— промежуточные объекты, облегчающие вычисление снл по сравне-. нию с физической моделью плазмы, как совокупности заряженных,' частиц с" обычным электростатическим потенциалом. В этой «чис-'1 той» модели взаимодейсгвующих частиц действуют известные зако-* ны, в частности закон сохранения импульса, являющийся следстви- 3 341 пгивлижзннов интгчеиговлннв гг»знания вллсовА 385 к р, „= Х ч» 0(1, т; Х,, У,).
(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).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
