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

Для примера будем считать, что полный поток образуется только за счет разности скоростей ионов и электронов, что дает нам право положить

Пі « пе = п.

Разность VI — \ е состоит из ѵ? (~ mg X В) и \р (~ mdE/dt). Поскольку в dE/dt наибольшие затруднения связаны с членом dEidt, мы оставим \ Р ~ т {дЕІді) и отбросим \ g (обычный дрейф ѵ£>ѵР). Таким образом, приближенно

ние для дивергенции, мы получаем уравнение, куда входят п и ф:

-й-<ѵгф) = - ^ [ ^ - ( ѵ Ѵ ) ] - ^ [ | - ( Ѵ Ф ) - ѵ ге] + ... • (119)

В процессе решения мы вычисляем правую часть на предыдущем и текущем временных слоях и используем полученные величины для определения последующих значений. Как будет показано ниже, подробности зависят от выбора разностной схемы. Два чле­ на, содержащие д (Ѵ2Ф)/dt, не объединяются, так как в нашей двухжидкостной модели это уравнение не решается явно.

Для простоты предположим, что п является переменной нуле­ вого порядка и что ф зависит только от одной координаты х} —

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

где

В правой части можно использовать предыдущий и текущий вре­ менные слои, в левой части — предыдущий, текущий и последую­ щий.

Схема «с перешагиванием» дает уравнение (индекс фурьеразложения к опущен)

фі _ ф- 1 = ( _ 2К + 2iKD) (ф° - ф'1),

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

S = - 2К + 2iKD,

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

X2 - Sk — (1 - S) = 0.

Физическому решению отвечает корень

К= 1 ,

акорень, соответствующий нежелательной вычислительной моде, равен

Таким образом, эта мода оказывается безусловно неустойчивой. Разностная аппроксимация по схеме Адамса — Башфорта дает

ф 1 _ ф 0 = і _ 5 Ң ( ф 0 _ ф - 1 ) _ ^ . ( ф - 1 _ ф - 2 ^ ,

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

X3- X 2( l + j S } + X S — j S = 0.

В этом случае физическому решению отвечает корень X = 1, а ура­ внение

определяет два корня, ответственных за вычислительные моды. Если D пренебрежимо мало, то при К = 0,5 для одной из этих мод выполнено предельное условие устойчивости |Я| = 1 ; для дру­ гой моды |^| = 0,25, и, следовательно, она затухает.

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

Необходимо сделать некоторые замечания относительно воз­ можных значений коэффициента D 1см.(122) ]. Тригонометрическая функция изменяется между 1,0 и 0,0 при возрастании кАх от 0 до л, причем нижняя граница А, = 2Ах определяется условием раз­ личимости на пространственной сетке. Множитель1) Ѵп/кп при­ ближенно равен 1 /кб, где б — толщина пограничного слоя. При малых к, т. е. больших длинах волн, здесь возникают определен­ ные трудности. Действительно, если к мало, то даже схема Адамса — Башфорта может стать неустойчивой из-за некоторого длин­ новолнового воздействия поляризационного дрейфа на электро­ статическое поле. Источник неприятностей заключен в разност­ ном аналоге оператора Лапласа V2. Подобные проблемы возникали и в предшествующих исследованиях устойчивости [2 2 ], где было отмечено, что длинноволновые явления, по-видимому, не вызывают нежелательных последствий. В проделанных расчетах мы не обнаружили неустойчивости, которую нельзя было бы предсказать основываясь на предположении, что D < 1.

Если условие К<^1 не выполнено, то существует дополнитель­ ное затухание физической моды, обусловленное неточностью при­ нятой разностной схемы для dcp/dt при t = 0 , так как аппроксима­ ция (ф° — ф_1)/Д£ имеет погрешность ~ d2q>/dfi. Это затухание уменьшается с уменьшением ©А£, но при ©Аt = 0,3 оно еще достаточно велико, если К ~ 1. Более точное разностное выра­ жение имеет вид

~КГ ( Т Ф° 2ф- 1 + Y Ф-2) >

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

Таким образом, при К <С1 можно разработать различные устой­ чивые разностные схемы для поляризационного члена. Однако данная ситуация не вполне удовлетворительна, так как хоте­ лось бы иметь такие схемы для любого К в области от К<^_1 до К^>1. В следующем подпункте описана модификация рассматри­ ваемой модели для случая К > 1.

Выше мы показали, что затруднения, связанные с вычисли­ тельной модой в нашей двухжидкостной модели, обусловлены поляризационной скоростью ионов ѵр = mil (еВ2) (dTH/dt), кото­ рая является реакцией плазмы на изменение поля Е. Чем больше величина К = Пі/Пі/е^Б2, тем сильнее плазма отзывается на изменения Е. С увеличением К до больших значений результиру­ ющее разделение зарядов становится очень малым (|н; — гсе|<^дг;). Даже если двухжидкостная модель численно устойчива, все рав­ но при /Г$>>1 будет происходить возрастание ошибок округления, связанное с вычислением разности приблизительно равных вели­ чин.

Для того чтобы устранить как ошибки округления, так и чис­ ленную неустойчивость, свойственные двухжидкостной модели, перейдем от двухжидкостного описания к одножидкостному. Вместо раздельного рассмотрения электронной и ионной жидко­ стей введем нейтральную жидкость (п = nt та пе) и заряженную жидкость (пс = Пі — Пе) и будем учитывать только ЕІВ -дрейф; остальные виды дрейфа ионов, vg, ѵр, служат только для разделе­ ния зарядов. (Эти приближения означают, что справедливо усло­ вие Vz^Vg, которое обычно хорошо выполняется при К > 1.)

Можно найти изменение условия устойчивости вычислений. Уравнение (118) (дополненное членами, которые были опущены в предшествующем рассмотрении) в случае одножидкостной модели решается явно, поэтому можно объединить слагаемые с Ѵ2ср, пере­ писав уравнение (119) в форме

(123)

Уравнение (122) переходит в

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

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