Метод предиктор-корректор для решения задач магнитной гидродинамики*

Вычислительные технологии

Том 11, часть 1, Специальный выпуск, 2006

МЕТОД ПРЕДИКТОР-КОРРЕКТОР ДЛЯ РЕШЕНИЯ ЗАДАЧ МАГНИТНОЙ ГИДРОДИНАМИКИ*

В. М. Ковеня

Институт вычислительных технологий СО РАН, Новосибирск, Россия

e-mail: kovenya@ict.nsc.ru т. в. козлинская Новосибирский государственный университет, Россия e-mail: ktv-sun@gorodok.net

An economical difference scheme of predictor-corrector type for numerical modeling of magnetohydrodynamie equations is proposed. A numerical simulation of scattering of plasma in a non-uniform magnetic field is conducted.

Введение

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

В настоящей работе для решения задач магнитной гидродинамики предлагается разностная схема типа предиктор-корректор, где на этапе предиктора строится схема расщепления по физическим процессам и пространственным направлениям, обладающая большим запасом устойчивости, а на этапе корректора восстанавливается консервативность. Заметим, что идеология расщепления предложена в работах [1, 2]. Построенная схема безусловно устойчива (в линейном приближении), что позволяет эффективно использовать ее для расчета нестационарных задач магнитной гидродинамики и физики плазмы. Проведены расчеты разлета высокотемпературной плазмы в магнитном поле.

* Работа выполнена при частичной финансовой поддержке Российского фонда фундаментальных исследований (грант № 05-01-00146) и Интеграционного проекта СО РАН (№ 116).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.

1. Исходные уравнения

Движение плазмы в неоднородном магнитном поле будем рассматривать в одножидкост-ном магпитогидродинамическом приближении в предположении, что эффекты вязкости и теплопроводности малы и ими можно пренебречь. Тогда система уравнений магнитной гидродинамики включает уравнения неразрывности, движения, полной энергии и уравнения магнитного поля [3]. В цилиндрической системе координат (г, г) в осесимметричном приближении уравнения магнитной гидродинамики в безразмерном виде могут быть представлены в векторной форме:

ди

— + 1¥г + УУ, = О,

дЬ

/

где \УГ = -

д_

дг

\

д др дВ

— (ГПУ2г) + к0г^ + Ш\Вг —-дг дг ^ дг

д дВ

— (гпугУх) - М дг \ дг

дд Т- [гугр) + 1р—ГУг дг дг

0

д

—г {ухВг - угВх)

\

(

д_

дг

\

ТГ {пугух) - М\в^ я дг „ дг

д , , др л дВг

— (пу1) + кот^ + м

дг

\

дг

д ду

— {УГВХ - ухВг) дг

0

II = (п, пуг, пуг, р, Вг, Вг)т-,1 = 7— 1; Ма = £>о/\/ 4:7тпоЩгПг — альфвеновское число Маха; к0 = Т0/(и0mj), mj — масса иона. Система уравнений (1) замыкается уравнением состояния в виде

р = пТ.

При преобразованиях уравнений (1) используется условие

ЛуВ = 1А (гВг) + ^

г дг дг

0.

При переходе к безразмерному виду пространственные координаты нормировались на характерную длину время — на Ь/и0, а плотность, температура, скорость и магнитное поле — па их невозмущенные значения п0, Т0, и0, В0. При построении разностной схемы на этапе предиктора использованы уравнения в недивергентной форме. Выбирая в качестве искомых функций вектор / = (п, уг, ух, р, Вг, Вх)т, перепишем (1) в виде

д/

дЬ

+ С / + С / = 0.

Здесь

С г —

( д

УгТГ

дг 0

0

0

0

0

пд

г дг д

УгТГ

дг 0

7р д

г дг 0

г дг

г —-

0

0 д дг 0

0

г дг

0

код_

п дг 0

д

УгТГ

дг 0

0

0 0 0 0

уг д г дг 0

п дг

0

0

д_

дг '

\

0

у

С

(Уг 0 п 0 0 0

0 Уг 0 0 -МА Вг /П 0

0 0 Уг ко/п М\ВГ /п 0

0 0 IV Уг 0 0

0 -Вг Вг 0 Уг 0

0 0 0 0 0 Уг

д_

дг

С

д_

дг'

При записи матричных операторов Сг и Сг использовано условие (3).

Целью настоящей работы является построение экономичных схем, требующих для своей реализации минимального числа арифметических операций на решение задачи. Это может быть достигнуто на основе метода предиктор-корректор, где на этапе предиктора используется метод расщепления, обладающий большим запасом устойчивости, а на этапе корректора повышается точность метода и восстанавливается консервативность схем. Известно, что расщепление операторов приводит к появлению дополнительных членов в разностной схеме, что, как правило, вызывает понпженпе точности расчета [4]. Поэтому при построении схемы будем выбирать такие расщепления, которые приводили бы к минимальному числу этих членов при сохранении свойств безусловной устойчивости схемы и ее экономичности, т. е. скалярной разрешимости разностных уравнений на дробных шагах.

2. Разностная схема для одномерных задач

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

д/ „д/ п ди дШ п

— система одномерных уравнений магнитной гидродинамики в виде (1), (4), векторы искомых функций / и и совпадают с двумерным случаем, а дШ/дг = Шг, Сд/дг = Сг. Решение задачи (5) будем находить в области Q = {0 < £ < Т, 0 < г < X}, в которой введем равномерную сетку с шагами т и Н. Аппроксимируем первые производные д/дг разностным оператором с порядком 0{Кк). Разностная схема предиктор-корректор

/п+1/2 _ /п

--— + СгаЛ/га+1/2 = о,

ип+1 _ ип

---—|- = о

т

аппроксимирует исходные уравнения (5) с порядком 0(тт + Нк), оде т = 2 при а = 0.5 + 0(т), и па этапе предиктора реализуется векторной прогонкой, как следует из вида оператора С. При а > 0.5 она безусловно устойчива. Для получения схем, реали-

С

С = С1 + С2,

где расщепление зададим в следующей форме:

С1 =

Здесь к1 = к0/п, к2 = Мд/п.

0 0 0 0 0 0 /и* 0 п 0 0 0

0 0 0 0 0 0 0 V* 0 0 _к2 В 0

0 0 у* к1 0 , С2 = 0 0 0 0 0 0

0 0 7Р 0 0 0 0 0 0 V* 0 0

0 0 Вг 0 0 0 0 0 0 V* 0

0 0 0 0 0 0 0 0 0 0 0

Разностная схема

/ п+1/4 _ / п

та

+ С 1Л/г+1/4 = 0, и п+1 _ ип

/

п+1/2 _ /п+1/4

та

+ Л^п+1/2 = 0

+ С2 Л/п+1/2 = 0,

т

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

(пп+1/4 _ пп)/(та) = 0, К+1/4 _ <)/(та) = 0,

(^п+1/4 _ <)/(та) + ^п+1/4 + к1Лрп+1/4 + к2ВгпЛВп+1/4 = 0, (рп+1/4 _ рп)/(та) + 7рпЛ^п+1/4 = 0, (Вгп+1/4 _ Вгп)/(та) + £гпЛ^+1/4 = 0, (Вп+1/4 _ Вп)/(та) = 0.

Исключая рп+1/4 и Вп+1/4 го уравнения движения для компоненты Уп+1/4, получим относительно нее разностное уравнение

[I + та<Л _ т2а2

Его решение находится скалярной прогонкой, после чего явно определяются остальные функции па слое п + 1/4. Подобным образом вычисляются значения функций на слое

, 1 /П п+1/2 „п+1/2

п+ 1/^^^шючая IV из уравнения для Вг , приходим к разностному уравнению вида (10) и решаем его скалярной прогонкой. Все остальные уравнения решаются независимо друг от друга скалярными прогонками. Наконец новые значения ип+1 находятся явно из последнего векторного уравнения схемы (8).

Таким образом, решение системы из шести уравнений при переходе с п-го на (п + 1)-й слой находится шестью скалярными прогонками. Введение расщепления (7) приводит в отличие от нефакторизованной схемы (6) к дополнительным членам лишь в уравнении движения для у*, так как матрица С1ЛС2Л содержит ненулевые элементы лишь в третьей строке. Все другие формы расщепления (7) приводят к большему числу дополнительных членов.

Покажем, что разностная схема (8) для линейных уравнений, полученных из (5) "замораживанием" коэффициентов, безусловно устойчива. Для линейных уравнений разностная схема (8) может быть переписана в виде

/п+1 _ / п

т

где индекс "0" в матричных операторах означает, что коэффициенты этих матриц постоянны. Пусть для определенности Л — симметричный оператор. Будем искать численное решение схемы

[(1 + гу¿)\ - 1]2 {[(1 + %ув)Л - 1]4 - %В\е2<1ъ [(1 + гу¿)\ - 1]2 Л2+

+ [(1 + %^)Л - 1]2 [(1 + %ус[)Л - %уд\ (В2 + В2 + с2) ¿2Л +

+В?В^ [

где у = у2, с2 = 7р/п, ^ = (т/К)вт кК. Очевидно Л1 = Л2 = 1/(1 + %уд) и |Л51 < 1 при 5 = 1, 2. Оставшиеся корпи могут быть получены из решения уравнения четвертого порядка (выражение { } в (12)). В общем виде получить его решение затруднительно, однако все предельные случаи (М >> 1, М << 1, МА >> 1, МА << 1) показывают, что |Лв| < 1 (^ = 3,6). Можно ожидать, что разностная схема безусловно устойчива и в остальных случаях. Схема (8) является обобщением схемы [5], предложенной для решения уравнений газовой динамики.

3. Схема для многомерного случая

Идеология расщепления при построении разностной схемы для одномерных уравнений легко может быть обобщена и на многомерный случай при совместном расщеплении уравнений по пространственным направлениям и расщеплении каждого одномерного оператора. Пусть т, Кг, К2 — шаги сетки по времени и пространству соответственно. Введем расщепление каждого из операторов Сг и С2 па сумму операторов таким образом, чтобы соответствующая разностная схема реализовалась скалярными прогонками, была безусловно устойчивой и содержала минимальное число дополнительных членов. Тогда по аналогии со схемой (8) разностная схема для решения двумерных задач может быть записана в виде

£п+1/8 _ £п £п+1/4 _ £п+1/8

--— + С1 г/га+1/8 = о, -±-+ с2 г/га+1/4 = о,

та ' та '

£п+3/8 _ £п+1 /4 £п+1/2 _ £п+3/8

1-1-+ сип+3/8 = 0, 1-1-+ С2 Г^/2 = 0,

та

та

ип+1 — ип

т

+ w;+1/2 + жп+1/2 - ^ г+1/2 = 0,

где С/г'Г = С\т + С2Г, С/г2 = + С2 2; Wh,г, Wh2 — разностные операторы, аппроксимирующие соответствующие дифференциальные операторы с порядком ), % = г, г:

С1

( уг Лг (п/т)Лгт 0 0 0 0 \

0 уг Лг 0 (ко /п)Лг 0 (М 1/п)ВЛг

0 0 уг Лг 0 0 0

0 (7р/г)Лг г 0 0 0 0

0 0 0 0 (уг/г)Лг г 0

0 (В2/г)Лг г 0 0 0 0

C1

0 0 0 0 0 0 \

0 0 0 0 0 0

0 0 0 0 0— (M\/n)BrAr

r 0 0 0 Vr Лг 0 0

0 0 0 0 0 0

0 0 — (Br/г)Л rr 0 0 Vr Лг 1

/Vz Лг пЛ z 0 0 0 0 \

0 Vz Л z 0 (ко/п)Л z0 -(Mi/n)BzAz

0 0 Vz Лz 0 0 0

0 yMz 0 0 0 0

0 0 0 0 Vz Л; 0

V0 Bz Лz 0 0 0 0 /

0 0 0 0 0 0 ^

0 0 0 0 0 0

C 2 = 0 0 0 0 0 0 0 Vz Лz 0 {M\/n)BrAz 00

0 0 0 0 0 0

0 0 — Br ^^-z 0 0 Vz Лz /

Здесь Aj, Aj - разностные операторы, аппроксимирующие первые производные с порядком O(h'j), j = г, z. В качестве Aj, Aj могут выбираться симметричные трехточечные операторы (k = 2), подобно тому, как представлено в разд. 2, пли несимметричные операторы с учетом знака скорости (k = 1), задаваемые по формулам [4, 5]:

Vj > 0 : Aj = А3_, Л; = V; Vj < 0 : Aj = Aj+, Л,- = Л; [j = г, z).

Разностная схема (13) аппроксимирует исходные уравнения с порядком 0(т + hk), где h = max(hr, hz), m = 2 при a = 0.5 + 0(т) и m = 1 иначе, и безусловно устойчива при a > 0.5. Как следует из вида расщепленных операторов, подобно одномерному случаю она реализуется на дробных шагах скалярными прогонками, а на этапе корректора — явно.

2

4. Численные расчеты

По предложенной разностной схеме проведены расчеты задачи о распространении облака плазмы (7 = 5/3). В начальный момент в области \Jr2 + z2 <0.1 задан скачок давления, превышающий фоновое давление на три порядка (в силу симметрии задачи рассчитывалась только четверть области):

(га, vr, vz, р, Д., BZ)J=0 = (по, 0, 0, 500, 0, 1)т при л/г2 + z2 < 0.1, иначе (n, vr, vz, p, Br, Bz)T=0 = (1, 0, 0, 0.5, 0, 1)T.

Изучалось распространение облака с течением времени при различных значениях магнитного поля, плотности и температуры при следующих краевых условиях:

z = 0, 0 < r < 1: dn/dz = dvr/dz = vz = dp/dz = dBr/dz = dBz/dz = 0,

z = 1, 0 < r < 1 : n = 1, vr = 0, vz = 0, p = 1, Br = 0, Bz = 1,

r = 0, 0 < z < 1: dn/dz = vr = dvz/dz = dp/dz = dBr/dz = dBz/dz = 0,

r = 1, 0 < z < 1 : n = 1, vr = 0, vz = 0, p = 1, Br = 0, Bz = 1.

Численное решение задачи находилось по схеме предиктор-корректор (13) с первым порядком при а = 0.505 на расчетных сетках, содержащих 100 х 100, 200 х 200 узлов. На рис. 1 приведены распределения плотности и давления на момент времени £ = 0.3

по

п0=1000

93.4п

0.5 хС' 0.5

г 0.4 V 0.4 г

0.3 х^ 0.3

0.2 х^ 0.2

0.1 0.1 00

2.51 р

1.5-

11

' -С" "Ч •• 1

0.9

. 0.8

^^0.5 0.60.7 г 0.4 Чх х<^ 0.4 г

0.3 V. Х^ 0.3

0.2 ^/х 0.2

0.1 0.1 00

П0=100

4 32-

ш

1

0.5 ',,. :: х'Г 0.5

г 0.4 'V 0.4 Г

0.3 0.3

0.2 0.2 0.1 0.1 00

0.1 ту 0.1

00

«„=10

1.5] »<>-';>

1-< 0.50 1

0 0

0.1 0.1 00

Рис. 1. Разлет плазмы при Ма = 0-

расширяться. Формируются две волны плотности. Внутри облака образуется область с пониженной плотностью порядка плотности фоновой плазмы и давлением, которое на несколько порядков ниже давления фоновой плазмы. Чем меньше п0, тем быстрее происходит разлет плазмы и быстрее уменьшается плотность в центре области.

«0=1000

200-

100-

п

Ль:

1

0.9

0.8

0.7

0.9

0.8

0.7

Г 0.5 0.4 Г

0.3

0.2 'V 0.2 0.1 ^ 0.1 00

р

11

0.9

0.9

0.8

0.6

0.5 ^е!:', '' 0.5

г 0.4 к 0.4 Г

0.3 0.3

0.2 0.2 0.1 ^ 0.1 00

0.8

«0=100

1 1

о о

о о

1 1

1.5] Р(г'г) 10.5-

о.

о о

0.5 X,,. шЖк п1 0.5

г 0.4 Г х^' 0.4 г

0.3 ^ШШш0.3 0.2 0.2 0.1 0.1 00

10

Рис. 2. Разлет плазмы в магнитном поле.

Рис. 3. Разлет плазмы при различных значениях Ма-

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

На рис. 2 представлены распределения плотности и давления для различных значений п0 на момент времени £ = 0.15 при МА = 5 Видно, что чем больше п0, тем медленнее происходит разлет плазмы и уменьшается плотность плазмы в центре области.

На рис. 3 приведены распределения плотности и давления на момент времени £ = 0.2 при п0 = 1000 для МА = 5, МА = 30 соответственно. С возрастанием магнитного поля (увеличением Ма) разлет плазмы в радиальном направлении становится меньше, чем в продольном, при МА = 30 он практически прекращается, облако разделяется па две части, которые разлетаются вдоль оси г. При этом амплитуды плотности и давления сильно возрастают. Разлет в направлении оси г не зависит от магнитного поля.

Результаты проведенных расчетов правильно передают физику исследуемого явления, качественно и количественно совпадают с полученными в [6]. В отличие от [6] предложенная разностная схема безусловно устойчива, позволяет проводить расчеты с большими шагами по времени.

Полученные результаты позволяют сделать вывод об эффективности предложенного алгоритма и возможности его применения для решения более сложных задач с учетом реальных эффектов.

Список литературы

flj Яненко H.H. Математика, механика: Избранные труды. М.: Наука, 1991. 416 с.

[2] Ковеня В.М., Яненко H.H. Метод расщепления в задачах газовой динамики. Новосибирск: Наука, 1981. 304 с.

[3J Брагинский С.И. Явления переноса в плазме // Вопросы теории плазмы. М.: Госатомиздат, 1963. Т. 1. С. 183-212.

[4] Ковеня В.М. Разностные методы решения многомерных задач. Новосибирск: Изд-во НГУ, 2004. 145 с.

[5] Ковеня В.М., Козлинская Т.В. Об алгоритме расчета нагрева плазмы электронным пучком // Вычисл. технологии. 2004. Т. 9, № 6. С. 59-67.

[6] Астрелин В.Т., Бурдаков A.B., Губер H.A., Ковеня В.М. Моделирование движения и нагрева неоднородной плазмы // ПМТФ. 2001. Т. 42, № 6. С. 125-140.

Поступила в редакцию 4 апреля 2006 г.