14.13. НЕЯВНЫЙ МЕТОД ДЛЯ БОЛЬШИХ ВРЕМЕННЫХ ШАГОВ

Характерные временные масштабы коллективных плазменных явлений охватывают много порядков по величине. На наиболее быстрых временных масштабах, где решающую роль играют кинетические эффекты, для изучения нелинейной эволюции плазмы с большим успехом применялись методы частичного моделирования. Однако как для приложений, так и для фундаментальных исследований все больший интерес представляет обобщение техники моделирования кинетических явлений на более медленные масштабы времени. Один подход к моделированию медленного поведения состоит в изменении исходных уравнений для исключения не интересующих нас высокочастотных мод. Например, в частичном моделировании применяется электростатическое приближение или приближение Дарвина. Другими подходами являются уже упоминавшиеся в § 9.7 организация подциклов, усредненные орбиты и неявное интегрирование по времени.

Самые большие частоты в плазме — это электронная ленг-мюровская и электронная циклотронная сосе. Наименьшее время — это время пробега электронов или света на характерное расстояние. Но большие масштабы времени могут определяться инерцией ионов, электромагнитными эффектами и большими пространственными размерами. Соотношение электронных и ионных плазменных и циклотронных частот, а также магнитогидродинамических и электронных времен определяется малым параметром где -кратность ионизации. В тех случаях, когда преобладают силы, обусловленные магнитным полем и текущими по плазме токами, частоты уменьшаются по сравнению с по меньшей мере в раз, где скорость света и тепловая скорость электронов. В этом параграфе применительно к плазменным кодам описан способ нахождения неявных временных конечных разностей, при. котором продвижение по времени осуществляется непосредственно при помощи линеаризованных уравнений движения частиц без использования уравнений для жидкости (момента скорости) [Langdon, Barnes, 1985].

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

Рис. 14.21. Координаты, используемые для перемещения частицы в цилиндрических координатах

Анализировались и применялись различные варианты конечных разностей для частиц [Cohen, Langdon, Friedman, 1982]. Здесь мы обсудим только -алгоритм (см. 9.8), называемый также схемой 1 [Barnes et al., 1983], который записывается в виде

где

Аналогично § 9.8 для проверки точности этой схемы можно получить и разрешить дисперсионное соотношение для гармонических колебаний:

При получаем [Cohen, Langdon, Friedman, 1982]

И постороннюю моду При моды сильно затухают:

Добавляя в уравнении совершая поворот и снова добавляя можно найти В результате получаем

где описывающий вращение на угол — оператор имеет следующий вид:

причем -единичный тензор. При малых

При больших

Прямой метод для электростатических полей: решение неявных уравнений. Вторая основная проблема при использовании неявных кодов заключается в более сложной организации временного цикла. При использовании явных разностей временной цикл разделен между продвигающимися частицами и полями — эти вычисления чередуются и проводятся независимо. Однако неявный код должен решать зацепленную систему уравнений (8.22), (8.23) и (3.5) совместно с (14.91), (14.92) или (14.93), (14.94). Из-за того, что используется электрическое поле во всех неявных схемах будущие координаты зависят от ускорений Но это поле еще не известно, так как оно зависит от плотности частиц с координатами Задача состоит в решении этой большой нелинейной системы связанных уравнений для частиц и полей.

Исторически в первом методе, применявшемся для решения данной задачи, поля на новом временном уровне предсказывались при помощи решения связанной системы уравнений для полей и жидкости, причем кинетический тензор напряжений вычислялся из известных в более ранние моменты времени скоростей частиц. После того как поля становятся известными, частицы передвигаются на новый временной уровень, при необходимости снова вычисляется тензор напряжений и процесс итерируется. Такой моментный подход детально описан в работах [Denavit, 1981 и Mason, 1981].

Описываемая здесь модификация этого метода заключается в непосредственном предсказании будущих полей при помощи линеаризованных уравнений движения для полей и частиц. Одна из форм этого метода, а также его улучшение и некоторые примеры, подтверждающие его действенность, приведены в работе [Friedman, Langdon, Cohen, 1981]. Другая форма метода описана в работе [Barnes et al., 1983]. Ленгдон, Коэн и Фридман [Langdon, Cohen, Friedman, 1983] существенно обобщили алгоритм и рассмотрели много важных деталей, таких как пространственная дискретизация и фильтрация, а также итеративное решение неявных уравнений.

Сущность прямого метода состоит в том, что он имеет дело непосредственно с уравнениями движения частиц и с уравнениями, связывающими частицы и поля. Последние линеаризуются в окрестности своих ожидаемых значений (экстраполяция) на новом, временном уровне. Будущие значения разделяются на две части: 1) приращения зависящие от (неизвестных) полей на временном уровне, и 2) экстраполяции объединяющие все прочие вклады в уравнения движения.

Плотность заряда собирается по координатам затем вычисляются коэффициенты в выражении для разности между плотностями, полученными после интегрирования с полем и с корректированным полем В результате получается источник в законе Гаусса

Это — линейное эллиптическое уравнение с переменными коэффициентами для определения или

Точность выражения для приращений определяется компромиссом между сложностью и быстрой сходимостью [Langdon, Cohen, Friedman, 1983; Barnes et al., 1983]. При необходимости для каждого сорта частиц можно вычислять непосредственно как производные «прямая разность» [Langdon, Cohen, Friedman, 1983] или как упрощенное разностное представление (8.22) [Langdon, Cohen, Friedman, 1983; Barnes et al., 1983].

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

где контролирующий неявность параметр, для алгоритма равен 1/2; -координата, получающаяся из уравнений движения при пренебрежении ускорением выражена через координаты, скорости и ускорения в моменты времени, предшествующие и равные Исключая из (14.91) и получаем

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

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

Если разложить в (8.22) по координате, то

при этом Ускорение частицы выражается через при помощи соотношения (8.23), взятого в точке

Записывая (14.103) в виде

получаем, что соотношение (14.103) является конечно-разностным представлением выражения

Следовательно, мы получаем эллиптическое уравнение в частных производных

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

При экстраполяции плотности заряда и разумном конечно-разностном представлении линеаризованного неявного вклада одномерные уравнения поля принимают вид

Здесь используются два представления

или

где

— сумма по всем сортам частиц [Langdon, Cohen, Friedman 1983]. Как в (14.108), так и в определяется выражением (14.102).

Из при полуцелых координатах образуется

При помощи (14.104) ускорение частицы в точке вычисляется через Этот алгоритм — самая короткая схема, которую мы когда-либо видели, и он был опробован Ленгдоном и Барнсом [Langdon, Barnes, 1985].

Общий электростатический случай. Возвращаясь к многомерному случаю, включающему, возможно, создаваемые внешними токами магнитные поля, продемонстрируем, начиная с (14.91) и (14.92), необходимые этапы вычислений. Начнем с формулировки метода в наиболее общем виде. Координаты из которых при помощи (8.22) вычисляется экстраполированная плотность получаются из уравнений движения, в которых определяется полем что является предсказанием для Плотность заряда при этом не совсем точно соответствует полю т. е. Мы хотим вычислить подправленное поле и переместить частицы в положения так, чтобы плотность заряда удовлетворяла соотношению

Перепишем для этого (14.112) в виде

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

При помощи (14.107) и уравнений движения выразим как линейный функционал В общем случае поправки получаются при помощи линеаризации каждого уравнения движения координат х [Langdon, Cohen, Friedman, 1983; Barnes et al., 1983]:

и

без магнитного поля или

в случае замагниченной плазмы, где, как следует из (14.95),

С учетом (14.116) неявный член в (14.113) принимает вид

Суммирование проводится по сортам частиц, а не по каждой частице.

Если неявно рассматриваются только электроны, то только они и дают вклад в (14.112). В этом случае для вычисления суммы нужна электронная плотность [помимо полной плотности которая используется в правой части (14.113)]. В общем случае достаточно накапливать отдельно для каждых сортов частиц с различным отношением Это требует большей памяти, но не больших вычислений, чем при использовании явного кода.

Неявная восприимчивость

из-за вращения является тензором.

Теперь у нас есть все записи уравнения Подставляя выражения в уравнение (14.113), получаем электростатическое неявное уравнение поля

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

Этот формализм позволяет успешно использовать пространственное сглаживание [Langdon, Cohen, Friedman, 1983; Barnes et al., 1983]. Если задаваемое оператором сглаживание применяется к на сетке, то для учета влияния сглаживания на поле следует включить Это оказалось важным в некоторых приложениях. Несогласованное сглаживание существенно влияет на линейную устойчивость.

При перемещении частиц к их положениям поля в (8.23) вычисляются в точках Возникающая при этом погрешность может приводить к ограничениям на связанным с градиентами поля и плотности [Langdon, Barnes, 1985 и § 9.7 настоящей книги]. Тем не менее полезные результаты получаются при значениях достигающих .

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

Неявные коды допускают использование сложных граничных условий для частиц. Поглощение или испускание частиц поверхностью зависит от тем самым граничные условия для частиц оказываются связанными с неявными уравнениями поля. Оказывается, что для электромагнитных полей можно приспособить методы, используемые в явных кодах; например Адам и Гурден использовали неявный подход для граничного условия с излучением волны [Lindman, 1975].

Гораздо больше деталей читатель может найти в работе [Langdon, Barnes, 1985] и в приведенных там ссылках.