(второй индекс у диагональных элементов матриц опускаем). Пара рекурсивных соотношений (76) и (74) используется по оче реди. Правое граничное условие определяет Xj, а левое гранич ное условие определяет Этот метод оказывается достаточно быстрым, так как он содержит небольшое количество умножений на каждом шаге. В общем же виде рассматриваемая задача нели нейна, S = S (ип+1), и для получения самосогласованного реше ния приходится использовать итерации.
5. Неявный метод Хейна
а. Одномерный случай
Система уравнений одномерной магнитной гидродинамики может быть решена неявным методом путем сведения к трехдиаго нальной форме, изученной выше [1]. Представим эту систему в модифицированном, несколько более точном варианте, который, возможно, и не совпадает с обычно применяемым на практике, а также предположим, что используется лагранжева сетка, так что переносом можно пренебречь. Достаточно рассмотреть следую щие простые уравнения:
|
dv |
т, |
dB |
|
(77) |
|
РИ Г = - В -эХ' |
|
|
|
dp |
P |
дѵ |
’ |
(78) |
|
dt |
dx |
|
|
|
dB |
D dv , |
|
dzВ |
(79) |
|
dt ~~ |
B |
|
dx2 " |
|
|
|
Так как эти уравнения нелинейны, придется применить итера ционный процесс; поэтому мы предположим, что и”, ип+1>р изве стны (р — порядковый номер итерации) и что un+1>p+1 должно быть определено. Сущность метода заключается в формальном решении уравнения (71) и подстановке результата в уравнение (77), после чего получается уравнение параболического типа, которое решается методом прогонки при е == 1/2, так что пространственные производные центрированы по времени. Обозначим чертой сред нее между точками (п) и (п + 1 , р), а тильдой — среднее между
(п) и (п + 1, р + 1). Тогда
В At / дѵ \ |
г)Аt / д2В \ |
2 ~ \~ дХ ) |
) ’ |
§ 3. Разностные методы |
361 |
или (группируем первый и третий члены в правой части равен ства)
|
в = ё |
В At |
I |
дѵ |
\ |
|
(80) |
|
2 |
\ |
дх |
) |
' |
|
|
|
Дифференцируя (80) и подставляя результат в (77), находим
|
искомое параболическое |
уравнение |
|
|
|
|
— dv |
и 9В |
ВAt д |
( п 9ѵ |
\ |
(81) |
|
р ч г = ~ в дх |
2 дх |
V дх |
/ ' |
|
|
Исследование устойчивости показывает, что множитель перехода по модулю равен единице и что диффузионный член в уравне нии (81) полностью гасит неустойчивость, возникающую из-за явного представления градиента магнитного давления в правой части уравнения.
На практике удобно вычислять скорость в точках простран ственной сетки с целыми индексами, а все остальные переменные в полуцелых точках, чем достигается надлежащее центрирование первых производных в пространстве. Уравнение магнитного
поля (79) решается с использованием известной величины дѵ/дх во втором члене, а уравнение непрерывности (77) — по формуле
|
At |
дѵ |
At |
\ |
T |
(82) |
|
2 |
дх |
2 |
/ |
J |
|
|
Постоянный шаг по времени подбирается так, чтобы все перемен ные изменялись не слишком сильно за один шаг; для этого его подчиняют общему условию
дѵ At
дх 1 Г < 1 ,
которое определено уравнением (82) и другими описывающими адиабатическое сжатие членами. Обычно бывает необходимо не
более двух |
или трех |
итераций. Уравнения для электронного |
и ионного |
давлений |
записываются аналогично уравнению (80) |
и после дифференцирования подставляются в (81), а если урав нения для полей В е и В z связаны анизотропным сопротивлением, то эти уравнения записываются в форме матриц 2 X 2 и решаются совместно.
б. Случай двух и трех измерений
Хейн [26] обобщил свой метод на случай двух пространствен ных измерений, а позднее и на случай трех измерений (не опуб ликовано). Рассмотрим простейшую систему гидродинамических уравнений:
dp
3 62 |
|
|
Гл. 9. МГД-методы |
|
|
|
|
|
|
|
|
|
dvx ____ j9p |
|
|
|
|
|
(84) |
|
|
|
|
dt |
|
dx |
’ |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
dVy |
|
dp |
|
|
|
|
|
(85) |
|
|
|
|
dt |
|
dy |
’ |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
) |
— ypV .у . |
|
|
|
|
(86 ) |
|
|
|
|
|
|
|
|
|
По |
аналогии с (80) |
представим решение уравнения (86) в виде |
|
|
|
р = рп- ^ р Ѵ - ѵ , |
|
|
|
|
(87) |
что |
после подстановки в (84) |
и (85) |
дает |
|
|
|
|
|
|
- ddvt x _ |
d p n |
, |
УрАі |
( |
d 2v x |
|
d*Vy |
|
(88) |
|
d x |
\ |
d x * |
1 |
d x d y |
|
|
P |
|
d p n |
1 |
y p 2 |
â Z vx |
|
|
|
|
|
n |
d t |
|
h t |
,t |
, ö 4 |
|
' |
|
|
|
1 |
|
( |
|
|
|
(89) |
|
d^v |
|
|
|
|
|
|
2 |
, |
|
P |
|
d y |
2 |
d x d y |
|
d y |
|
|
|
|
|
' |
|
|
|
|
Члены с д21дх2 и d2/ % 2 вычисляются методом прогонки, а сме шанные производные — путем итераций. Уравнения (83) и (86 ) решаются с помощью формул, аналогичных уравнению (82). Компоненты скорости определены в смещенных пространственных
ячейках, т. |
е. |
если р и р заданы в узлах (г, ;'), то ѵх задано в точке |
(і + Ѵ2, /), |
а |
ѵу — в (г, ] + Ѵ2). |
6. Эйлеровы и лагранжевы координаты
Одна из основных проблем в МГД-расчетах возникает из-за преобладающего влияния переносного члена ѵ *Ѵ/, который при сутствует во всех уравнениях и определяет перенос от точки к точке при движении среды как целого таких физических пере менных, как плотность, температура, магнитное поле и скорость, оставляющий их значения неизменными. На этот перенос нало жены и другие кинематические эффекты, связанные с движением среды как целого, такие, как адиабатическое сжатие (div ѵ) и враще ние [например, член (В - V) ѵ в уравнении поля], которые, как и процессы диффузии, могут быть физически малы, но которые тем не менее важно вычислять точно. Очевидно, что лагранжево представление, при котором используется сетка, движущаяся вместе со средой, дает возможность вообще исключить движение среды как целого, что в свою очередь позволяет вычислять остаю щиеся члены более точно. В случае одного измерения такой под ход оказался вполне успешным.
В случае двух или трех измерений применение лагранжевой сетки оказывается делом значительно более сложным, так как она сразу же становится неортогональной и в конце концов может
§ 3. Разностные методы |
363 |
чрезмерно исказиться. Для специфических задач Хейн [42] и Хертвек и Шнайдер [43] разработали псевдолагранжево представление, но обычно на практике используется фиксированная эйлерова сетка, которая требует осторожности в аппроксимации перенос ного члена, чтобы численные эффекты не скрыли физическую диффузию.
а. Точная аппроксимация переносного члена
Численная точность различных аппроксимаций переносного члена изучалась Робертсом и Вайсом [24]. Рассмотрим уравнение
ËL |
df n |
(90) |
V-г—= 0 . |
dt |
дх |
|
Приписываемый Лелевьеру метод ([25], гл. 12) в зависимости от того, V < 0 или V > 0, использует соответственно правую или левую разностную аппроксимацию производной dfldx. Эти раз ности не центрированы и потому аппроксимируют комбинации производных
|
д} |
Дх |
Ö2/ |
|
|
|
ІйГ + |
~2~ Их2 ’ |
|
|
так что уравнение (90) заменяется уравнением |
|
<9/ , |
df |
I |
у I Дх |
92/ |
(91) |
dt |
дх |
|
2 |
9х2 ’ |
|
|
где правая часть отвечает диффузионному или затухающему процессу.
Численную диффузию довольно легко исключить подходящим смещением обеих производных, и существует ряд схем, выполняю щих эту операцию. Те из них, которые будут рассмотрены в этой главе, используют трехслойную схему, так что уравнение (90) решается в виде
fT l-h 1 |
г П — 1 |
V j \ t , .71 |
гП . |
/ П О \ |
і і |
— Ji |
= — ~і^г ш +i —/;-i)- |
(УЛ |
Хотя диффузия и устранена, численная дисперсия все еще проис ходит и коротковолновые моды распространяются со скоростью, отличной от точного значения ѵ. Для уменьшения этого эффекта Робертс и Вайс [24] разработали разностные схемы более высокого порядка точности, однако полностью его устранить не удается.
б. Устранение нефизических значений
Одно из преимуществ схемы Лелевьера состоит в том, что функции, которые всегда положительны, такие, как плотность или температура, остаются положительными в течение всего вре мени вычислений. Это происходит в силу предписанного вклада каждого из узлов в линейную интерполяцию между значениями
364 Гл. 9. МГД-методы
/” и fi± 1 - Иначе можно сказать, что, несмотря на численное рас сеяние, коротковолновые моды так сильно затухают, что исчезают полностью, прежде чем распространятся достаточно далеко, и по тому не вызывают какой-либо ряби в решении.
Если удалить затухание, то короткие волны смогут свободнораспространяться, а поскольку их скорость (а возможно, и ее направление) отлична от физической скорости ѵ длинноволновых мод, они могут дойти до «спокойных» участков области расчетов, где их амплитуда достаточно велика, чтобы дать заведомо невер ный результат. В частности, плотность или температура могут стать отрицательными, что в нелинейных расчетах приводит к опас ным последствиям. В методе Лакса — Вендроффа вводится числен ный диффузионный член четвертого порядка, который оказывает диссипирующее воздействие ([25], гл. 12), но и он не исключает возможности появления нефизических эффектов.
Понадобилось некоторое время, прежде чем эта возможность была осознана и были приняты дополнительные меры [38], впро чем, недавно Хейн *) нашел более изящный подход. Рассмотрим только переносный член в трехмерном случае. Введя обозна чения
/ м а к с — |
ПЗЭХ ( / i j f t |
, |
/ і ± 1 , j ± l , h ± l ) , |
|
/ м и н = |
min {fijh |
, |
fi±l, j + 1, |
f t ± l ) , |
|
заменим величину /ц^ 1 |
на |
|
|
|
|
/* =m ax ( / м и ш min ( / м а к с , |
f m 1))- |
(93) |
Простой смысл этого равенства в том, что новое значение /* не должно выходить из области семи известных значений, исполь зуемых для вычисления /*. Такая процедура препятствует воз никновению ложных максимумов и минимумов, но она вызывает некоторое размытие истинных максимумов и минимумов при их движении по сетке. Другие члены, такие, как сжатие, должны остаться неизменными. Если условие (93) запрограммировать наиболее экономично, то, по нашим оценкам, к затратам машин ного времени для І-кодов типа описанных в § 7 добавится около
30%.
в. Псевдолагранжев метод
В таких экспериментах, как восьмиметровый Ѳ-пинч [9, 10], создается длинный тонкий шнур плазмы, помещенный в почти параллельное магнитное поле, причем длина этого шнура в не сколько сот раз превышает его диаметр. По отношению к радиаль ным смещениям плазма почти равновесна, в то время как вдоль
0 К. Hain, частное сообщение.
