Файл: Вычислительные методы в физике плазмы..pdf

ПРЕДИСЛОВИЕ К АНГЛИЙСКОМУ ИЗДАНИЮ

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

Метод частиц является гибким, но требует много машинного времени. Одномерный вариант этого метода дается.» в главе Доусона, двумерные варианты описаны в главах Хокни, Морза

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

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

в частности, справедливо в задачах, в которых определяющую роль играют столкновения, как показано в главах Робертса и Пот­ тера, а также Киллина и Маркса. В общем случае метод расчета, который используется в непрерывных моделях, представляет собой метод конечных разностей для уравнений типа гидродинамических,

В бесстолкновительном случае (случай Власова) в одномерных задачах были использованы также два специальных метода, опи­ санные в главах Армстронга с соавторами и Бэрка и Робертса.

Дальнейшее совершенствование численных методов и вычисли­ тельных машин должно привести к еще более правдоподобному моделированию плазмы, включая трехмерные расчеты. К сожале­ нию, применение описанных ниже численных*методов к близким по характеру задачам динамики звезд не могло быть включено из-за ограниченного объема книги.

Январь 1970

Берни Олдер Сидни Фернбах Мануэлъ Ротенберг

ГЛАВА 1

ЭЛЕКТРОСТАТИЧЕСКАЯ МОДЕЛЬ ПЛОСКИХ ЛИСТОВ ДЛЯ ПЛАЗМЫ И ЕЕ МОДИФИКАЦИЯ ДЛЯ ЧАСТИЦ КОНЕЧНОГО РАЗМЕРА

Дж. Доусон*

§ 1. В ведение

Модель плоских листов — одна из наиболее старых и гибких одномерных моделей для численного исследования плазмы. Впер­ вые она была применена для исследования электростатических эффектов, плазменных колебаний, кинетики одномерной плазмы и двухпотоковой неустойчивости [1—5]. В этой главе мы обсудим электростатическую модель плоских листов. Основное внимание будет уделено методам решения этих проблем на ЭВМ и ряду трудностей, которые при этом встречаются (например, шуму). Однако реальные программы приводиться не будут.

Позднее модель плоских листов была модифицирована для того, чтобы допустить наряду с перпендикулярным движением движе­ ние в плоскости самих листов [6—8]. Эта модификация позволяет включить магнитные эффекты и эффект соударений в приближении Фоккера — Планка (используя метод Монте-Карло). На основе таких модификаций были изучены циклотронные волны [61, эффект излучения 18] и неустойчивости типа Харриса (см. [9]), связанные с моноэнергетическими распределениями по скоростям. Хотя эти последние модификации весьма интересны и разнообразны, мы не будем их здесь рассматривать, поскольку это увело бы нас слишком далеко.

§ 2. Э л е к т р о с т а т и ч ес к а я модель п л о с к и х л и с т о в

Имеются две электростатические модели плоских листов: однокомпонентная и двухкомпонентная. В однокомпонентной модели (Доусон [2]) рассматривается плазма, состоящая из боль­ шого числа тождественных заряженных плоских листов, поме­ щенных в неподвижный однородный нейтрализующий фон (фиг. 1). Предполагается, что все плоские листы всегда перпендикулярны какой-то оси X, но могут свободно двигаться вдоль нее. Считается, что они свободно проходят друг через друга.

Двухкомпонентная модель [4, 5] включает как положительно, так и отрицательно заряженные листы. Снова предполагается,

* John М . Dawson, Princeton University, Plasma Physics Laboratory, Princeton, New Jersey.

Легко видеть, что (1) есть уравнение гармонического осцилля­ тора, так что

(2) (3)

Каждый лист осциллирует с плазменной частотой независимо от всех других.

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

положения равновесия.

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

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

имогут осциллировать только вдоль этой линии центров.

Спомощью этой модели можно продемонстрировать ряд свойств

одномерной плазмы. Например, если первый маятник отвести в сторону и освободить так, чтобы он смог ударить второй, то он передаст свою скорость второму, второй — третьему и т. д. Таким образом, импульс движется через маятники. Когда какой-то маят­ ник ударяет соседа, он теряет свою скорость, но не свое смещение. Таким образом, импульс оставляет маятник после себя в смещен­ ном состоянии, и маятники начинают колебаться. Это эквивалент­ но возбуждению колебаний плазмы при движении быстрого листа через плазму.

2. Численный метод описания движений однокомпонентной одномерной плазмы

Будем искать решения уравнения (4) для листов, корректируя орбиту для пересечений с соседними листами. Когда лист пересе­ кает соседний, ускорение претерпевает скачок,

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

Допустим, что частицы пронумерованы в соответствии с их положением вдоль оси х. Если не происходит пересечение, то реше­

Код вычисляет эти непересекающиеся положения и проверяет, не произошло ли какое-либо пересечение, т. е. он разыскивает ситуации, когда

(Напомним, что частицы пронумерованы вдоль оси х.) Если он находит, что такое пересечение произошло, он вычисляет время

Ф и г. 2. Схема пересечения.

пересечения, определяя сначала хорды двух орбитальных кривых в плоскости X, t, а затем вычисляя времена пересечений для них (фиг. 2). Первое приближение для времени пересечения £с1 запи­ сывается в виде

Далее код вычисляет, исходя из непересекающихся орбит, положения этих двух частиц в предсказанный момент пересечения. Затем он вычисляет хорды, проведенные через эти точки в момен­ ты t и tc1 , и определяет их время пересечения tc2. Это время исполь­ зуется как правильное время пересечения. Вклад от более ранних пересечений (если таковые были) в течение Аt не включается. Если две частицы пересеклись до момента t из-за поправок к их орбитам от пересечений в предшествующий временной интервал, то tc2 будет отрицательным и метод найдет поправку. Орбиты частиц корректируются путем добавления постоянного ускорения + ©|б к их движению на протяжении остальной части временного шага, причем знак зависит от порядка пересечения.

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

( 10)

где Лt — временной шаг. Эта скорость довольно велика, кроме того, точность имеет только первый порядок относительно At. Благодаря включению первой поправки для пересечения закон сохранения энергии выполняется с большей точностью, ~ At2, и это примерно на три порядка лучше для <врДі = 0,05 и rikD~ 10. Введя вторую по­ правку для времени пересечения, мы получили дальнейшее улучше­ ние выполнения закона сохранения энергии на два порядка для (OpAt = 0,05 и rikD = 10. Увеличение рабочего времени из-за этой последней поправки оказалось незначительным (~1 %) и с избытком окупалось возможностью использовать большие временные шаги при заданной точности.

Наряду с проверкой сохранения энергии ранний вариант кода без второй поправки для времени пересечения проверялся на обратимость во времени. Движенйе системы из 9 листов было обращено, и мы обнаружили повторение траектории системы с точ­ ностью до ІО-3 (все орбиты имели эту точность) в течение 6 колеба­ ний (сdpt — 36). В этом случае на дебаевской длине находилось примерно 2,5 частицы.

Если частицы пронумерованы в том порядке, в котором они встречаются на оси х, то при проверке на пересечения нам нужно

при котором частицы еще могут сталкиваться на протяжении вре­ менного шага. Это расстояние определяется формулой

Д^макс ( 0 — — макс “I“ ^ (0 ^ ^ * (И)