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

ду, что распределение нейтронов по энергиям сильно анизотропно, а скорость центра масс составляет 100 см/мкс. Таким образом, заслуживает внимания вопрос о применимости как МГД-модели, так и разностной схемы.

Для ионной компоненты параметры плазменного фокуса соот­ ветствуют режиму, промежуточному между бесстолкновительным и режимом с определяющей ролью столкновений, при этом ®ciTa ~ 10. (Электронная компонента является заведомо бесстолкновительной.) Так как характерное время ион-ионных столк­ новений, %ц 30 нс, имеет тот же порядок, что и макроскопиче­ ское время развития плазменного фокуса, условия гидродинамиче­ ского приближения выполнены лишь частично. Кроме того, ион­ ный ларморовский радиус at имеет тот же порядок, что и масштаб длины (радиус ~ 0,06 см), и, следовательно, условие Чу — Голдбергера — Лоу, е = aJL <К 1, строго говоря, не выполнено. Од­ нако, пока эти условия не нарушены на порядок или больше, мы можем рассчитывать на то, что гидродинамическая модель будет качественно справедлива.

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

§ 7. Т р е х м е р н ы й код

Существует много физических явлений, для моделирования которых интересны трехмерные МГД-расчеты. В исследованиях по термоядерному синтезу хотелось бы иметь возможность рас­ считать максимальное давление плазмы в таких трехмерных уста­ новках, как стеллараторы или ловушки с магнитной ямой, и про­ следить появление трехмерных неустойчивостей в таких двумер­ ных установках, как токамак или пинчи. Согласно теореме Каулинга [91, 92] в геофизике, только трехмерные движения в ядре могут генерировать магнитное поле Земли. Другие трехмерные задачи возникают в космофизике, физике Солнца и астрофизике.

На сетке 64 X 64 можно выполнять полезные двумерные расче­ ты, пользуясь такими машинами среднего быстродействия, как IBM 7090 или установленная в Калэме ICL KDF 9. При этом один шаг по времени занимает 10—15 с, и потому можно ожидать, что-

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

с помощью какой-либо из новых больших машин типа CDC 7600 или IBM 360/91, имеющих примерно в 64 раза большее быстро­ действие [80], удастся провести трехмерные расчеты на сетке 64 X 64 X 64. Однако гораздо более жесткие требования предъ­ являются к объему оперативной памяти, поскольку нужно пом­ нить восемь функций в 256 тысячах узловых точек.

В настоящем параграфе рассмотрен простой трехмерный яв­ ный код, разработанный Борисом и Робертсом, и изложены мотивы некоторых из принятых решений.

1. Выбор языка и машины

На KDF 9 стоимость составления одной команды на автокоде примерно равна стоимости выполнения ІО7 машинных команд. Написать программу на Фортране или Алголе значительно де­ шевле, чем на автокоде, но последний работает в 2—4 раза быстрее, поэтому естественно программировать на автокоде лишь те команды, которые выполняются более чем ІО7 раз, а осталь­ ные программировать на языке более высокого уровня. Возможно, что для какой-либо из новых быстродействующих машин лучшим критерием будет ІО8. В случае трехмерных вычислений на сетке 64 X 64 X 64 ІО8 операций соответствуют всего 400 шагам по вре­ мени, так что внутренний цикл программы, почти наверное, сле­ дует писать на автокоде, но двумерные расчеты и тем более одно­ мерные, по-видимому, следует программировать на более слож­ ных языках.

Интересно, что уменьшение количества узлов сетки в а раз в каждом измерении позволяет использовать гораздо менее быстро­ действующие ЭВМ и программирующую технику, так как при этом расчет ускоряется в а 4 раз. Например, расчет на KDF 9 с ис­ пользованием Фортрана и сетки 16 X 16 X 16 отнимет столько же времени, сколько расчет на IBM 360/91 с применением авто­ кода и сетки 64 X 64 X 64, причем для первого из них достаточно 32 тысячи ячеек, если в каждую кодовую группу объемом 48 бит записывать два числа. Это дает возможность отлаживать про­ грамму на машинах среднего быстродействия и проводить на боль­

из уравнений (1), (14), (19) и (20) с тензором потока импульса, определяемым (5) и (6). Если коэффициенты переноса постоянны, то время работы программы можно значительно уменьшить, а про­ цесс программирования на автокоде значительно упростить, так как в этом случае диффузионный член принимает вид Д/, что

позволяет применить непосредственно схему Дюфора — Франкела. По аналогичным причинам расчет следует проводить на равномерной прямоугольной сетке по явной схеме.

3. Граничные условия

Граница может быть либо периодической, либо жесткой прово­ дящей стенкой (vj_ = Е II = 0). В последнем случае используется граничное условие свободного проскальзывания (отсутствие натя­ жения) или нулевого проскальзывания (ум = 0). Стенки могут быть либо термоизолирующими (нулевой поток тепла), либо иметь постоянную температуру. Все эти условия легко запрограм­ мировать. Интересно закрепить поток у одной из границ и при­ вести ее в движение по заданному закону, моделируя тем самым закручивание протуберанцев в солнечной хромосфере, связанное с движениями в фотосфере. В уравнение (19) легко добавить гра­ витационный член.

4. Разностная схема «с перешагиванием»

и Т определены в одних и тех же узлах трехмерной кубической решетки типа решетки NaCl (фиг. 4), причем атомам Na соответ­ ствуют четные слои по времени, 2nAt, а атомам С1 — нечетные

фициента avd (р), точно центрированы. Фигурирующие справа первые производные divd, deld и rotd содержат только функции, определенные в наружных узлах FC на промежуточных слоях, а оператор Лапласа можно аппроксимировать по методу Дюфо­ ра — Франкела, связывая тем самым вместе все восемь простран­ ственных ячеек. Коэффициент avd (р) можно вычислять либо как

по времени.

Трудность в вычислении члена Г; заключается в том, что выра­ жение div V не центрировано должны образом в узлах FC. Можно было бы использовать формулу, построенную на 1 2 ближайших

к FC узлах того же слоя (т. е. на серединах ребер куба с центром

вFC), но это заняло бы много машинного времени. Борис пред­ ложил разбить Гг на члены двух сортов, например d^vjdx2, и д2ѵ!дх ду. Первый член можно вычислить по схеме Дюфора — Франкела, а второй — с использованием в качестве узлов четырех угловых точек центрального сечения основного куба (х , г/)-плос- костью, причем два узла следует взять на старом временном слое, а два других — на новом. Такая процедура вносит небольшие ошибки, которые можно устранить, изменяя направление пере­ счета.

Разностная схема «с перешагиванием» описана здесь потому,

что она, видимо, является простейшей из всего, что можно при­ думать. С другой стороны, лучше использовать в (19) дивергенцию тензора, чем члены вида (ѵ-Ѵ) ѵ или j X В; есть также некоторые преимущества в употреблении поля В вместо векторного потен­ циала А. Описанная схема была запрограммирована на Фортране и в настоящее время находится в стадии отладки, но некоторые вопросы остаются нерешенными, например введение искусствен­ ной вязкости для рассмотрения ударных волн и подавление числен­ ных колебаний плотности, которые могут возникнуть из-за отсут­ ствия в (1 ) диффузионного члена.

§ 8. З а к л ю ч и т ел ьн ы е за м е ч а н и я

Мы обсудили численные МГД-расчеты в случаях одного, двух и трех измерений. В случае одномерных цилиндрических задач построены реалистические коды для полностью ионизованной плазмы и достигнуто хорошее согласие с экспериментом — сначала для основной динамики разряда, а затем и для диффузии плазмы поперек магнитного поля. Для того чтобы объяснить ранние ста­ дии разряда Ѳ-пинча, пришлось ввести в код повышенное удель­ ное сопротивление, однако тщательное сопоставление расчета с экспериментом показало, что на более поздней стадии диффузия становится классической [1 0 ].

Двумерные расчеты проводятся уже в течение нескольких лет, но только недавно получены результаты, сопоставимые коли­ чественно с экспериментом. Наиболее продвинуты расчеты пн плазменному фокусу, поскольку в них магнитное поле описывается скаляром B q и имеет простой закон изменения в вакууме, B q ~ 1 /г, и уже первое сопоставление с экспериментом было обнадеживаю­ щим. Главная трудность обобщения на случай произвольной

анизотропии теплопроводности.

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

9.Bodin H. A . B., McCartan J ., Newton A . A . , Wolf G. ff., «Diffusion and stability of high-ß plasma in an 8-metre theta pinch». Third IAEA Conf. Plasma Phys. and Controlled Nucl. Fusion Res. (Новосибирск, август

1968), Paper, CN-24/K-1.

10.Beach A . D ., Bodin H. A . B., Bunting C. A ., Dancy D . J ., Hey wood G. С . I I ., Kenward M . R ., McCartan J., Newton A . A . , Pasco I . K., Peacock R.,

12.Paul J . W. M ., Parkinson M . J ., Sheffield J., Holmer L. S., Nature, 208, 133 (1965).

13.Paul J . W . M ., Goldenbaum G. C ., Iiyoshi A ., Holmes L. S ., Hardcastle R. A .,

Nature, 216, 363 (1967).

14.Roberts К. V., Journ. Nucl. Energy, Pt. C., 5, 365 (1963).

15.Duchs D., «Three-fluid model for a partially ionized plasma in 0-pinch discharges». Sixth Intern. Conf. Ionization Phenomena in Gases (Paris,