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

пенсирующие количество жидкости, протекающее через ячейку j [см. (ИЗ)]. Такая схема имеет довольно неприятную особен­ ность: соседние ячейки оказываются разъединенными (четные ячейки компенсируют поток жидкости в других четных ячейках, а нечетные ячейки компенсируют поток жидкости в других нечет­ ных ячейках). Для одной специальной модели было показано, что отмеченная особенность является причиной нелинейной вычисли­ тельной неустойчивости на коротких волнах [171. В работе 1181 показано, что эта особая неустойчивость может быть подавлена применением специальных пространственно-разностных схем. Воз­ можно применение и многих других разностных аппроксимаций по пространственным переменным. Известно, что расчеты для не­ линейных уравнений жидкости часто наталкиваются на серьезные проблемы устойчивости вычислений (проявляющиеся в виде большой величины ложной энергии в коротких волнах), которые возникают из-за слагаемых с разностными производными по координатам. Во многие модели для подавления коротких волн вносится искусственное затухание. Один из таких методов основы­ вается на специальной разностной схеме, определяемой уравне­ нием (17). Для этой схемы множитель роста, получаемый при анализе устойчивости, имеет вид

I К\ = 1 — О (А;Ах)4.

Знак минус означает, что схема будет иметь тенденцию к зату­ ханию волн, а зависимость вида (&Ах)4 означает, что короткие вол­ ны будут подавляться наиболее сильно.

Результаты коротких машинных вычислений, приемлемые для простых «частных случаев», нельзя переносить на другие пробле­ мы. Другими словами, полный расчет должен быть по возможно­ сти кратким, чтобы за это время погрешности с малой длиной волны не выросли до опасных уровней. Отдельные кратковремен­ ные вычисления для обсуждаемой модели были успешно проведе­ ны с незначительным влиянием или при полном отсутствии про­ блемы коротких волн [13, 141. Другие расчеты для этой же моде­ ли, проведенные в более сложных случаях, требовавших большей продолжительности машинного счета, поставили проблему устой­ чивости вычислений как раз такого типа. Полностью успешные вычисления требовали определенной формы контроля за корот­ кими волнами.

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

в. Процедуры дифференцирования но времени. Использование комбинированных схем

Рассмотрим вычислительную неустойчивость конечно-разност­ ных уравнений, являющихся аналогами уравнения

(114)

которое есть общая форма уравнений непрерывности. Погреш­ ность аппроксимации находится из сравнения уравнения (114) с разложением в ряд Тейлора конечно-разностной формы этого уравнения. Полная погрешность, возникающая за счет накопле­ ния погрешности аппроксимации для каждого шага по времени, должна оставаться малой. Практически все расчеты нелинейных жидкостей, когда в схемах с дифференцированием по времени для получения новых значений используется информация на двух временных уровнях, ограничиваются памятью машины. Однако до сих пор существует большое число схем такого класса, и они сильно отличаются по точности и устойчивости. «Наилучший» выбор схемы обычно зависит от фактического типа ожидаемого решения.

Изучим устойчивость (накопление погрешности вычислений) отдельных схем,, применяемых для решения уравнения колебаний,

В работах [18, 19] приводится сравнительное описание раз­ личных схем, основанное на применении их для решения урав­ нения (115). В этом подпункте мы построим и проанализируем комбинированные схемы, которые сохраняют положительные свой­ ства отдельных схем. Аналитическое решение уравнения (115) есть, конечно,

и (t) = ц0 ехр (іозі).

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

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

294 Г л . 7 . М о д е л и п л а з м ы без с т о л к н о в е н и й

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

Схема «с перешагиванием» (LF), записанная для общего урав­ нения (114), для временных уровней, отмеченных верхними индек­ сами, имеет вид

2М — ’

причем погрешность аппроксимации пропорциональна d3u/dt3. При получении решения конечно-разностного уравнения мы сле­ дуем обычной методике [1 ], когда разностное решение и (nAt) представляется в виде

ип = Хпи°,

где величины X носят название множителей роста. Таким образом, схема «с перешагиванием» записывается для уравнения (115) в виде

и1 — и- 1 = 2 ibu°,

где b = (£>At; характеристическое уравнение дает

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

Схема Адамса — Башфорта (AB) в применении к уравнению (114) имеет вид

и1 — ио

Аt

ипогрешность аппроксимации также пропорциональна d3u/dt3. Записывая эту схему для уравнения (115), мы получаем

и1— и° = у іЪи°— ibu-1,

что приводит к дисперсионному уравнению

Xz- X ( i + ^ i b ) + ± - i b = 0.

Когда Ь<^ 1, оно дает

|М = 1 + 0 (Ь * )+ ...,

і ч = 4 < і ,

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

Комбинированные схемы, в которых используются различные комбинации одношаговых схем для полного шага по времени, испытывались по отношению к уравнению (115). Требовалось получить частичную компенсацию погрешностей аппроксимации и подавить вычислительные моды. Суммарный эффект, однако, является более тонким, чем простое сложение погрешностей противоположного «знака». Например, в применении к уравне­ нию (115) схема «с перешагиванием» (нет роста) схема Адамса — Башфорта (медленный рост) дают вместе медленное затухание.

Подробности анализа устойчивости комбинированных схем можно найти в работе Байерса [20]. Было проанализировано много различных комбинированных схем и открыто несколько превос­ ходных схем. В качестве примера получаемых улучшений сравним схему Адамса — Башфорта с комбинированной схемой LF AB AB, т. е. схемой, использующей «перешагивание» для 1-го, 4-го и т. д. шагов по времени и схему Адамса — Башфорта для 2-го, 3-го, 5-го, 6-го и т. д. шагов по времени. Мы сравнивали квадрат ампли­ туды, R2, в момент времени оit = 50,0, т. е. после примерно восьми полных волновых периодов. Если Ъ = 0,2, то для схемы Адамса — Башфорта в этот момент R 2 — 1,24, а для комбинированной схемы R 2 = 0,99. Таким образом, видно, что схема Адамса — Башфорта увеличивает энергию волны на 24%, тогда как комбинированная схема уменьшает энергию волны всего на 1 %.

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

г. Разностно-временная неустойчивость из-за вычислительных мод

Хорошо известно, что схема «с перешагиванием» неустойчива для уравнения *

В этом случае характеристическое уравнение имеет вид

№ + 2 ЬХ — 1 = 0 , где b = соАі,

откуда

и = — ъ ± (і + б2)1/2.

Заметим, что |Я_| > 1, поэтому вычислительная мода нарастает и вскоре по своей величине превосходит решение. Сопоставим это со схемой Адамса — Башфорта, для которой характеристиче­ ское уравнение записывается в виде

Х2- Я ( 1 - | б ) - 1 б = 0,

или

* * - Н 1- т г,) ± - Н 1- |,+ т г'! Р '

и заметим, что |А, | <£ 1 для малых Ъ. Это означает, что использова­ ние схемы Адамса — Башфорта приводит к сильному затуханию вычислительной моды.

Для наших уравнений существует опасность появления неус­ тойчивости из-за нарастающих вычислительных мод при исполь­ зовании схемы «с перешагиванием» или комбинированной схемы, в которой «перешагивание» преобладает. Пользуясь методикой, описанной Рихтмайером и Мортоном 11], легко сформулировать проблему полного анализа устойчивости задачи в целом, однако решение ее затруднительно. Наилучший результат для большин­ ства сложных систем достигается с помощью разбиения всей системы на отдельные части, после чего определяются условия устойчивости каждой подсистемы. Хотелось бы надеяться, что суперпозиция таких отдельных условий устойчивости будет доста­ точным условием устойчивости всей системы. В дальнейшем будем использовать именно этот метод. (Стоит, однако, указать на то, что Касахара [2 1 ] привел пример, когда подобная процедура разбиения на подсистемы оказалась непригодной.)

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