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

370 Гл. 9. МГД-методы

Точность и устойчивость этой схемы можно исследовать на простом уравнении переноса

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

Поэтому схема устойчива, если At удовлетворяет условию Куран­ та — Фридрихса — Леви, и дает небольшое, но полезное затуха­ ние четвертого порядка для всех мод, для которых аЯ + ßp Ф 0.

Однако из фиг. 4 видно, что двумерная схема имеет две несвя­ занные группы ячеек (обозначенных на фиг. 4 цифрами 1 и 3), а трехмерная схема — четыре такие несвязанные группы. Чтобы избавиться от нежелательных численных мод, нужно какимнибудь способом связать эти ячейки.

д . С хем а « с п е р е ш а г и в а н и е м »

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

по модулю равен единице. Здесь в случае N измерений 2N ячеек оказываются несвязанными 155].

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

линий, так как в этом случае и явный, и эйлеров методы, по-ви­ димому, непригодны.

Описанный выше метод не связывает ячейки 1 и 3 на фиг. 4, но он рассмотрен здесь потому, что, периодически применяя неко­ торую дополнительную фильтрующую процедуру, можно уничто­ жить нежелательные моды. При использовании для этой цели физической диффузии можно вычислять производные по четырем угловым точкам NE — N W на слоях t2m (фиг. 5). При этом мы получим сложную консервативную схему, в которой динамиче­ ские члены связывают ячейки с центрами в узлах С, FN — FW (фиг. 5). Коэффициенты диффузии нужно вычислять в четырех точках M NE (mid-northeast) — MNW, например попарно усред­ няя значения в узлах N, Е, S, W. Однако для плотности все же необходимо ввести некоторые дополнительные построения.

б. Схема Дюфора — Франкела

Если для аппроксимации динамических членов используется схема «с перешагиванием», то для описания диффузии естественно

Эта схема использует сетку типа NaCl, аналогичную показанной на фиг. 4, однако функ­ ция магнитного потока А и скорость ѵ расположены в шахматном порядке относительно остальных переменных. Каждый из основных членов связывает различные комбинации ячеек сетки, а якобиан вычисляется в повернутой системе координат (а, ß).

чин, средних между слоями і2ГПи /2,71+2. Если коэффициенты пере­ носа зависят от пространственных координат и времени, то обеспе­ чить выполнение условия консервативности схемы довольно труд­ но, и возможно, что метод «с перешагиванием» имеет преимущест­ во перед методом Лакса — Вендроффа только в случае постоянных коэффициентов. Схема Дюфора — Франкела устойчива при всех At и обладает достаточной точностью, если мало отношение y\At/(Ax)2.

374 Гл. 9. МГД-методы

9. Специальные разностные методы

а. Область внешнего вакуума

Во время пинчевого разряда плазма отрывается от стенок возрастающим магнитным полем и обычно окружена областью низкой плотности, физические свойства которой трудно рассчи­ тать с помощью МГД-уравнений. Большинство лагранжевых раз­ ностных схем предсказываает, что вначале прилегающий к стенке слой частиц плазмы будет двигаться внутрь, образуя у стенки область вакуума, но этому результату нельзя всецело доверять, так как у границы плазмы точное аналитическое решение диф­ ференциальных МГД-уравнений обладает неприятными сингуляр­ ностями, которые уничтояшются конечно-разностной сеткой. Например, Розенблют и Кауфман [77] изучили случай Те = T t и установили, что температура становится сингулярной. Физиче­ ски это происходит из-за того, что электропроводность а не зави­ сит от плотности, и потому при р — 0 джоулево тепло оЕ2 распре­ деляется между все уменьшающимся числом частиц. Более того, а увеличивается с увеличением Т, еще больше увеличивая таким образом нагрев, в то время как поперечная теплопроводность падает и отвод тепла затрудняется. Эта сингулярность становится еще острее, если Те Ф Т ;. Однако в действительности МГД-опи- сание становится непригодным вблизи границы плазмы, и этот парадокс, вероятно, может быть разрешен введением ряда повы­ шенных коэффициентов переноса, которые сгладят изменение плотности и устранят другие сингулярности.

В реальном эксперименте плотность вне основного разряда, вероятно, невелика, но едва ли она обращается в нуль, посколь­ ку некоторые заряженные частицы должны уходить во внешнюю область вследствие различных неклассических процессов; кроме того, из стенки выделяются нейтральные частицы, часть которых ионизуется. Поскольку альфвеновская скорость велика, эта область будет находиться в почти равновесном состоянии, Ѵр ^ « j X В, что в свою очередь означает та 0, поскольку р прене­ брежимо мало. Для слабых продольных токов эта область может быть хорошим проводником, но сильные токи должны привести к аномальной диффузии поля, направленной к восстановлению вакуумной конфигурации поля. В отсутствие сколько-нибудь определенных экспериментальных данных наиболее удовлетвори­

где Ѵ24е = 0. Таким образом, для установок с плазменным фоку­ сом, которые имеют только Б ѳ-составляющую поля, рассмотрение вакуумной области не представляет трудностей, но для аксиально­ симметричной системы общего вида на каждом шаге по времени приходится решать уравнение Лапласа в сложной изменяющейся области между границей плазмы и внешними проводниками.

Внекоторых установках эта область может быть неограниченной.

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

Ѵ2А = 0, где В = rot А.

б. Решение методом последовательной верхней релаксации

Хейн [26] предложил использовать функции Грина, но на практике уравнение Лапласа можно решить достаточно быстро методом итераций, не слишком увеличивая полное время расчета. В методе последовательной верхней релаксации (SOR) [78] ска­ лярное уравнение решается на прямоугольной двумерной сетке с помощью алгоритма:

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

а спектральный радиус (декремент наиболее устойчивой моды) — равенством

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

В практических кодах векторный потенциал А п+1 вычисляется по методу Лакса — Вендроффа там, где р ^ р Мин> а затем методом

области R используется как краевое условие, а старое значение А п