CC BY
53
СОЛИТОНЫ В ДВУХЖИДКОСТНОЙ МАГНИТНОЙ ГИДРОДИНАМИКЕ С УЧЕТОМ ИНЕРЦИИ ЭЛЕКТРОНОВ
М.Б. Гавриков, В.В. Савельев, А.А. Таюрский
В работе аналитически и численно исследуется взаимодействие уединенных волн в модели двухжидкостной магнитной гидродинамики с последовательным учетом инерции электронов. Рассматриваются волны с линейной поляризацией магнитного поля. Главным отличием настоящей работы является использование «точных» уравнений, а не того или иного приближенного подхода (модельного уравнения). Численно показано, что уединенные волны действительно являются солитонами, то есть их взаимодействие подобно взаимодействию сталкивающихся частиц.
Ключевые слова: Двухжидкостная магнитная гидродинамика, плазма, солитоны, инерция электронов, дисперсия волн, схема Лакса-Вендрофа.
Введение
Исследование уединенных волн - актуальная задача механики сплошных сред, имеющая важные приложения. Для существования уединенных волн необходимо, чтобы тензор внутренних напряжений среды имел достаточно сложную структуру, гарантирующую помимо нелинейности уравнений динамики среды также наличие дисперсии. В плазменных сплошных средах кроме нелинейности и дисперсии для существования нелинейных волн необходима нелокальная зависимость электрического поля от других параметров плазмы, выражаемая обычно сложной формой обобщенного закона Ома.
Особый интерес представляют уединенные волны, называемые солитонами, которые при взаимодействии с другими такими же уединенными волнами сохраняют свои форму, скорость и другие характеристики.
Исследование солитонов традиционно базируется на модельных уравнениях [1-3], список которых возглавляют уравнение Кортевега-де Фриза (КдФ), нелинейное уравнение Шредингера (НУШ), синус-уравнение Гордона (СГ). Для плазменных сред этот перечень значительно расширился за счет уравнений Захарова [4], Кадомцева-Петвиашвили [5], нелинейных альфвеновских и магнитозвуковых волн [6] и др.
В настоящей работе численно исследуются уединенные волны в двухжидкост-ной полностью ионизованной квазинейтральной плазме с полным учетом инерции электронов. Эти волны были впервые найдены в работе [7], в которой получены точные формулы для нелинейных бегущих волн в холодной плазме с учетом инерции электронов. Они характеризуются линейной поляризацией магнитного поля -вектор магнитного поля в волне все время лежит в поперечной плоскости и меняется только по величине, сохраняя неизменным свое направление. Принципиальное отличие настоящей работы от аналогичных исследований состоит в отказе от идеологии модельных уравнений. Таким образом, точные уравнения двухжидкост-ной гидродинамики плазмы [8], выражающие фундаментальные законы сохранения массы, энергии, импульса электронов и ионов и законы электродинамики служат и для нахождения уединенных волн, и для исследования их взаимодействия. В этом смысле приведенные в работе результаты являются полезными при анализе результатов, полученных другим путем.
Как показано в работе, уединенные волны с поляризацией магнитного поля являются солитонами - при столкновении уединенных волн, набегающих друг на друга, они подобно материальным частицам, сохраняют форму, скорость, амплитуду и т.д., при этом сам процесс столкновения имеет конечную длительность.
1. Основные уравнения двухжидкостной электромагнитной гидродинамики
Двухжидкостная гидродинамика плазмы исходит из представления об электронах и ионах как двух взаимно проникающих жидкостях, распределенных по всей области течения. В отсутствие диссипаций для полностью ионизованной двухкомпо-нентной квазинейтральной плазмы уравнения двухжидкостной гидродинамики плазмы [8] можно записать в следующей одножидкостной форме, составляющей содержание уравнений электромагнитной гидродинамики (ЭМГД) [9]:
^ + ё1ури = 0, ^ +Б1УП = 0, (1)
дЬ дЬ
^ + УРе^у-и = ^^ (Р£) , + ^1уи = — А,еРу-1,]У (Р с2КгКе ^ 1 Гтт „п ттт й д
(2)
Е +-—го^Е = — Н] + -Б1у1, — = — + - -V, (3)
4пр с р аЬ дЬ
где тензоры П и 1 вычисляются по формулам
П = П + ПР + Пс,
Пн = рии + р1з, П = Н2 1з — НН, Пс = ^, (4)
8п 4п р
1 = (Ке — Кг) (ПР + Пс) + (кеРг — КРе) 1з + ЬК (— + jU) .
В формулах (1)-(4) 13 - единичный трехмерный тензор; р = ре + рг - суммарная плотность; и = (р^г + р^е)/р - массовая гидродинамическая скорость плазмы; ке = ше/в; кг = шг/^в), где Z - кратность заряда ионов. Кроме того, для простоты предполагается, что электроны и ионы - идеальные политропные газы с общим показателем адиабаты у. Система (1)-(4) замыкается уравнениями электродинамики Максвелла для квазистационарного электромагнитного поля
1 9Н
—— + г^Е = 0, ёгеН = 0,
с аЬ
(5)
4л.
гоШ = — ^
с
Таким образом, (1)-(5) - замкнутая система уравнений относительно величин р, рг, ре, И, Н, Е. Используя ее решение, гидродинамические параметры электронов и ионов восстанавливают по формулам
V; = И + —= И - —j, Рг = кр, Ре = кер, к = ке + кг. (6) р р к к
Уравнения (1)-(3) математически эквивалентны [9] законам сохранения массы, энергии, импульса для электронов и ионов. Система (1)-(5) позволяет рассматривать плазму как сплошную среду, термодинамическое состояние которой задается тремя параметрами (скажем р, рг, ре). Согласно (3), (4), имеется два принципиальных отличия от МГД-теории [10]. Во-первых, тензор плотности потока импульса П содержит дополнительное слагаемое Пс (его можно назвать «токовым»), а в классической МГД - П = Пн + Пр. Во-вторых, кардинально изменился обобщенный закон Ома (3), где появилась добавка пропорциональная г^г^Е. В результате, в ЭМГД поле Е нелокально зависит от остальных параметров плазмы; в частности, для его нахождения необходимо решать краевую задачу для системы (вырожденных) эллиптических уравнений (3). Наконец, тензор W в обобщенном законе Ома (3), состоящий из так называемых холловских членов, также отличается от классического и хол-ловского случаев дополнительными слагаемыми (в классической МГД [10] W = 0, в холловской МГД [11] W = —кг(Пр + ре13)). Уравнения классической МГД формально получаются из системы (1)-(5), если в ней удалить все слагаемые, содержащие р в знаменателе, а уравнения холловской МГД - если в (1)-(5) формально положить ке = 0.
Уравнения динамики холодной плазмы получаются из (1)-(5) при рг = ре = 0. Тогда уравнения (2) выполнены тождественно и упрощаются выражения для тензоров П и W
П = рИИ, W = (ке — кг) (ПР + Пс) + кекг (И + jИ) . (7)
2. Уравнения бегущих волн в холодной ЭМГД-плазме
Рассмотрим решения ЭМГД-уравнений (1)-(5) в случае холодной плазмы, зависящие от Ь, г в комбинации 9 = гк — аЬ, где к- единичный вектор, а- константа. Такие решения обычно называются плоскими бегущими волнами, а- фазовой скоростью волны, к- направлением распространения волны. Нас интересует взаимо-
действие уединенных бегущих волн, все параметры которых имеют конечные и равные предельные значения при 9 ^ Параметры произвольной бегущей волны и(6), р(9), Н(9), Е(9) при фиксированных к и а находятся подстановкой в систему (1)-(5), где положено = ре = 0. После несложных преобразований получаются следующие уравнения для поперечного магнитного поля Н^(9) и продольной скорости и(9) = Пц (9) — а в системе отсчета движущейся волны:
тт2
Ju + Hi = D, 8п
u —
Н\\ 4 4nJ
H
1
С2 AjAe d ( H \ сН\\ )
-——u— uH I + -—7 (Aj — Ae) u 4nJ d6 V xJ 4nJ v j '
(8)
k, H
i
+ q = 0,
где точка над буквой означает дифференцирование по 6, каждый вектор раскладывается вдоль (У) и поперек (±) направления распространения волны k, q±k, а J = 0, D > 0 - произвольные константы интегрирования. Кроме того, в бегущей волне всегда Ну (6) = const. Остальные параметры волны выражаются через u(6), Hx(6) по формулам
р = '—, U = u + a, u
u± = q« +
1 J 4nJ
H
,
E
=
H
k, q + -j qo — aHi
E
—^[Ux, Hi] — ^Ae H k,
с 1 1 1 8пр d6 '
=0' ji = 4n
k, H
_L
(9)
1
С
где qo^k еще одна константа интегрирования. Из = 0, в частности, следует, что в бегущей волне электроны и ионы вдоль направления волны к двигаются с общей гидродинамической скоростью (9), которая, как и поперечное электрическое поле Е±(9), зависит от фазовой скорости а. Система (8), очевидно, сводится к автономному дифференциальному уравнению 2-го порядка относительно Н^(9). В частном случае в двухжидкостной форме система (8) получена в [12], в общем случае - в работе [13].
В этой работе исследуются уединенные волны специального типа, являющихся решением системы (8) для Нц = 0 вида Н^(9) = Н(9)ео, где во^к - единичный вектор, q = дво. Иными словами, ниже рассматриваются только бегущие волны, в которых вектор магнитного поля Н меняется только по величине, имея при этом фиксированное направление в поперечной плоскости. Для таких волн, согласно системе (8), функции Н(9), и(9) ищутся из уравнений
Н
Ju + — = D, 8п
с2 A Ae dfdH
Hu---—— u— u^— + q = 0
4nJ d6 V d6 4
где q - произвольная константа.
(10)
3. Уединенные волны в покоящейся плазме
Найдем все уединенные волны, являющиеся решением системы (10). Не теряя общности, достаточно ограничиться волнами в покоящейся плазме. Приведем систему (10) к безразмерному виду, нормируя магнитное поле и плотность на их предельные значения на бесконечности в неизвестной пока уединенной волне: Н = Н—, ро = р—. Примем альфвеновскую скорость у0 = Н— /^/4лр— за характерное значение скорости и считаем Ь0 = Ь0/у0, где Ь0 - характерное значение длины. Тогда из (9), (10) имеем, учитывая неподвижность плазмы на бесконечности (и— = 0)
Н 2 Н 2
Т п т I Н^> 2 I Н0
п— = —а, . = р—п— = —р0а, V = ..и— + —— = р0 а + —.
8п 8п
Тогда в безразмерном виде система (10) перепишется так:
1 \ Н2
п = — а +--+--,
2а/ 2а
„2 а / ан )
?2п0ё па9 +аНп+9 = 0'
(11)
где ^ = (ск1/2ке/2)/^\/4лр0Ь0) - безразмерный параметр, а - безразмерная фазовая скорость, д - безразмерное значение константы д из (10).
Поскольку величина п(9) = . /р = —ар0/р(9) сохраняет знак, то можно перейти к независимой переменной в, d/ds = п(9)й/й9 и свести систему (11) к автономному уравнению 2-го порядка на прямой
+ д + Н — И {а2 + 1) =0. (12)
Уравнение (12) несложно проинтегрировать [14]. Введем потенциальную функцию
(13)
/Н3 1 Н4 Н 2
[д + -у- — Н (а2 + 2)]аН = дН + — — —(1 + 2а2).
Тогда решение имеет формальный вид [14]
в = аН , (14)
V ^2(г — и(Н))' ( '
где е- произвольная константа. Фактически решение собирается из дуг, задаваемых формулами (14). Тип решения и правило «сборки» зависят от вида потенциальной ямы функции (13), содержащей решение. Уединенная волна получается, если один (и только один) из краев потенциальной ямы является локальным максимумом. В принятой системе единиц, таким краем должна быть точка Н = 1. Таким образом, должны быть выполнены условия
е = и(1), и'(1) = 0, и"(1) < 0.
Отсюда находится константа д = а2, получается неравенство |а| > 1 и вычисляется е = а2/2 — 1/8 = е0. Потенциальная функция и(Н) изображена на рис. 1.
При е = е0 множество и(Н) < е0 распадается на две потенциальные ямы с общим краем Н = 1. Другие края этих ям обозначим через Н1, Н2. В принципе каждая из этих ям может доставлять уединенные волны. Однако уединенные волны в левой яме [Н1,1] не годятся по физическим ограничениям. Действительно, для безразмерных величин ри = —а и значит
1 = — и = (1 + ±_ — Н2
р а V 2а2/ 2а2 Но плотность - положительная величина, так что для всех 9 должно быть справедливо неравенство
1 Н2 (9)
Рис. 1. График потенциальной функции U(H)
1 + 2а2 -
2а2
> 0.
Откуда Н(9) < Нмакс = у/1 + 2а2. Следовательно, необходимо выполнение условий Н2 < Нмакс, |Н1| < Нмакс. Величины Н1, Н2 - это корни уравнения и(Н) = ео, и они легко ищутся, учитывая что один двукратный корень Н = 1 этого уравнения известен. Имеем
(Н — 1)2
U(H) - ее =
[(Я + 1)2 - 4а2] .
Отсюда Н1 = —1 — 2 |а|, Н2 = —1 + 2 |а|. Условие Н2 < Нмакс равносильно неравенству |а| < 2, а условие |Н1| < Нмакс всегда не выполняется. Итак, уединенная волна существует только при
1 < а < 2,
расположена в правой потенциальной яме [1, Н2] и имеет амплитуду Н2 — 1 = = 2(|а| — 1). Явное выражение для уединенной волны получается вычислением интеграла (14). Удобнее сразу найти зависимость 9(Н)
е =
u(H )dH
у/2(ее - U(Я))
= Т2&/
_ Я2 1 + 2а2) - 2а2
dH
(Я - 1)д/4а2 - (H + 1)2
= ±у/4а2 - (H + 1)2-
y/O2-!
ln
2а2 - 1 - H + у/а2-1у/4а2 - (H + 1)2
H1
+ const.
Подберем const так, чтобы при е = 0 выполнялось равенство H(е) = H(0) =
= H2 = 2 |а| - 1. Легко найти значение const = "%а 1па^у/а2 - 1. Откуда получаем окончательное выражение
е=±1
а | у/а-
ln
2а2-1 - H + у/а2-Т^4а2 - (H + 1)2 |а| (H - 1)
-у/4а2 - (H + 1)2
2
а
Полученное выражение (15) для 9(Н) позволяет найти ширину уединенной волны. Воспользуемся тем, что четная функция Н(9) (рис. 2) имеет с точностью до знака единственную (см. ниже) точку перегиба 9о > 0. Проведем в точках перегиба ±9о касательные к графику функции Н = Н(9) до их пересечения с прямой Н = 1, к которой при 9 ^ асимптотически приближается график Н = Н(9). Тогда под шириной уединенной волны А естественно понимать длину отрезка, высекаемого из прямой Н = 1 указанными касательными. Из этого определения, очевидно, следует формула для ширины уединенной волны
Рис. 2. Зависимость Н(9) для различной фазовой скорости а
А = 2
(Нщ + 9о) =2 ((1 — Но)9'(Но) + 9о) , Но = Н(9о).
Точки перегиба ищутся из уравнения Н"(9) = 0. В данном случае проще решить равносильное уравнение 9"(Н) =0 и найти сначала Но. Дифференцируя выражения
для
9'(Н) = П(Н)
У ' ^^2(ео — и (Н))
и приравнивая его нулю, находим единственный корень Но = (2а2 + 1)/3 (в частности, есть только одна, с точностью до знака, точка перегиба функции Н(9)). Теперь из (15) находим 9о = 9(Но) > 0 и
Но = (2а2 + 1)/3, 9о = ^
1 _
|а| 1 ^а2—1
2 , 2 + ^4—О2 2
1п
|а|
— у/ (а2 — 1)(4 — а2)
В итоге получается следующая явная зависимость ширины уединенной волны от фазовой скорости:
21
2 + а2
А = 2((1 — Но)9'(Но) + 9о) = , , ГЬ2 _ <| л/4—а2 + а21п ~ ' " }. (16)
1а1^ а2 — 1
|а|
В частности, А = А(а) монотонно убывает от при |а| = 1 до 0 при |а| =2 и пропорциональна 1 .
4. Методика численного моделирования уединенных волн
Выше было установлено, что в покоящейся однородной с плотностью плазме в однородном магнитном поле с напряженностью Н= 0 поперек этого поля могут распространяться уединенные волны (15) с фазовой скоростью а^л, где
уа = Н^/альфвеновская скорость с амплитудой 2(|а| — 1)Нте и шириной А(а), вычисляемой по (16) для любого 1 < |а| < 2. В частности, волны с большей амплитудой более узкие и двигаются быстрее более широких волн с меньшей амплитудой. В уединенных волнах продольное магнитное поле равно нулю, а поперечное меняется только по величине, сохраняя неизменным свое направление. Тогда, как следует из соотношений (9), поперечное электрическое поле будет тоже меняться только по величине, причем Е^,те = 0, а поперечная скорость постоянна.
Рассмотрим четыре задачи о взаимодействии волн. Нас интересует, что происходит при:
• движении двух уединенных волн одинаковой амплитуды навстречу друг другу;
• движении двух уединенных волн разной амплитуды навстречу друг другу;
• набегании уединенной волны с большей амплитудой на уединенную волну с меньшей аплитудой в предположении, что волны двигаются в одну сторону;
• распаде начального возмущения, локализованного в пространстве и, в частности, может ли начальное возмущение порождать пакет уединенных волн.
При этом принципиально важно, что взаимодействие уединенных волн определяется не модельными уравнениями (типа КдФ, НУШ, СГ и пр.), а полной системой уравнений двухжидкостной гидродинамики плазмы в ЭМГД-форме, выражающей фундаментальные законы сохранения массы, энергии, импульса электронов и ионов и законы электродинамики Максвелла. Для математического исследования взаимодействия уединенных волн ЭМГД-система решается численно.
Направим ось х вдоль направления распространения волны к, ось г - вдоль направления поперечного магнитного поля ео, тогда поперечное электрическое поле направлено вдоль оси у, а функции их, Нх, Еу, р подчинены ЭМГД-уравнениям (1)-(5), которые с учетом плоской геометрии и приближения холодной плазмы примут вид
др + д рЦх дЬ дх
= 0,
дри. + д (ри2 + НЛ =о
дь + дх ^рих + 8п у
1 дНг + дЕу „
с дь
дх
(17)
Е _ с% Хе д2Еу
у 4пр дх2
ихНг сАДе д / дНх
и х
4пр дх
дх
Система (17) решается на прямой —то < х < для Ь > 0 с начальными условиями при Ь = 0, которые уточняются ниже, и граничными условиями
р(±то) = р^, Нх (±то) = НЕу (±то) = 0, их(±ю) = 0.
(18)
Уравнения (10) можно получить более прямым путем, если искать решения системы (17), зависящие от Ь и х в комбинации 9 = х — аЬ. Для существования в плазме нелинейных, в частности, уединенных волн необходимо наличие дисперсии. В общем
с
случае наличие дисперсии вытекает из дисперсионного соотношения для ЭМГД-системы (1)-(5), которое имеет вид [15]
1+ -I -
Юр
k2
2k2
ю2 — c2k2
+
ю2
í ck\' \юр)
1+ — -
k2
ю2
ck
юр
- - Л2
v\\k\'
ю
Л = (АДе)1/2 - (^e/li)1/2 , ю2 =
4про
lile
Но
v = VA =
\/4Про'
le С2 + kc2
d =
le + l
■гЧ
ce2 =
e
YP0
p0 :
c2 = Ypi c p0 ■
(19)
Здесь р0, р0, ро = р0 + р0, Но, и0 = 0 - невозмущенное однородное состояние плазмы, к - направление распространения линейной волны ехр[г(гк — ю¿)] и все векторные величины раскладываются по компонентам к. В частности, из (19) с учетом приближения холодной плазмы = р0 = 0 и поперечного характера волны = 0 получается дисперсионное уравнение для системы (17)
1+1 2 -ю
VAk
ю
0,
2
v2 _ Нс
VA
4лрс
которое приводит к нелинейной зависимости (то есть к дисперсии) ю от к
ю^) = ±
VAk
\Jl + (ck/Юр)2
(20)
Отклонение от линейного закона ю(к) = ±^лк, справедливого в МГД, определяется скиновой или дисперсионной длиной с/юр. Для длинноволновых возмущений имеем (ск/юр)2 ^ 1 и дисперсией, как следует из (20), можно пренебречь. Заметим, что наличие дисперсии является необходимым, но не достаточным условием существования уединенных волн. Например, в холловской МГД [11] дисперсия есть, но плоских нелинейных волн нет, поскольку обобщенный закон Ома в холловской МГД определяет локальный характер зависимости электрического поля от остальных параметров плазмы.
Для численного решения задачи (17), (18) перепишем ее в безразмерном виде (дополнительно полагая Е0 = у0Н0/с и обозначая их = и, И = И, Еу = Е)
ди Ю (и, И) п . . . гг.
К + =0' и = (и1,и2) = (р, ри),
дН дЕ u2 И2. . тт тт2 И2.
-Ж + аХ = 0' f = (u2, ui + ^) = (р^' pu2 +
(21)
Е - b(u)
д2Е дх2
= D(u,H), -ж <х< ж, t > 0
с граничными условиями
р(±ж) = 1, И (±ж) = 1, и (±ж) = 0, Е (±ж) = 0.
(22)
2
V
V
V
0
2
Здесь b(u) = u1 = р > 0, а оператор D(u, H) имеет вид
D(u, H) = U2H/U! - (?2/u!) U^ H = UH — ^ иН .
dx \u1 dx ) р dx \ dx )
В записи (21) исследуемая система распадается на «гидродинамическую» часть, записанную в дивергентном виде, и «электродинамическую» часть, содержащую обобщенный закон Ома и закон Фарадея. Рассмотрим численный метод решения системы (21), (22), основанный на двухшаговой схеме Лакса-Вендроффа [16] для гидродинамических уравнений. Выберем две равномерные разностные сетки на прямой с целыми Xk = kh, —ж < k < +ж и дробными (полуцелыми) xfc+1/2 = (k + 1/2)h, —ж < k < +ж узлами, где h > 0 - шаг сетки, k - целое. Аппроксимируем значения неизвестных функций р, U, H, E в узлах целой сетки сеточными функциями, которые и подлежат нахождению в каждый момент дискретного времени. Переход с временного слоя t = t0 на слой t = t0 + т, {pk, Ц0, Hj0, E^} ^ {pk,U1 ,H1,E1} происходит в два этапа.
1. Вычисляем вспомогательные сеточные функции |р1+21/2, H1^1/2,
<2/2}, которые интерпретируются как приближенные значения неизвестных функций р, U, H, Е в дробных узлах в момент времени to + т/2:
u1/+21/2 — Uk+1/2 + f (uk+1 ,H0+1) — f (u0 ,H0) =o т/2 + h '
H1/2 _ H0 E0 E0 u0 + u0 tT0 + TT0
Hk+1/2 Hk+1/2 + Ek+1 — Ek = 0 0 = uk + uk+1 я 0 = Hk + Hk+1
т/2 + h ' uk+1/2 = 2 ' Hk+1/2 = 2 :
E1/2 _ 2E1/2 + E1/2
—Ek+21/2 — b(uki21/2)^-h+21/ k"1/2 = D(u1/2, H 1/2)k+1/2
с простейшей аппроксимацией дифференциального оператора D(u1/ ,
H1/2)
D(u1/2 Н 1/2)k+1/2 = Uk+21/2Hk+21/2 —
fc2 1 и1/2 +и 1/2 я1/2 — я1/2 и 1/2 + TT1/2 H1/2 — я1/2
§2 l JUk+3/2^Uk+1/2 Hk+3/2 Hk+1/2 Uk+1/2^ Uk—1/2 Hk+1/2 nk-1/2\
р1/2 h 1 2 h 2 h
^+1/2
2. Вычисляем искомые сеточные функции {рк ,Н1,Е1}, зная вспомогательные сеточные функции, найденные на первом этапе:
и1 — и0 + ГК^/И^) — £(и1-21/2' Н-1/2) = 0
Н
Н — И + Е1/21/2 — Ек2-1/2 =0,
т Н '
Е1 — 2Е1 + Е1 Е1 Ь(и1) Ек +1 2Е к + Е к-1 п(и1 Н 1) Ек — Ь(ик)-Н2- = П(и , Н
с естественной аппроксимацией дифференциального оператора ^(и1, И1)
п(„1 Н1) и 1Н1 1 / и1+1 + и1 Н1+1 — Нк и1 + и1—1 Н1 — И-1 \ я(и ,Н )к = ^^ — рк^-2--Н---2--Н-/ •
пРи этом ь(ик+21/2) = ?2/рк+21/2, ь(ик) = ¥/р1 .
Реальный расчет ведется на конечной сетке в конечной области [0,£] в целочисленных узлах хк = кН, Н = £/N N > 1 - натуральное) с учетом граничных условий (22). На первом этапе вспомогательные величины и1/2, Н1/2, Е1/2 вычисляются во внутренних дробных узлах хк+1/2, 0 < к < N — 1. На втором этапе расчет и1, Н1 ведется во всех целых узлах хк, 0 < к < N, а расчет Е1 - во всех внутренних целых узлах хк , 0 < к < N. Недостающие (при расчете Е1/2 на первом этапе и и1, Н1 - на втором) значения р1/2, и1/2, Н1/2 в заграничных узлах считаются равными константам в соответствии с граничными условиями (22). Сеточные функции Е1/2 и Е1 ищутся прогонкой с нулевыми, согласно (22), условиями на границе: Е-/2/2 = Е1/_+1/2 = 0, Е0 = Е^ = 0 . Наконец, шаг по времени т выбирается из условия Куранта
т = к-
тах 0< г<м
«0
+ £
где £ > 0 - малая добавка, 0 < к < 1 - коэффициент запаса, «Л г = Н0^ ^4лр0.
5. Результаты расчетов
Рассмотрим эволюцию двух солитонов равной амплитуды, двигающихся навстречу друг другу с безразмерной скоростью а = 1.5. Для этого система (21), (22) (здесь и ниже ^ = 1) решается численно двухшаговым методом Лакса-Вендроффа с начальным условием для магнитного поля
Н(0,х) = Н1.б(х — 45), 0 < х < 75,
Н(0,х) = Н—1.5(х — 105), 75 < х < 150
на отрезке [0,£] с £ = 150, где На(х) - функция, получающаяся обращением ± ветвей (15) и заданная для всех вещественных х. Начальные значения остальных параметров вычисляются по формулам раздела 3
гт Н2 ~ 1
и =-, р =
2а ' к
Е = а(Н - 1).
1 + 2а2 / - 2а2
-1
(23)
Рис. 3. Взаимодействие уединённых волн с одинаковыми амплитудами и фазовыми скоростями: а = 1.5 и а = —1.5 в различные моменты времени: ¿0 - до взаимодействия, ¿1 - во время взаимодействия, ¿2 - после взаимодействия
На рис. 3 изображены результаты расчета профилей напряженности магнитного поля Н в различные моменты времени. Как видно, сначала (£ = ¿0) уединенные волны двигаются со скоростью а = 1.5 навстречу друг другу, не меняя форму. Затем начинается слияние двух волн в одну; в результате возникает единая волна с амплитудой большей, чем суммарная амплитуда уединенных волн до слияния (при Ь = ¿1 процесс слияния еще не достиг своего пика - на горбе волны видна впадина). Далее образовавшаяся единая волна распадается и рождает два солитона, расходящиеся в противоположные стороны с теми же скоростями и амплитудами, что и до слияния (Ь = ¿2). Численно находится время взаимодействия двух волн тс = 0.36. Как показал расчет, в хвосте каждой из расходящихся уединенных волн возникают затухающие колебания малой амплитуды, заполняющие пространство между разбегающимися после столкновения солитонами.
Пусть навстречу друг другу двигаются уединенные волны с разными амплитудами и фазовыми скоростями: слева направо двигается уединенная волна с амплитудой 1 и фазовой скоростью а = 1.5, справа налево - волна с амплитудой 0.2 и фазовой скоростью а = -1.1. Начальное магнитное поле для этих волн
Н(0, х) = Н1.5(х - 80), \х - 80\ < 35,
Н(0, х) = Н-1.1(х - 155), \х - 155\ < 35
и Н(0, х) = 1 в остальных точках счетной области [0, £] с £ = 300. Начальные значения остальных параметров вычисляются по (23).
На рис. 4 изображены расчетные профили напряженности магнитного поля Н в различные моменты времени ¿о < ¿1 < ¿2. Как видно, основные закономерности взаимодействия уединенных волн такие же, как и в предыдущем случае равных амплитуд. Только уменьшились время столкновения тс = 0.23 и амплитуды затухающих колебаний в хвосте каждой из расходящихся после взаимодействия уединенных волн. Однако максимальная амплитуда при слиянии волн, по-прежнему, больше суммарной амплитуды волн до взаимодействия.
н
н
-1.6
-1.2
о
100
150
х
300
400
500 х
Рис. 4. Волны с разными амплитудами и фазовыми скоростями: а = 1.5 и а = —1.1 , движущиеся навстречу друг другу
Рис. 5. Волны с разными амплитудами и фазовыми скоростями: а = 1.5 и а = 1.1, движущиеся друг за другом
Далее рассмотрим случай, когда волна с большей амплитудой 1 догоняет волну с меньшей амплитудой 0.2. Тогда начальное условие для магнитного поля Н имеет вид
и Н(0, х) = 1 в остальных точках счетной области [0, £] с £ = 600. Начальные значения р, и, Е вычисляются по (23). На рис. 5 представлены профили напряженности магнитного поля Н в разные моменты времени. Из рисунка следует, что уединенная волна с большей амплитудой догоняет (Ь = ¿о) волну с меньшей амплитудой и сливается (Ь = Ь\) с ней, образуя единую волну с амплитудой меньшей, чем суммарная амплитуда уединенных волн до взаимодействия. Затем возникшая единая волна раздваивается (Ь = Ь2) и порождает две уединенные волны с такими же параметрами, которые были у волн до взаимодействия. Однако теперь волна с большей амплитудой располагается правее волны с меньшей амплитудой. Иными словами, волна с большей амплитудой обогнала волну с меньшей амплитудой. Как показал расчет, затухающие колебания в хвосте уединенных волн после столкновения не возникают.
Наконец, рассмотрим случай, когда начальное возмущение магнитного поля имеет вид гауссова распределения Н = Н0 ехр[—(х — х0)2/О] + 1, для которого начальные возмущения и, Е вычисляются по (23), а р = 0.5а2[(1+2а2)—Н2]_1+0.75 (хотя можно и р вычислять по (23)). На рис. 6 представлены расчетные профили магнитного поля в разные моменты времени для
и Н(0, х) = 1 в остальных точках счетной области [0, £] с £ = 300, причем начальные условия для р, и, Е вычисляются при а = 1.5.
Как следует из рис. 6, происходит рождение нескольких (в данном случае пяти) уединенных волн разной амплитуды, которые разлетаются вдоль оси х, так что расстояние между соседними уединенными волнами увеличивается. Любопытно,
Н(0,х) = Н1.5(х — 200), \х — 2001 < 25, Н(0,х) = Н1.1 (х — 250), \х — 250\ < 25
Выше приведены расчетные профили только для магнитного поля Н. Однако аналогичная картина профилей имеет место и для остальных параметров плазмы.
что пики возникающих уединенных лежат на одной прямой. Из расчетов следует, что далеко не каждое начальное возмущение распадается на соли-тоны. Критерий подобного распада, по-видимому, отсутствует, и этот вопрос требует дальнейшего исследования.
100 200 х
Рис. 6. Эволюция начальных данных в виде гауссова импульса
Заключение
Общий качественный вывод проведенного исследования состоит в том, что рассмотренные в работе уединенные волны подобно материальным частицам могут взаимодействовать друг с другом, сохраняя после взаимодействия свои характеристики, а некоторые плазменные возмущения, локализованные в пространстве, могут распадаться на уединенные волны.
Подобные эффекты типичны для решений модельных уравнений (КдФ, НУШ и т.д.) и обычно объясняются вполне интегрируемостью и наличием бесконечного числа первых интегралов уравнений [2] и последующим применением достаточно рафинированной математической техники (метод обратной задачи рассеяния, преобразования Бэклунда, уравнения Лакса, цепочки Тода и пр. [2, 3]). Однако в настоящей работе эти явления были обнаружены численно на основе более строгой полной системы уравнений двухжидкостной гидродинамики холодной плазмы, то есть законов сохранения массы, импульса и уравнений Максвелла. Данная система не является полностью интегрируемой, чем обусловлены неупругие эффекты при взаимодействии попутных волн, а также ограниченный диапазон амплитуд и скоростей, с которыми могут распространяться уединенные волны.
Важно подчеркнуть, что аналогичные явления были обнаружены и в некоторых других гидродинамических моделях в физике плазмы. Это, например, ионно-звуковые волны [17,19], волны пространственного заряда [18,20]. В указанных работах отмечался упругий характер взаимодействия попутных волн, хорошо описывающийся уравнением КдФ. Встречные же волны всегда взаимодействуют неупруго, но с образованием одиночного пика. Осцилляции между расходящимися после взаимодействия встречными волнами наблюдались и экспериментально в [21]. Некоторые из этих особенностей (неупругие эффекты, образование одиночного пика для встречных волн) в настоящей работе обнаружены численно на основе полной системы уравнений двухжидкостной гидродинамики холодной плазмы, то есть законов сохранения массы, импульса и уравнений Максвелла.
Авторы выражают благодарность проф. Г.Г. Малинецкому за интерес и внимание к работе. В проведении численных расчетов принимал участие студент МИФИ Д.В. Ларионов.
Работа выполнена при поддержке РФФИ (проекты № 09-01-00181, № 08-01-00299).
Библиографический список
1. Захаров В.Е., Манаков С.В., Новиков С.П., Питаевский Л.П. Теория солитонов. Mетод обратной задачи I Под ред. С.П. Новикова. M.: Наука, 1980.
2. Тахтаджян Л.А., Фадеев Л.Д. Гамильтонов подход в теории солитонов. M.: Наука. Гл. ред. Физматлит, 198б.
3. Солитоны в действии I Под ред. К. Лонгрена и Э. Скотта. M.: M^, 1981.
4. Захаров В.Е. Коллапс ленгмюровских волн II ЖЭТФ. 1972. Т. б2. С. 1745.
5. Кадомцев Б.Б., Петвиашвили В.И. Об устойчивости уединённых волн в слабо-диспергирующих средах II Докл. АН СССР. 1970. Т. 192. С. 753.
6. Mio K., Ogino T., Minamy K., Takeda S. Modified nonlinear Scrodinger equation for AlfVen waves propagating along the magnetic field in cold plasma II J. Phys. Soc. Japan. 197б. Vol. 41. P. 2б5.
7. Adlam J.H., Allen J.E. II Phil. Mag. 1958. Vol. 3. P. 448.
8. Брагинский С.И. Явления переноса в плазме II Вопросы теории плазмы. Вып. 1. M.: Атомиздат, 19б3. С. 183.
9. Гавриков М.Б., Сорокин Р.В. Однородные деформации двухжидкостной плазмы с учётом инерции электронов II Изв. РАН MЖГ. 2008. Т. б. С. 15б.
10. Куликовский А.Г., Любимов Г.А. Mагнитная гидродинамика. M.: Физматгиз, 19б2. 2-е изд., M.: Логос, 2005. 328 c.
11. Морозов А.И., Соловьёв Л.С. Стационарные течения плазмы в магнитном поле ИВопросы в теории плазмы (Под ред. MA. Леонтовича). M.: Атомиздат, 1974. Вып. 8. С. 3.
12. Saffman P.G. Propagating of a solitary wave along a magnetic field in a cold collision-free plasma II J. Fluid Mech. 19б1. Vol. 11. P. 1б.
13. Гавриков М.Б. Апериодические колебания холодной плазмы. Препринт № 33. M.: ИПM им. M3. Келдыша АН СССР, 1991. 28 с.
14. Арнольд В.И. Обыкновенные дифференциальные уравнения II Ижевск: ИРТ, 2000.
15. Гавриков М.Б. Линейные волны в нерелятивистской магнитной гидродинамике. Препринт № 199. M.: ИПM им. M3. Келдыша АН СССР, 1988. 28 с.
16. Роуч П.Вычислительная гидродинамика. M.: M^, 1980.
17. Кадомцев Б.Б. Коллективные явления в плазме. M.: Наука, 1988.
18. Рыскин Н.М. Уединенные волны пространственного заряда II Изв. вузов. Прикладная нелинейная динамика. 1994. Т. 2, № 5. С. 84.
19. Рыскин Н.М., Трубецков Д.И. Нелинейные волны. M.: Наука, 2000. C. 150.
20. Трубецков Д.И. Уединенные волны в электронном потоке и юбилей одного уравнения II Соросовский образовательный журнал. 2000. Т. б, № 4. С. 103.
21. Ikezi H., Taylor R.J., Baker R.D. Formation and interaction of ion-acoustic solitons II Phys. Rev. Lett. 1970. Vol. 25, № 1. P. 11.
Институт прикладной математики Поступила в редакцию 24.11.2009
им. М.В.Келдыша РАН После доработки 1.07.2010
SOLITONS IN TWO-FLUID MAGNETOHYDRODYNAMICS WITH NON-ZERO ELECTRON INERTIA
M.B. Gavrikov, V.V. Savelyev, A.A. Tayurskiy
The interaction between solitary waves in two fluid magnetohydrodynamic model of plasma with electron inertia taken into account is investigated analytically and numerically. Waves with linear polarization of a magnetic field are considered. A principal difference of the present work is the using of the «exact» equations, instead of the modeling equations. It is numerically shown, that solitary waves really are solitons, i.e. their interaction is similar to interaction of colliding particles.
Keywords: Two-fluid magnetohydrodynamic, plasma, soliton, electron inertia, dispersion, Lax-Wendroff difference scheme.
Гавриков Михаил Борисович - родился в 1952 году в Москве, окончил Московский университет (мех-мат) в 1975 году и аспирантуру МГУ (1985). В 1990 году защитил диссертацию на соискание ученой степени кандидата физико-математических наук по теории дифференциальных уравнений. Опубликовано свыше 50 научных работ по математической физике и вычислительной математики. Автор (совместно с О.В. Локуциевским) монографии «Начала численного анализа», М.: ТОО «Янус», 1995. 580 с.
125047, Москва, Миусская пл., д.4
Институт прикладной математики им. М.В. Келдыша РАН E-mail: nadyap@cognitive.ru
Савельев Вячеслав Владимирович - родился в 1949 году в Куйбышеве, окончил Московский инженерно-физический институт в 1973 году. Затем учеба в аспирантуре Института прикладной математики АН СССР. В 1977 году защитил диссертацию на соискание ученой степени кандидата физико-математических наук. Имеет более 80 публикаций, посвященных аналитическому и численному исследованию различных задач физики плазмы. В их числе обзоры (совместно с А.И. Морозовым) «О Галатеях - ловушках с погруженными в плазму проводниками». УФН, т. 168, № 11, с. 1153-1194, 1998 г. и «Fundamentals of Stationary Plasma Thruster Theory». Reviews of Plasma Physics, v.21, p.203-391, 2000, Kluwer Academic/Plenum Publishers, N.Y. Consultant Bureau. Лауреат Премии им. И.В.Курчатова (2003).
125047, Москва, Миусская пл., д.4
Институт прикладной математики им. М.В. Келдыша РАН E-mail: ssvvvv@rambler.ru
Таюрский Алексей Александрович - родился в 1984 году в пос. Хандыга, республики Саха (Якутия), окончил Московский государственный технический университет им. Н.Э. Баумана в 2007 году. В настоящее время работает и учится в аспирантуре в ИПМ им. М.В.Келдыша. Результаты дипломной работы опубликованы в 2007 году.
125047, Москва, Миусская пл., д.4
Институт прикладной математики им. М.В. Келдыша РАН E-mail: tayurskiy2001@mail.ru
