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

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

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

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

всоответствии с положением частиц. Если координатами частицы, используемыми при решении уравнения Пуассона, являются, допустим, величины при t = 0, то хранимые с ними скорости опре­

делены фактически при t — —Д£/2, где Д£ — временной шаг, а скорости и координаты изменяются на временном шаге по сле­ дующей простой схеме:

причем формулы (8) так естественны для чисто электростатических проблем, что вопросов почти не возникает. Однако допустим, что здесь имелось бы также постоянное магнитное поле В, так что Е 0 заменилось бы на Е0 + v0 X В. В этом случае скорость ѵ0 непо­ средственно не известна, а должна определяться по схеме, подоб­ ной схеме Рунге — Кутта, справедливой до второго порядка по At и вытекающей из подстановки оценки первого порядка v0

ж V—1/2 + Е0 (Дг/2).

Можно пойти и дальше, требуя, чтобы разностная схема была обратимой, т. е. ѵ0 = х/2 (ѵ+і/2 у ѵ_і/2), и давала бы еще некоторое удобное, явное алгебраическое выражение для ѵ_щ/2 вместо первого из выражений (8). В плоскости, перпендикулярной вектору В,

это дает следующий результат:

Если поле В не постоянно, а должно определяться самосогласо­ ванно по току в плазме, то положение значительно ухудшается. Для такой проблемы Шонк и Морз [10] применили неявную, но необратимую схему; в обратимых схемах трудно придумать быст­ рую обработку полных уравнений электромагнитной волны.

Необходимость центрирования по времени связана с тем, что перемещение частицы должно быть центрированным по времени

ЬО

по крайней мере до величины порядка (Д£)2, а обратимые схемы сами достаточно хороши, чтобы использовать их всегда, когда можно. Это подтвердилось двумерными моделированиями электро­ статических пучковых неустойчивостей, подобных рассмотренным здесь, но при наличии постоянного магнитного поля, направлен­ ного перпендикулярно скорости сталкивающихся потоков. Соот­ ветствующие расчеты, выполненные Нильсоном и Морзом, про­ демонстрировали постоянное возрастание полной энергии в слу­ чаях, когда член ѵ X В не центрировался должным образом. (Напоминаем, что впервые Оскар Бунеман настойчиво потребовал выполнения обратимости при моделировании плазмы с помощью частиц.)

На фиг. 6 представлены: а) начальное распределение в фазовом пространстве (х, ѵх) случайной подгруппы модельных частиц для двух одинаковых теплых пучковых распределений с L = 100 KD и б) функция /о (ѵх).

На фиг. 7 представлены результаты: а) одномерного, 1-D, б) дву­ мерного, 2-D, и в) трехмерного, 'VD, сопоставимых расчетов этого

случая, который был рассмотрен в l-D-случае с большим L Морзом и Нильсоном [12,13] и продемонстрировал образование, смеши­ вание и длительное существование BGK-мод, видных на фиг. 7. Каждый столбец на фиг. 7 содержит два графика фазового про­ странства (X, ѵх) [первый относится к моменту насыщения энергии поля, а второй — к концу интервала в 10 плазменных периодов

Координата х, t- ^ ,0

—4/71----------- —I 4/7 г -

100 О 100 о

Координата х, t = 10,0

-4/7

Ф и г . 7. Диагностики сопоставимых машинных расчетов для одного (1-D),

двух (2-D) и трех (3-D) измерений.

(не 10 Юр1)], функцию / (ѵх) в конце интервала и график е (t) — изменение во времени полной энергии электростатического поля в единице объема.

Во всех трех расчетах энергия поля сначала плавно нарастала до теплового уровня, а затем резко нарастала до уровня насыще­ ния, причем сперва экспоненциально со скоростью нарастания самой быстрой линейной моды и с большей частью энергии поля в соответствующей интенсивности | Еь |2. (Спектр энергии рас­ считан, но не приведен здесь.) Такая концентрация энергии в наи­ более сильно нарастающих модах является естественным следст­ вием экспоненциального нарастания интенсивности этих мод от малого, однородного уровня спектра начальных шумов [13].

Многочисленные проверки ошибок, которые были выполнены [12, 13], показывают, что число модельных частиц должно быть

достаточно большим, чтобы уменьшить начальный тепловой шум до уровня гораздо ниже уровня насыщения энергии е (t) и соот­ ветственно получить превышение энергии е, т. е. коэффициентов I E k I2 , мод со средними длинами волн и с наибольшими инкре­ ментами над спектральной энергией теплового шума, которая про­ порциональна 1 /к2. Если эти критерии не выьолняются, то правиль­ ные бесстолкновительные результаты не накладываются просто на уровень тепловых шумов, а вместо этого разрушаются шумами.

Высота первоначального теплового плато в (t) в нашем случае увеличивалась при переходе от 1-D к 3-D из-за возрастающей необ­ ходимости экономии числа модельных частиц; в l-D-расчете исполь­ зовалось 20 000 частиц, что было более чем достаточно; в 2-D- расчете использовалось 80 000; а в З-Б-расчете использовалось 332 750 частиц, и этого хватало в обрез. Необходимое число частиц можно было бы снизить, уменьшая L , но ввиду ограничения на число неустойчивых мод используемая здесь длина L = 100 уже достаточно мала. При счете использовался временной шагг равный 0,04 плазменного периода, как и в работах Морза и Ниль­ сона [12, 13], по причинам, изложенным там при анализе ошибок.

Числа ячеек в 1-D-, 2-D- и З-Б-сетках были соответственно 64, 64 X 64 и 32 X 32 X 32. Читатель, который хочет заняться модели­ рованием, должен иметь в виду, что нужно очень осторожно выби­ рать физическую длину L и число ячеек в L. Если целью модели­ рования является изучение турбулентности плазмы, то длина L должна бытъ достаточно большой, чтобы содержать по крайней мере несколько длин волн наиболее неустойчивой в линейном при­ ближении моды. Нелинейное поведение системы, в которой укла­ дывается одна или две длины, может быть интересным, но оно мало что говорит о турбулентности. Использование большего числа ячеек увеличивает плотность фурье-мод в ^-пространстве, которая важна в некоторых случаях турбулентности, и создает возмож­ ность для рассмотрения более коротковолновых мод. В связи с этим недавние исследования показали, что коротковолновые моды с относительно малыми энергиями поля, | Е^ |2, оказываются гораздо более важными при описании некоторых нелинейных про­ цессов, чем это можно было бы предположить на основании рас­ пределения энергии.

Наиболее заметное физическое различие существует между случаями 1-D и 2-D. В случае 1-D завихрения в фазовом прост­ ранстве сливаются до тех пор, пока на длине L не останется один вихрь, а е после стадии насыщения уменьшается до уровня, соот­ ветствующего этой стационарной структуре. Если L велико, то происходит дальнейшее перемешивание, а е выходит на стацио­ нарное соответственно более низкое значение, однако одно завихре­ ние остается и е не возвращается к тепловому уровню. В 2-Б-слу- чае завихрения начинают формироваться в момент насыщения,

исчезают в конце. Соответственно е (t) достигает насыщения в 2-D- случае при более низком уровне, чем в l-D-случае (заметим, что масштаб на графиках разный), и затем быстро спадает почти до начального теплового уровня. Отметим, что насыщение происхо­ дит быстрее во времени при переходе от 1-D- к 2-D- и к З-Б-случаям, по-видимому, из-за (по крайней мере частично) последовательно более высоких уровней начальных шумов. В З-Б-случае эти раз­ личия с l-D-случаем становятся настолько более резко выражен­ ными по сравнению с 2-Б-случаем, что для количественного срав­ нения с соответствующим экспериментом необходимы результаты З-Б-вычислений. Частично сформированные завихрения менее отчетливы, чем в 2-Б-случае, а уровень насыщения е еще ниже. Модельные частицы, показанные на графиках в фазовой плоско­ сти (X, Уд.) для 2-D- и З-Б-случаев, брались из плоского слоя, ограниченного по у для уменьшения влияния мод с конечными ку, которые маскируют появление завихрений.

Конечные графики / (их) на фиг. 7, которые получены с помо­ щью простой гистограммной процедуры без какого-либо сглажи­ вания, оказываются более гладкими для более многомерных расче­ тов из-за улучшенной статистики. Все три кривые / (ѵх) пред­ ставлены на одном и том же интервале скоростей (±40), измерен­ ных в единицах дебаевской длины, поделенной на плазменный период, и показывают, что 2-D- и З-Б-кривые на хвостах, т. е. вне интервала ±30, по существу не изменяются в сравнении с /о (уж), в то время как в l-D-случае хвосты расширяются из-за устойчивой вихревой структуры с большой амплитудой. Эта же самая вихревая структура дает конечную l-D-функцию / (ѵх) с заметным провалом в середине, тогда как в более многомерных численных экспериментах, в которых плазма в это время прост­ ранственно однородна, функция / (ѵх) является по существу пло­ ской на вершине, особенно для З-Б-случая.

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

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

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

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

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