CC BY
49
Вычислительные технологии
Том 11, часть 1, Специальный выпуск, 2006
МЕТОД ПРЕДИКТОР-КОРРЕКТОР ДЛЯ РЕШЕНИЯ ЗАДАЧ МАГНИТНОЙ ГИДРОДИНАМИКИ*
В. М. Ковеня
Институт вычислительных технологий СО РАН, Новосибирск, Россия
e-mail: [email protected] т. в. козлинская Новосибирский государственный университет, Россия e-mail: [email protected]
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¥г + УУ, = О,
дЬ
(1)
/
где \УГ = -
д_
дг
(гПУг )
\
д др дВ
— (ГПУ2г) + к0г^ + Ш\Вг —-дг дг ^ дг
д дВ
— (гпугУх) - М дг \ дг
дд Т- [гугр) + 1р—ГУг дг дг
0
д
—г {ухВг - угВх)
\
(
д_
дг
(пух )
\
ТГ {пугух) - М\в^ я дг „ дг
д , , др л дВг
— (пу1) + кот^ + м
дг
\
дг
д ду
— {УГВХ - ухВг) дг
0
II = (п, пуг, пуг, р, Вг, Вг)т-,1 = 7— 1; Ма = £>о/\/ 4:7тпоЩгПг — альфвеновское число Маха; к0 = Т0/(и0mj), mj — масса иона. Система уравнений (1) замыкается уравнением состояния в виде
р = пТ. (2)
При преобразованиях уравнений (1) используется условие
ЛуВ = 1А (гВг) + ^
г дг дг
0.
(3)
При переходе к безразмерному виду пространственные координаты нормировались на характерную длину время — на Ь/и0, а плотность, температура, скорость и магнитное поле — па их невозмущенные значения п0, Т0, и0, В0. При построении разностной схемы на этапе предиктора использованы уравнения в недивергентной форме. Выбирая в качестве искомых функций вектор / = (п, уг, ух, р, Вг, Вх)т, перепишем (1) в виде
д/
дЬ
+ С / + С / = 0.
(4)
Здесь
С г —
( д
УгТГ
дг 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 = о,
(6)
ип+1 _ ип
---—|- = о
т
аппроксимирует исходные уравнения (5) с порядком 0(тт + Нк), оде т = 2 при а = 0.5 + 0(т), и па этапе предиктора реализуется векторной прогонкой, как следует из вида оператора С. При а > 0.5 она безусловно устойчива. Для получения схем, реали-
С
С = С1 + С2, (7)
где расщепление зададим в следующей форме:
С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,
(8)
т
(9)
аппроксимирует исходные уравнения (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(к Л7рпЛ + ^2ВгпЛВгпЛ)] <+1/4 = _ та(^Лрп + А^ЛВ^. (10)
Его решение находится скалярной прогонкой, после чего явно определяются остальные функции па слое п + 1/4. Подобным образом вычисляются значения функций на слое
, 1 /П п+1/2 „п+1/2
п+ 1/^^^шючая IV из уравнения для Вг , приходим к разностному уравнению вида (10) и решаем его скалярной прогонкой. Все остальные уравнения решаются независимо друг от друга скалярными прогонками. Наконец новые значения ип+1 находятся явно из последнего векторного уравнения схемы (8).
Таким образом, решение системы из шести уравнений при переходе с п-го на (п + 1)-й слой находится шестью скалярными прогонками. Введение расщепления (7) приводит в отличие от нефакторизованной схемы (6) к дополнительным членам лишь в уравнении движения для у*, так как матрица С1ЛС2Л содержит ненулевые элементы лишь в третьей строке. Все другие формы расщепления (7) приводят к большему числу дополнительных членов.
Покажем, что разностная схема (8) для линейных уравнений, полученных из (5) "замораживанием" коэффициентов, безусловно устойчива. Для линейных уравнений разностная схема (8) может быть переписана в виде
/п+1 _ / п
(/ + таС10К)(1 + таСдА)--+ С0Л/га = 0, (11)
т
где индекс "0" в матричных операторах означает, что коэффициенты этих матриц постоянны. Пусть для определенности Л — симметричный оператор. Будем искать численное решение схемы (11) в виде £п = = ¡0Лпегк^н. Подставляя его в (11), получим характеристическое уравнение (при а = 1)
[(1 + гу¿)\ - 1]2 {[(1 + %ув)Л - 1]4 - %В\е2<1ъ [(1 + гу¿)\ - 1]2 Л2+
+ [(1 + %^)Л - 1]2 [(1 + %ус[)Л - %уд\ (В2 + В2 + с2) ¿2Л +
+В?В^ [(1 + %у^)Л - 1] [(1 + %у^)Л - ъуд\ Л(Л - 1)} =0, (12)
где у = у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,
(13)
где С/г'Г = С\т + С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 г.
