"Описание крупномасштабных процессов в бесстолкновительной космической плазме и численное моделирование тонких токовых слоев." тема диссертации и автореферата по ВАК РФ 01.03.03, доктор наук Мингалев Олег Викторович

• Мингалев Олег Викторович
• доктор наук
• 2020

Оглавление диссертации доктор наук Мингалев Олег Викторович

Оглавление

Введение

Глава I. Система уравнений Власова в приближении квазинейтральности

1.1 Введение главы I

1.2 Система Власова-Максвелла, моделирование электростатических эффектов и обобщенный закон Ома

1.3 Система уравнений Власова-Дарвина

1.4 Модификация системы уравнений Максвелла для плазмы в приближении квазинейтральности

1.5 Система уравнений Власова в приближении квазинейтральности

1.6 Схема численного интегрирования

1.7 Выводы главы I

Глава II. Система уравнений Власова в случае наличия замагничен-

ных компонент

11.1 Введение главы II

11.2 Условия замагниченности и их следствия

П.3 Уравнение Власова в дрейфовом приближении

П.4 Система уравнений Власова в случае замагниченных электронов и

незамагниченных ионов

П.4.1 Формальная запись системы уравнений

П.4.2 Основные изменения в схеме численного интегрирования по

сравнению со случаем незамагниченных электронов

П.5 Система уравнений Власова в дрейфовом приближении в случае

полностью замагниченной плазмы

П.5.1 Формальная запись системы уравнений

П.5.2 Схема численного интегрирования системы уравнений для

полностью замагниченной плазмы

П.6 Система уравнений Власова в случае замагниченных электронов и

частично замагниченных ионов

П.7 Заключение главы II

Глава III. Газодинамическое и гибридное описание бесстолкновитель-ной космической плазмы в приближении силового равно-

весия электронов вдоль магнитного поля

111.1 Введение главы III

111.2 Дивергентная форма уравнения потока импульса для протонов в плазме из протонов и замагниченных электронов

111.3 Анализ системы уравнений одножидкостной МГД

111.4 Система уравнений магнитной газовой динамики для замагниченной компоненты

111.5 Система уравнений многокомпонентной магнитной газодинамики для бесстолкновительной космической плазмы

111.5.1 Формальная запись системы уравнений

111.5.2 Схема численного интегрирования

111.6 Гибридное описание бесстолкновительной плазмы в приближении продольного силового равновесия электронов

111.6.1 Формальная запись системы уравнений

111.6.2 Схема численного интегрирования системы уравнений гибридного описания

111.7 Заключение главы III

Глава IV. Модель стационарного пространственно двумерного токового слоя с нормальной компонентой магнитного поля и

кинетическим описанием замагниченных электронов

IV.1 Введение главы IV

IV.2 Постановка задачи

IV.3 Общий вид функции распределения замагниченных электронов . . . 133 IV.4 Система уравнений модели стационарного пространственно двумерного токового слоя

IV.5 Случай распределения Максвелла-Больцмана для электронов

IV.6 Заключение главы IV

Глава V. Численное моделирование стационарного тонкого токового

слоя в хвосте магнитосферы Земли

V.1 Введение главы V

V.2 Силовой баланс и ток изотропных электронов в тонком токовом слое

V.3 Уравнения модели с плоским магнитным полем без шира

V.4 Уравнения модели с широм магнитного поля

V.4.1 Случай изотропных электронов

V.4.2 Случай анизотропных электронов

V.5 Основные особенности численных моделей

У.6 Моделирование симметричного токового слоя без шира магнитного

поля

У.6.1 Конфигурации с изотропными электронами

У.6.2 Конфигурации с анизотропными электронами

У.7 Моделирование токового слоя с широм магнитного поля и изотропными электронами

У.8 Заключение главы У

Приложение. Новые методы численного решения стационарного уравнения Власова

1. Обоснование численных методов

2. Дискретизация в координатном пространстве и в фазовом пространстве

3. Схема итерационного процесса 1-го метода

4. Схема итерационного процесса 2-го метода

5. Экономичный метод расчета фазовых траекторий

Заключение

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

Введение

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

Введение диссертации (часть автореферата) на тему «"Описание крупномасштабных процессов в бесстолкновительной космической плазме и численное моделирование тонких токовых слоев."»

Актуальность проблемы

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

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

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

Однако в современной мировой науке численное моделирование указанных процессов в плазме околоземного космического пространства столкнулось с ря-

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

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

1) адекватность постановки задачи и используемой системы уравнений, которые должны правильно описывать рассматриваемые в модели процессы;

2) корректность используемой в модели методики численного решения в том смысле, что, во-первых, решение используемой в модели "дискретной" системы уравнений является некоторым приближением к искомому "непрерывному" решению исходной системы уравнений модели, и, во-вторых, получаемое в модели решение действительно является точным или приближенным решением "дискретной" системы уравнений модели;

3) правильное задание входных параметров модели, в частности, начальных и граничных условий.

Качество численной модели характеризуется уровнем выполнения всех трех перечисленных условий.

Описывающую плазму систему уравнений можно разделить на две тесно связанные между собой части: на систему уравнений переноса, которая описывает движение каждой компоненты плазмы, и на систему уравнений для полей.

Отметим, что во всех системах уравнений, применяемых для описания достаточно крупномасштабных медленных плазменных процессов, в уравнениях Максвелла отбрасывается ток смещения, что означает переход к мгновенному дальнодействию, в рамках которого пространственное распределение магнитного и электрического полей в области моделирования должно однозначно определяться текущим распределением гидродинамических параметров плазмы и граничными условиями из систем уравнений эллиптического типа, не содержащих частных производных по времени. Именно так происходит определение полей в успешно работающих моделях на основе системы уравнений Власова-Дарвина, в которых уравнения для полей преобразованы к эллиптической форме. Система уравнений Власова-Дарвина получается из исходной системы Власова-Максвелла в результате отбрасывания соленоидальной части тока смещения, что соответствует ис-

ключению электромагнитного излучения.

Главной теоретической проблемой является правильное описание электрического поля. Эта проблема разделяется на две. Первая состоит в правильном описании продольного электрического поля, которое появляется в областях неоднородности плазмы за счет относительно очень малого разделения зарядов и обеспечивает квазинейтральность плазмы, а также существенно влияет на ее движение. Вторая проблема состоит в правильном описании соленоидальной части электрического поля, которая согласно уравнению Фарадея определяется изменениями магнитного поля по времени.

Первым типом численных моделей, которые стали использовать для моделирования рассматриваемых плазменных процессов начиная с 70-х годов прошлого века, являются модели, основанные на численном решении различных вариантов системы уравнений одножидкостной магнитной газовой динамики — так называемые МГД-модели.

В этих системах уравнений указанное выше продольное электрическое поле за счет относительно очень малого разделения зарядов вообще не учитывается. Как будет показано в работе в разделе III.3 главы III, вывод системы уравнений одножидкостной МГД основан на следующих 4-х предположениях: плазма является замагниченной, электроны являются холодными по сравнению с протонами, то есть их давление много меньше протонного, электрическое поле считается ортогональным магнитному, и ортогональная часть скорости ионов считается равной скорости электрического дрейфа.

Анализ характерных масштабов крупномасштабных процессов в магнитосфере и в солнечном ветре приводит к выводу, что условия применимости МГД-подхода в магнитосфере и в солнечном ветре заведомо нарушаются, особенно сильно в наиболее важных областях с сильной пространственной неоднородностью — в различных токовых слоях.

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

В результате в системе уравнений МГД правильное распределение полного тока в плазме не рассчитывается, а для расчета магнитного поля вместо имеющих эллиптический тип уравнений Ампера и Гаусса используется эволюционное уравнение индукции - вмороженности, в котором в качестве коэффициента фигурирует неполная часть ортогональной составляющей скорости плазмы.

Такой подход является аналогом приближения невязкой акустики в газовой ди-

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

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

Вторым типом моделей, которые также применяются для моделирования рассматриваемых крупномасштабных плазменных процессов в очень большом числе работ начиная с 70-х годов прошлого века, являются модели, в которых численно решаются различные варианты системы уравнений гибридного описания плазмы. В этих системах протоны описываются уравнением Власова, причем часто с добавлением модельной "силы трения с электронами", а замагниченные электроны описываются в рамках различных вариантов упрощенной газовой динамики. Чаще всего электроны считаются нейтрализующим безмассовым газом с постоянной температурой по пространству и времени. Электрическое поле определяется из условия силового равновесия электронов с учетом модельной "силы трения с протонами".

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

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

гтч «->

Третьим типом моделей, которые используются для моделирования рассматриваемых крупномасштабных плазменных процессов также в очень большом чис-

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

В этих моделях скорость света занижается не менее чем на 2 порядка, используется очень малое среднее число модельных макрочастиц в дебаевской ячейке ~ 10 — 50, и рассматриваются модельные "тяжелые электроны", которые в 16128 раз легче протонов, в то время как реальные электроны легче протонов в 1836 раз. Численные эксперименты показывают, что для минимально приемлемого уровня воспроизведения электростатических эффектов в самых совершенных неявных версиях метода крупных частиц необходимо использовать в расчетах порядка 100 - 1000 модельных электронов в дебаевской ячейке в пространственно одномерном случае, а в пространственно двумерном и трехмерном случаях требуется ~ 104 — 105 модельных частиц. При этом пространственное разрешение модели в самом грубом из допустимых вариантов должно быть примерно равным дебаев-скому расстоянию электронов. Поэтому при изложенных упрощениях в численной модели крупномасштабных процессов очень хорошим уровнем относительного отклонения от электронейтральности считается 10_1, в то время как в реальности эта функция имеет порядок 10_8 — 10_6, то есть плотность заряда в численной модели завышена как минимум на 5 порядков. В результате в численной модели потенциальное электрическое поле, а значит и скорость электрического дрейфа, завышены на не меньшее число порядков. Это сильно искажает глобальную динамику как ионов, так и модельных "тяжелых электронов", и вносит искажение в распределение плотности тока, которое приводит к искажению магнитного поля, а также к искажению соленоидальной части электрического поля.

Также отметим, что замена реальных электронов на 1-2 порядка более тяжелыми "модельными электронами" существенно меняет свойства и поведение системы на качественном уровне. В результате такие модели лишь на отдаленном качественном уровне могут воспроизводить поведение реальной системы. Поэтому результаты расчетов по этим моделям всегда приводятся в безразмерном виде.

Из сказанного выше следует, что полная система уравнений Власова-Максвелла подходит только для численного моделирования локальных быстрых "излучатель-ных" процессов с пространственным разрешением порядка дебаевского расстояния, которое в магнитосфере имеет характерные значения порядка сотни метров, и с временным разрешением порядка плазменного периода, который в магнитосфере имеет характерное значение ~ 10_4 — 10_3 секунды. Для реалистичного

численного моделировании в области с размером в десятки радиусов Земли крупномасштабных медленно изменяющихся полей с пространственным масштабом изменения на 4 порядка большим характерного дебаевского расстояния и временным масштабом изменения на 5 порядков большим характерного плазменного периода эта система требует огромный объем вычислительных ресурсов, заведомо недоступный в обозримой перспективе.

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

Отметим, что результат осреднения по плазменным колебаниям вдоль магнитного поля известен с 70-х годов прошлого века и хорошо проверен на практике в численном моделировании ионосферных процессов при помощи частично самосогласованных моделей ионосферы. Он состоит в условии квазинейтральности и условии силового равновесия электронов вдоль магнитного поля. В этом случае плазма формально в узлах сетки считается электронейтральной, а объемная плотность продольной силы, действующей на электроны со стороны усредненного по продольным плазменным колебаниям электрического поля, равна продольной компоненте дивергенции тензора давления электронов.

Также отметим, что условие продольного силового равновесия электронов с начала 2000-х годов успешно использовалось в численно-аналитических моделях стационарного пространственно одномерного тонкого токового слоя в хвосте магнитосферы.

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

Цели и задачи данной работы

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

Вторая цель состоит в выводе системы уравнений для бесстолкновительной плазмы в указанном приближении в результате замыкания системы уравнений для полей соответствующей системой уравнений переноса для плазмы в следующих наиболее актуальных случаях:

1) кинетическое описание незамагниченной плазмы, когда системой уравнений переноса является система уравнений Власова для каждой компоненты;

2) кинетическое описание в случае замагниченных электронов и незамагничен-ных ионов, когда система уравнений переноса состоит из уравнений Власова для ионных компонент и уравнения Власова в дрейфовом приближении для замагни-ченных электронов;

3) кинетическое описание в случае полностью замагниченной плазмы, когда си-

и и и 1 > и 1

стемой уравнений переноса является система уравнений Власова в дрейфовом приближении для каждой компоненты;

4) кинетическое описание в случае замагниченных электронов и частично замаг-ниченных ионных компонент, когда система уравнений переноса состоит из уравнений Власова для незамагниченных ионных компонент и уравнений Власова в дрейфовом приближении для замагниченных ионных компонент и электронов;

5) газодинамическое описание в случае полностью замагниченной плазмы;

6) гибридное описание плазмы в случае замагниченных электронов и незамагниченных ионов.

ГП «-> «->

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

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

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

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

Для достижения поставленных целей решались следующие задачи.

1. Вывести систему уравнений для магнитного и электрического полей из уравнений Максвелла, условия квазинейтральности и условия силового равновесия электронов вдоль магнитного поля, а также вывести схему её замыкания системой уравнений переноса плазмы.

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

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

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

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

6 1 > и и и и

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

7. Вывести систему уравнений гибридного описания бесстолкновительной плазмы из незамагниченных ионов и замагниченных электронов, с учетом продольного силового равновесия электронов, а также вывести схему численного интегрирования по времени этой системы уравнений.

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

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

10. При помощи этой теории построить численную модель указанного тонкого токового слоя, в которой учитываются электростатические эффекты, электроны имеют распределение Масвелла-Больцмана и их вклад учитывается аналитически, а для протонов численно решается уравнение Власова.

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

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

Методы исследования

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

моделирование тонких токовых слоев, а также анализ результатов численного моделирования и сравнение этих результатов с данными измерений. Используемая в диссертации модель реализована в виде комплексов программ на языке Fortran.

Научная новизна

1. Показано, что система уравнений Власова-Максвелла в крупномасштабных задачах для медленных плазменных процессов в пространственно неоднородной плазме не позволяет моделировать электростатические эффекты, вызванные относительно очень малым отклонением от электронейтральности.

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

3. Впервые получены замкнутые системы уравнений для описания безызлу-чательных процессов в бесстолкновительной плазме с учетом продольного силового равновесия электронов для 6-ти наиболее актуальных для физики космической плазмы вариантов системы уравнений переноса:

1) в случае, когда все компоненты плазмы не замагничены и описываются уравнениями Власова;

2) в случае, когда электроны замагничены и описываются уравнением Власова в дрейфовом приближении, а все ионы не замагничены и описываются уравнениями Власова;

3) в случае, когда все компоненты плазмы замагничены и описываются уравнениями Власова в дрейфовом приближении;

4) в случае замагниченных электронов и частично замагниченных ионов, когда электроны и замагниченные ионные компоненты описываются уравнениями Власова в дрейфовом приближении, а незамагниченные ионные компоненты описываются уравнениями Власова;

5) в случае, когда все компоненты плазмы замагничены и описываются системой уравнений многокомпонентной магнитной газовой динамики для за-магниченной плазмы;

6) в случае гибридного описания плазмы, когда ионы описываются уравнениями Власова, а замагниченные электроны описываются в рамках газовой

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

4. Впервые для бесстолкновительной плазмы из протонов и замагниченных электронов получена новая форма уравнения силового баланса протонов, в которой исключено электрическое поле. Это уравнение необходимо в моделях токовых слоев для правильной постановки граничных условий и для контроля силового баланса. Из этого уравнения получена векторная форма условия силового баланса в виде закона сохранения для стационарного пространственно одномерного токового слоя с учетом анизотропии давления электронов.

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

Список литературы диссертационного исследования доктор наук Мингалев Олег Викторович, 2020 год

- -

- B(z) -

- в (2) У —

— ------- о —

_ ---Biz) _

: :

- —

-9 -10 -11 -12 -13 -14 -15 -16 -17 -18 -19 -20 -21

-1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Е

Рис. У.17. Профили компонент магнитного поля Вх(г), Ву(г) и его величины В (г) в конфигурации ТТС с "колоколообразным" Ву (г) при Вг = 2 нТл

и 8 = Ув/ УТр 0 =2 .

Для таких же значений входных параметров Вг = 2 нТл и 8 = У^ Утр о = 2 , как и в представленной на рисунках У.15 и У.16 симметричной конфигурации ТТС, на рисунке У.17 представлены профили самосогласованных компонент Вх(г) и Ву(г) магнитного поля в нТл, на рисунке У.18 представлены профили компоненты плотности тока 3у(г) = з (г) и 3х(г) = 3рх(г) в нА/м , а на рисунке У.19 представлен профиль концентрации п(г) в ем-3 .

Сравнение рисунка У.17 с рисунками У.1 и У.15, а также рисунка У.18 с рисунками У.2 и У.16 показывает, что в полученной конфигурации при одинаковом полном токе через слой

Г Яв

Зу = 3у) ^ ,

то есть при одинаковом изменении при переходе через ТС тангенциальной компоненты магнитного поля АЕХ = Ех(Ь^ — Ех(—Ь^ , токовый слой в направлении оси У стал примерно в 3 раза шире, а максимальное значение компоненты плотности тока ]у(г) соответственно примерно в 3 раза уменьшилось.

./„00, (нА/м2)

24 22 20 18 16 14 12 10

-2 -4 -6 -8 -10 -12 -14

.......

----/%)

■ , Ч . (5) , . (Щ ■ ] (г) =] +7

у^ у у

-0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

/(г), (нА/м2)

•'х

■ X X -

]х

1

20 18 16 14 12 10

-2

-10

-12 -14 -16

-18 -20

-0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

г/.Я

Е

Рис. У.18. Профили плотности тока уу(г) — верхняя панель, и уХ(г) — нижняя панель, в конфигурации ТТС с "колоколообразным" Еу(г) при Ег = 2 нТл и

5 = У^ УТр 0 = 2 . Плотности тока от южного источника у (г) и уХ^ (г) обо-

• N)/ ч )/ ч

значены синеи линиеи, плотности тока от северного источника у у (г) и уХ (г) обозначены черноИ линиеИ, полные плотности тока у у (г) = у ^ (г) + у ) (г) и 3х(г) = 3Х^(г) + 3ХМ)(г) обозначены красноИ линиеИ.

Из рисунка У.18 видно, что попадающие в слой сверху протоны от северного источника дают отрицательный вклад у) (г) < 0 в плотность тока ]у(г), кото-

рый обозначен черной линией, а попадающие в слой снизу протоны от южного источника дают положительный вклад 3 ^ (г) > 0 , который обозначен синей линией. При этом максимальное значение тока от южного источника примерно в 1.5 раза больше модуля минимального значения тока от северного источника.

"(г)

о.з

0.29 0.28 0.27 0.26 0.25 0.24 0.23 0.22 0.21 0.2 0.19 0.18 0.17 0.16 0.15 0.14 0.13 0.12 0.11 0.1 0.09 0.08 -1

( см"3 )

- -

(S

Xz) -

- (N гг

(s) + (N) n J

-n(z ) = n -

- -

-

- -

- -

- -

-0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

z/R„

Рис. V.19. Профиль концентрации n(z) в конфигурации ТТС с "колоколообраз-ным" By (z) при Bz = 2 нТл, S = VD/ VTp 0 = 2 . Концентрация от южного источника n(S) (z) обозначена синей линией, концентрация от северного источника n(N) (z) обозначена черной линией, полная концентрация n(z) = n (S ) (z ) + n(N ) (z ) обозначена красной линией.

Отслеживание и анализ рассчитанных фазовых траекторий протонов показал, что фазовые траектории основной части тепловых частиц, на которых сосредоточены функции распределения в обоих источниках, и которые поддерживают текущий в слое ток, достаточно быстро пересекают слой либо без осцилляций в слое, либо с небольшим числом таких осцилляций. То есть в смысле терминологии, введенной в работах [Zeleniy et al., 2000; Sitnov et al., 2000] фазовые траектории основной части тепловых частиц являются пролетными. Эта картина сильно отличается от таковой для рассмотренных выше симметричных слоев, в которых фазовые траектории значительной части переносящих ток тепловых частиц являются квазизахваченными, причем существенная их доля "отражается" от слоя, то есть после осцилляций в слое вылетает из слоя обратно в сторону источника.

1) Баланс по X (нПа): lljz) - В BJ,z)///()~0, 2) по Y : lljz) - В Ву(г) / /«0 ~ О

г J

0.03 0.025 0.02 0.015 0.01 0.005 О

-0.005 -0.01 -0.015 -0.02 -0.025 -0.03

/

/ * /

t

* /

/ »

/ ! 1 T

/ * t

!

/ * / ■ П (z)

/ _ ■ В В (z) /«

/ „X / ■ Я (г) - XZK 7 В В (z) z 7 г 0

0.04 0.035 0.03 0.025 0.02 0.015 0.01 0.005 О

— п (z) yzv

в в (z) lu z у^ ' г 0 п (:)-в в (z) /« yzK z у 0

"V

s

/ / \ к.

/ \ \

/ / / \ \

\ \ \ \

s / /

/ ✓ \ у V

-0.9-0.8-0.7-0.6-0.5-0.4-0.3-0.2-0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 -0.9-0.8-0.7-0.6-0.5-0.4-0.3-0.2-0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

3) Баланс по оси Z (нПа): Я (г) + ( (г)|2+ (г)|2 ) / (2цЛ ) +po(z) « Const

0.38

1 1 1 1 1 1

00 —

.... я Z2

-------(|В<_~)|2 + IВ (z) 2)/(2А0) + \By{zf)!(7

(г)-К \Bx(zf А0) -

ZZ

-pew

-----я (z) + (\Bx(z)\2 + \By(zf)l(2„Q)+pe(z) -

Z2

-1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

(*=2нТл, ¿=(VD/VTpo) = 2) z/Re

Рис. V.20. Силовой баланс в конфигурации ТТС с "колоколообразным" By(z) при Bz = 2 нТл = VD/VTp0 = 2 .

1) Верхняя панель — баланс по оси X : Пх z (z) — фиолетовая линия, BzBx(z)/ ц0 — синяя линия, полный баланс Пх z(z) — BzBx(z)/ д0 ~ 0 — красная линия.

2) Средняя панель — баланс по оси Y : Пу z (z) — фиолетовая линия, BzBy (z) j

синяя линия, полный баланс Пу z (z) — BzBy (z) j ¡л0 ~ 0 — красная линия.

3) Нижняя панель — баланс по оси Z : Пг z (z) — фиолетовая линия,

(^Bx(z) + B^(z{2^0) — синяя линия, pe)z) — коричневая линия, баланс без учета электронов Пг z(z) + ^B^(z) + B^z^(2^0) — зелёная линия,

, z(z) + (bX(z) + By^(z) {2^0) + pe(z) ~ Const — красная ли-

полный баланс П2 z(z

ния.

На рисунке V.20 показан силовой баланс в рассматриваемой конфигурации ТТС. Верхняя левая панель показывает баланс по оси X , который определяется уравнением (V.14). Как видно из рисунка, силовой баланс выполнен с очень хорошей точностью, причем постоянная в правой части формулы (V.14) равна нулю, то есть это уравнение принимает вид

Пх,2(z) - BzBx(z)/^ = 0 . Верхняя правая панель показывает баланс по оси Y , который определяется уравнением (V.15) и также выполнен с очень хорошей точностью, причем постоянная в правой части формулы (V.15) также равна нулю, то есть это уравнение принимает аналогичный вид /

Пу z z (z) - Bz By (z)/^ = 0 .

На нижней панели показан баланс по оси Z, который определяется уравнением (V.16). Это уравнение с учетом того, что магнитное поле имеет вид (V.2), можно представить в форме

nz , z (z) + (Bx2(z) + By2(z)) + pe (z) « Const .

Красная линия показывает выполнение этого уравнения с очень высокой точностью.

Отметим, что из рисунка V.17 видно, что величина магнитного поля в слое относительно его значений на краях немного повышается. Это значит, что также ведет себя и магнитное давление B2(z)j{j^o} , которое показано синей линией. Также отметим, что из рисунка V.19 видно, что концентрация n(z) в слое относительно значений на краях немного понижается. Это вызывает аналогичное понижение компоненты тензора напряжений протонов Пг z (z), которая отображается фиолетовой линией, а также аналогичное понижение показанного коричневой линией давления электронов pe(z), так как их температура постоянна.

При этом, как видно из нижней панели рисунка, давление электронов близко к постоянному и равно примерно « 0.02 нПа, что дает около 1/18 части от полного напряжения по оси Z , которое обозначено красной линией и примерно равно « 0.36 нПа.

Отметим, что из формулы (V.18) и близости давления электронов к постоянному вытекает, что в рассматриваемой конфигурации ТТС с "колоколообразным" By(z) электрическое поле близко к нулю во всем слое, потенциал электрического поля, определяемый первой формулой в (V.36), близок к постоянному, то есть электростатические эффекты очень малы.

1) Ppi(z), Pp2(z), Pp3(z), ( нПа )

0.54 0.52 0.5 0.48 0.46 0.44 0.42 0.4 0.38 0.36 0.34 0.32 0.3 0.28 0.26 0.24 0.22 0.2 0.18 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0

-1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1

2) Г ¿z), Тр 2(z) , Тр 3(z), (КэВ)

1 1 [ 1 1

--

- -

- _ _ PP2 fz) -

— (Z) —

_

— (z)

: yj

. ___ - —

_ _

— —

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

0.9

11.5 11

10.5 10

9.5 9 8.5 8 7.5 7 6.5 6 5.5 5

4.5 4 3.5 3 2.5 2 1.5 1

0.5 О

Vz)

J_

J_

J_

J_

_L

J_

J_

_L

_L

J_

-0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1

0.1 0.2

0.3

z/R

0.4 0.5 0.6 0.7 0.8

0.9

Рис. V.21. 1) Верхняя панель — собственные числа

Ppi(z) < Pp2(z) <Pp3(z)= ppll(z)

тензора давления протонов PPp(z) в нПа ,

2) нижняя панель — собственные числа Tpl(z) < Tp2(z) < Tp3(z) = Tp^(z) тензора

температуры протонов rTp(z) = ~Pp(z)j(enp(z)) в КэВ в конфигурации ТТС с "колоколообразным" By(z) при Bz = 2 нТл и S = VD/ VTp 0 = 2 . Наименьшие собственные числа Ppi(z) и Tpi(z) обозначены синей линией, вторые собственные числа P 2(z) и Tp2(z) обозначены красной линией, и максимальные собственные числа Pp3(z)= ppn(z) и Tp3(z) = Tp^(z) обозначены фиолетовой линией.

На рисунке V.21 для рассматриваемой конфигурации ТТС на верхней панели показаны в нПа, занумерованные в порядке возрастания собственные числа

Ppl(z) < Pp2(z) <Pp3(z)= Ppn(z)

тензора давления протонов Рр(г), который определяется формулами (1.9). На нижней панели показаны аналогично занумерованные еобственные числа

Тр1{г) < Тр2{г) < Тр3{г) = Т^

тензора температуры протонов Тр(г) = Рр(г) ^ (впр(г)) в КэВ , который определяется формулами (1.10).

При этом максимальные собственные числа равны, соответственно, продольному давлению и температуре, то есть

Рр3(г) = Рр \\(г) , Тр3(г) = Тр\\(г) ,

и магнитное поле является собственным вектором для этих максимальных собственных чисел. Из рисунков также видно слабое нарушение гиротропности рассматриваемых тензоров в токовом слое при | г/ ЯЕ | < 0.5 , поскольку в ортогональной магнитному полю плоскости имеется два собственных направления с хоть и близкими, но различными собственными числами

Рр1(г) < Рр2(г) , Тр1(г) < Тр2(х) .

Отсюда следует, что протоны в этой конфигурации не являются полностью за-магниченными.

/'' (-)- ' <-> (пТ)

1 1 -

г -

- ^VS '//У Ч ч -

- .Л У Члч •/// Ш/ -

- J' V/ -

- V У'/ -

- 'Л* // / '/ / "/ ' 4\vN /t. i *Ч -

- s/ & / \\\ -

- S * / % -

- // -

- 1 ni -

- / / / ---------li (_ j j . is \ -

- -------« (--)). a 1 lll -

- / .......... H (ri). 2 ni 2 ni -

- / / /, \ -

- d / ----------«.(--))• « -

- у7■ ----------вш В 3 ni 3 ni -

- y4 D D -

- Î* -------------- tS^Z)) , ts -

- ^ £ -By(z)), £ 4 il 1 -

,Î s 2 = --------------В (z)), В 4 ni -

-

20 18 16 14 12 10 8 6 4 2 0 -2 -4 -6 -8 -10 -12 -14 -16 -18 -20

-1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Z/Re

Рис. V.22. Профили компонент магнитного поля Bx(z) и By(z) в конфигурациях ТТС с "колоколообразным" профилем By (z) при Bz = 1, 2 , 3 , 4 нТл и S = Vd/ Vtp о = 2 •

Отметим, что аналогичные стационарные конфигурации ТТС с "колоколооб-разным" профилем By(z) были получена также и для других значений нормальной компоненты магнитного поля Bz = 1, 3 , 4 нТл . Эти конфигурации оказались полностью аналогичны рассмотренной конфигурации для Bz = 2 нТл.

На рисунке V.22 приведены профили самосогласованных компонент магнитного поля Bx(z) и By(z) в указанных конфигурациях. Из рисунка видно, что профили самосогласованных компонент магнитного поля относительно слабо зависят от значения нормальной компоненты магнитного поля Bz .

V.8 Заключение главы V

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

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

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

области моделирования рассчитываются не для всех узлов сетки в пространстве скоростей, а только для "засвеченных" на первом шаге узлов.

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

Впервые в численном моделировании получены симметричные конфигурации стационарного пространственно одномерного тонкого токового слоя с заданной нормальной компонентой магнитного поля, в которых очень точно выполнены условия силового баланса. При этом впервые рассчитана функция распределения протонов и детально исследованы электростатические эффекты с кинетическим описанием замагниченных электронов, а также исследовано влияние на конфигурацию слоя анизотропии давления электронов и гидродинамической скорости образующих токовый слой потоков протонов. Особо стоит подчеркнуть, что в этих конфигурациях из-за наличия и важной роли захваченных электронов получить внутри ТС правильную электронную концентрацию при помощи только численного расчета фазовых траекторий электронов, влетающих в ТС от границы области моделирования, в принципе невозможно.

Впервые в численном моделировании получены и детально исследованы конфигурации тонкого токового слоя с "колоколообразным" профилем сдвиговой компоненты магнитного поля и примерно постоянными величиной магнитного поля и концентрацией. Конфигурации такого типа часто встречаются в данных измерений на космических аппаратах в хвосте магнитосферы Земли и на разных участках ее магнитопаузы, а также в солнечном ветре и в хвосте магнитосферы Юпитера.

Результаты моделирования качественно и количественно очень хорошо согласуются с данными измерений на космических аппаратах, в которых регулярно встречаются токовые слои указанных выше типов. При этом конфигурации с "колоколообразным" профилем сдвиговой компоненты были впервые получены при помощи численного моделирования в работах [Мингалев и др., 2012; Ма1оуа,

Mingalev et al, 2015].

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

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

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

Второй случай относится к токовому слою в хвосте магнитосферы Юпитера. Одно из важных отличий магнитосферы Юпитера от магнитосферы Земли состоит в источниках плазмы: если магнитосфера Земли, по большей части, заполнена протонами солнечного ветра, то в магнитосфере Юпитера существенны вклады популяций тяжёлых ионов, происхождение которых связывают с вулканической деятельностью спутников Юпитера.

Так, в магнитосфере Земли в спокойных геомагнитных условиях в токовом слое хвоста имеется в основном один сорт ионов — горячие протоны с температурой примерно 4-10 КэВ, в то время как в магнитосфере Юпитера помимо немного более горячих протонов с температурой примерно 5-20 КэВ, имеются в разы более горячие ионы кислорода с температурой примерно 20-40 КэВ и еще более горячие ионы серы с температурой примерно 40-60 КэВ.

Л и и 1 О и и

Отметим, что токовый слой хвоста магнитосферы Земли и токовый слой хвоста магнитосферы Юпитера имеют близкие значения магнитного поля, но толщина ТС в хвосте магнитосферы Юпитера больше толщины слоя в хвосте магнитосферы Земли в число раз, примерно равное отношению их радиусов, которое равно Rj/ Re ~ 11.2 , то есть на порядок.

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

различных вариантов параметров плазмы, а затем сравнить полученные модельные конфигурации с данными измерений космического аппарата Juno.

По данным космического аппарата Juno в хвосте магнитосферы Юпитера наиболее часто наблюдаются ТС с "колоколообразным" профилем сдвиговой компоненты, с полутолщиной примерно 50000 км и величиной магнитного поля примерно 20 нТл.

В таком ТС протоны с температурой 5-20 КэВ будут иметь тепловой гирора-диус в 350-700 км и могут рассматриваться как замагниченные, то есть их вклад в поля и плотность тока может учитываться аналитически в рамках дрейфовой теории аналогично вкладу электронов.

Приложение. Новые методы численного решения стационарного уравнения Власова

Здесь будут изложены два связанных между собой новых метода численного решения стационарного уравнения Власова. Также будет описан разработанный автором и показавший очень высокую эффективность алгоритм численного решения задачи Коши для системы уравнений Ньютона-Лоренца, который имеет 2-й порядок точности, допускает переменный шаг по времени, а также явный и неявный варианты, и основан на точном решении этой задачи в случае постоянных полей. Этот новый прием аналитически учитывает вращение в магнитном поле, дает нулевую фазовую ошибку и снимает необходимость выполнять большое число шагов по времени на период циклотронного вращения частицы, то есть в случае плавно изменяющихся полей позволяет на несколько порядков увеличить временной шаг по сравнению с традиционными конечноразностными алгоритмами. Это позволяет экономично и с высокой точностью рассчитывать траектории как ионов, так и электронов.

1. Обоснование численных методов

Рассмотрим нерелятивистскую систему стационарных уравнений Власова, которая описывает плазму из К сортов ионов (нижний индекс а = 1,...,К) и электронов (нижний индекс а = е) в ограниченной области П пространства

К3 :

(V • + — ((Е(х) + [ V х В(х)]) • ^ = 0 , а = 1,...,К,е. (VI. 1)

\ Ох / та\ ду )

Допустим, имеется классическое решение краевой задачи для этой системы уравнений, и поля Е(х) и В(х) известны также и вне области П . Тогда через каждую точку

Z = (Х, V) е П х К3

(где через П обозначено замыкание области П) фазового пространства одной частицы проходит характеристика каждого уравнения Власова (VI.!), которая является фазовой траекторией автономной системы обыкновенных дифференциальных уравнений Ньютона-Лоренца. В системе единиц СИ в нерелятивистском случае стационарных полей система Ньютона-Лоренца вместе с начальными условиями в момент времени Ь = 0 имеет вид:

^ = у(Ь) , = — (Е(х(Ь)) + [ у(Ь) х В(х(Ь)) ] ) , ,

¿Ь ¿Ь та\ к у к> к к "Ч1 ) (VI.2)

х(0) = X , у(0) = V .

Решение задачи Коши (VI.2) для последующего изложения удобно обозначить как

x(t) = Ra(t, X, V) , v(t) = Ua(t, X, V) , (VI.3)

то есть эти векторные функции удовлетворяют следующим уравнениям и начальным условиям:

дRa(t, X, V) ( )

;; 1 = иЛ t,х, v),

д Uf X' V) = (E (Ra(t, X, V) + [ Ua(t, X, V) X В R(t, X, V) ] ) ,

Ra(0, X, V) = X , Ua(0, X, V) = V .

Как известно, классическое решение уравнения Власова (VI.1) — функция распределения fa(x, v) постоянна вдоль фазовой траектории (VI.3), то есть верно тождество

fa(Ra{t, X, V), Ua(t, X, V)) = fa(X, V) = Const . (VIA)

Обозначим через n(x) единичную внешнюю нормаль к границе дО. Граничные условия для функций распределения состоят в том, что для каждого сорта частиц на множестве

Г(+) = { (x, v) : x £ дО , (v • n(x)) < 0 } , (VI.5)

которое логично назвать областью влета на границе фазового пространства задана функция влета f^^ (x, v). Отметим, что функция распределения вылетающих из области О частиц задана на множестве

Г(- = { (x, v) : x £ дО , (v • n(x)) > 0 } .

То есть множества Г(+ и Г(-) являются разложением границы фазового пространства: дО x R = Г(+) U Г(-). Обозначим через

F(+) = { (Ra(t, x, v), Ua(t, x, v)): (x, v) £ Г(+) ,t > о}р|О x R (VI.6)

множество точек фазового пространства системы (VI.2), через которые проходят фазовые траектории, выпущенные из области влета на границе фазового пространства Г(+ . На этом множестве значения функции распределения fa(x, v) определяются значениями функции влета (x, v) в соответствии с формулой (VI.4):

fa(Ra(t, x, v), Ua(t, x, v) )= fj+)(x, v) , (x, v) £ Г(+) , t > 0 . (VI.7)

Остальную часть фазового пространства системы ^1.2), которую обозначим через

F(-) = ( П х К3 )\F(+), (VI.8)

составляют "внутренние захваченные" фазовые траектории, у которых пространственная координата Яа(Ь, х, V) не покидает множество П (замыкание области П) . Отметим, что "внутренних захваченных" фазовых траекторий может не оказаться совсем. В этом случае множество F(_) является пустым, и на значения функции распределения /а(х, V) на всем фазовом пространстве П х К3 определяются значениями функции влета (х, V) по формуле (VI.?). Если множество не пусто, то в ходе моделирования на нем значения функции распределения /а(х, V) должны задаваться согласно специфике конкретной задачи.

Идея первого метода состоит в том, чтобы в уже известных полях текущей итерации из узлов 2а[д] некоторой достаточно подробной сетки (где q —

вектор с целыми координатами, нумерующий эти узлы) на множестве влета Г(+ , которая "покрывает" носитель функции влета / (х, V), "выпустить вперед" при Ь > 0 фазовые траектории

(Ь, ^]), иа(Ь, Za )) , Ь > 0, Za ^] е С 4+) (V1.9)

в область П . По этим фазовым траекториям в узлах г[кх] некоторой регуляр-

и и и М

ной пространственной сетки при помощи специальной процедуры "взвешивания" рассчитываются используемые в модели моменты функции распределения /а(х, V), а также значения функции распределения / а(г [кх], Va [кх, в уз-

лах Va[кх, kv] регулярной сетки в пространстве скоростей для узла г[кх], которая ориентирована по магнитному полю в этом узле. Здесь кх и kv — векторы с целыми координатами, нумерующие узлы сетки.

При этом шаг регулярной сетки в пространстве скоростей частиц сорта а Ауа следует определять через характерную температуру этих частиц Та0 (в электрон-вольтах) по формуле Ауа = ^Ута0 = Ута0 / Nу0 , где Nу0 е N — натуральное число, 7 = 1/Nv0 ~ 0.1 и Ута0 = \!еТа0/та — их характерная тепловая скорость.

Далее по найденным моментам в результате численного решения краевых задач для нахождения стационарных полей рассчитываются поля следующей итерации.

Идея второго метода состоит в том, чтобы в уже известных полях текущей итерации из узлов {г[кх], Vа[кх, kv]) указанной выше сетки в фазовом пространстве "выпустить назад в прошлое" при Ь < 0 фазовые траектории

(Яа(Ь, г[кх], V[кх, К]), иа(Ь, г[кх], V[кх, К])) , Ь < 0 . (^.10)

Заметим, что для любой точки (X, V) е (П х К3) У Г(-) возможны только два случая:

1) при некотором конечном времени £ = (X, V) < 0 фазовая траектория (У1.10) (В,а(Ь, X, V), иа(Ь, X, V)) попадает в область влета Г(+), на которой задана функция влета (х, V), а при -£(+) (X, V) <£ < 0 фазовая траектория не выходит из области П , то есть выполнены следующие условия:

яа(г, X, V) е п при - , V) <г< о , яа(-£{+)(X, V), X, V) е дп, (п(х) • иа(г, X, V) ) < о .

В этом случае фазовая траектория лежит в подмножестве фазового про-

странства.

2) при любом конечном времени £ фазовая траектория (У1.10) не выходит из П , то есть На(г, X, V) е П и £(+ (X, V) = . В этом случае фазовая траектория лежит в подмножестве фазового пространства

На фазовой траектории 1-го типа из подмножества F(+) в силу тождества (У1.4) функция распределения равна функции влета в точке

'пЛ-£(+)X, V) X, V), иЛ-Ч+)X, V), X, V)) е г(+),

лежащей в области влета на границе фазового пространства:

Ш,V) = й+)(ПаН^,V),X, V),иЛ-Ч+)X,V\X, V)) . (п.11)

Как уже отмечалось выше, для фазовой траектории 2-го типа из подмножества F(-)

нет общих явных условий, которые бы определяли на ней функцию распределения. В процессе моделирования либо функция распределения на таких фазовых траекториях, либо вклад ее значений на этих траекториях в ее моменты должны определяться из специфики конкретной задачи. В некоторых случаях может быть явное описание таких захваченных траекторий, примером которого является условие (1У.30) для системы дрейфовой системы уравнений (1У.22) в плоском токовом слое.

Также отметим, что в ходе численного моделирования фазовые траектории вычисляются не точно, а с некоторой погрешностью. Кроме того, стационарное решение является лишь абстрактной идеализацией, с помощью которой аппроксимируют достаточно медленно изменяющееся состояние реальной системы. Поэтому у реальной системы всегда имеется характерное время Ос существования квазистационарного решения. Отсюда следует, что в ходе численного решения достаточно рассчитывать фазовые траектории (У1.10) в прошлое на время контроля Ок ~ 100с . Если за это время фазовая траектория не выходит из области П, то она считается траекторией 2-го типа.

2. Дискретизация в координатном пространстве и в фазовом пространстве

Отметим, что детали дискретизации сильно зависят от специфики конкретной задачи. Поэтому здесь мы рассмотрим только общую методику. Обозначим через {е1, е2 , е3} декартов базис в пространстве К3 , то есть

3 3

x = xi ei , v = J2 vi ei i=1 i=l

Будем считать, что область моделирования является прямоугольником с центром

3

в точке r0 r0i ei :

i=l

п = {x : Xi - roi I <LXi = NXiax, i = 1, 2,3 } ,

и и и и Л

который связан с регулярной пространственной сеткой с шагом ах и узлами

3

r[kx]= ro + кх Ах = 5>i + кхг Ах) ei , кх = (к^ ,kx2 ^3) е Z3 . (VI.12)

i=l

Отметим, что в задачах для космической плазмы обычно присутствует существенное всюду ненулевое магнитное поле. Для этого случая в пространстве скоростей удобно использовать сетку, которая аналогично сетке (V.52) в случае токового слоя связана с местным магнитным полем B(r[kx^ следующим образом.

Для каждого узла регулярной пространственной сетки x = r [kx] в пространстве скоростей используется локальная ортонормированная система координат, ориентированная по магнитному полю, у которой 3-й базисный вектор h33(x) направлен вдоль магнитного поля, а 1-й и 2-й базисные векторы hi(x) и h2(x) лежат в плоскости, ортогональной магнитному полю и образуют правую тройку векторов, то есть:

h3(x) = b(x), h1(x)= [h2(x) x b(x) ] , h2(x) = [b(x) x h1(x) ] . (VI.13)

Компоненты вектора скорости в этой системе координат обозначим через

Fi (x) = (v • hi (x) ) ,

то есть 3-я компонента является продольной скоростью: V:3(x) = (v • b(x)) = Vy , а вектор скорости представляется в виде

33

v = £ Vi (x) hi (x) = Y, Vi i=l i=l

Для данного узла пространственной сетки г[кх] узлы сетки в пространстве скоростей частиц сорта а, ориентированной по магнитному полю, нумеруются вектором

= (^1 ,kv2 ,kvз) е с целыми координатами и определяется формулами

^ [К , К ] =Ava (kvl - hi (r [kx]) , (VI. 14)

которые определяют "полуцелые узлы". Размер сетки по каждому измерению зависит от задачи и определяется из условия, что эта сетка должна "с запасом" содержать носитель функции распределения fa(r[kx], v) по скоростям.

Дискретизация в координатном пространстве состоит в замене полей и моментов функции функций распределения конечным набором их значений в узлах r [kx] пространственной сетки :

F(x) ^ Fh[kx] = F(r[kx]) , F = B, E, na , ja , ua , Пa . (VI.15)

Значение полей в произвольной точке x G О, как и в методе частиц, вычисляются по значениям в ближайших к точке x узлах r [kx^ пространственной сетки при помощи линейного взвешивания:

F(x) = Е ^Fh [kx] , (VI16)

где сумма берется по соседним с точкой x узлам пространственной сетки, и безразмерная линейная передаточная функция взвешивания 1-го порядка (по 2-м точкам на измерение) в пространстве R3 определяется формулами

¿3 (t! ,t2 ,t3 ) = W (t!) • W (t2) • W (t3), где W (t) = max{ 1 - \t\;0 } . (VI.17)

Функция распределения fa(x, v) заменяется конечным набором своих значений в узлах сетки в фазовом пространстве (r [kx], Va [kx , kv]) :

fax v) ^ fah [kx , kv ] = fa (r [kx] , Va[kx , kv ] ) . (VI.18)

При этом интегралы по пространству скоростей в формулах (1.2), (1.3) и (1.9), которые определяют концентрацию, плотность тока и тензор напряжений данной

компоненты плазмы заменяются соответствующими суммами:

[К] ^ [К , К] ,

к

3а1г [К] = )3 £ Уа [кХ , К ] а [кх , К ] ) ,

к

^ [К] = та Аа )' £ Уа [кХ , К ] ® Уа [К , К ] ¡а1г [К , К ]

(У1.19)

где сумма берется по всем kv , то есть по всем узлам сетки в пространстве скоростей для данного узла пространственной сетки.

3. Схема итерационного процесса 1-го метода

Рассмотрим кратко основные детали схемы итерационного процесса в терминах перехода к ^ к +1, к = 0,1,... от текущего приближения, обозначенного верхним индексом к, к следующему, которое обозначено верхним индексом к +1, для полей и моментов функций распределения. При этом в минимальной постановке, если рассчитываются только поля и гидродинамические параметры плазмы, сами функции распределения можно и не рассчитывать, поскольку в ходе моделирования это требует большого расхода оперативной памяти.

Пусть известно текущее приближение для полей В^ [кХ] и Е^ [кХ] на пространственной сетке Пh с узлами (У1.12), которая содержит кроме прямоугольника ^ еще несколько слоев.

Пусть — достаточно подробная сетка узлов ^г [кХ], У а [кХ , к^ ] ^ в определяемой формулой (У1.5) области влета Г(+ фазового пространства системы, которая покрывает носитель функции влета (х, V).

Из узлов этой сетки "выпускаем вперед" при Ь > 0 фазовые траектории (У1.3) системы (У1.2):

х(Ь) = Яа(1, г [кХх] , Уа[кХ , ]) , v(t) = Ц^Ь, Г [кХ] , У^Х , к? ]) , (У1.20)

и рассчитываем их вплоть до вылета из содержащей область ^ большей области контроля ^ . Значение функции распределения на этой траектории в соответствии с формулой (У1.4) равно значению функции влета в исходной точке:

/а (х(Ь), ^;)) = /<+> [кХ , к1 ] = /<+> (т[к°х], Уа [кХ, к1 ]) . (У1.21)

При этом в ходе расчета фазовой траектории "отслеживается" все ее пересечения в пространстве с границами ячеек сетки (У1.12), когда хотя бы одна из пространственных координат х,1 (Ь) фазовой траектории равна значению такой же координаты ^ [кх4\ = г? + кХ{ Ах для узла пространственной сетки (У1.12), то есть

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

t : 3 i g{1, 2, 3} : я.(t) = r. [kxi] = r0i + kxi Ax . (VI.22)

В эти моменты выполняется "взвешивание" фазовой траектории, то есть расчет ее вкладов для ближайших узлов пространственной сетки в сеточные массивы моментов функции распределения и в сеточный массив функции распределения, если она также рассчитывается. При этом в пространственно одномерном случае x = x G R траектория будет точно попадать в узлы сетки r [kx] = r0 + kx Ax , как это отмечалось в главе V. Поэтому вклады в сеточные массивы моментов функции распределения будут накапливаться "напрямую" в этих узлах без взвешивания по формулам

Snah [kx ]= f(h } [kX , К]

0

(VI.23)

ja h [kx ] = q a v(t) f (Ж , k°v] ,

^Паh [kx]= ma V(t) 0 V(t) f(+)[k°x , ^ ] . Вклад в сеточный массив функции распределения в узле пространственной сетки будет вносится в 8 ближайших к v(t) узлов сетки (VI.14) в пространстве скоростей по формулам линейного взвешивания

г [ ] Г+U 0 0 (Uv(t) - Va k , kv ]) • hi №x]))\ , N fa h [kx , k v ] = fh kx , k°v] £ П W -^-J- ) , (VI.24)

kv l=l \ a )

где функция W1 (t) = max{ 1 — \t\ ; 0 } определена в (VI.17), а векторы hl(x) связанного с магнитным полем ортонормированного базиса определены в (VI.13).

В пространственно двумерном случае фазовая траектория в пространстве будет пересекать отрезок, соединяющий два соседних узла сетки, в которых она и будет вносить вклады в сеточные массивы моментов функции распределения и в сеточные массив функции распределения в пространстве скоростей. Вклады в этих узлах будут накапливаться при помощи линейного взвешивания только по одной дополнительной к Xi координате, "не попавшей" в значение для узла сетки по формулам, аналогичным формуле (VI.15). То есть формулы (VI.23) и (VI.24) изменятся соответствующим образом, и для программирования их надо выписать для каждой из двух координат. Ввиду их громоздкости, а также не использования их в этой работе, здесь их приводить не будем.

В пространственно трехмерном случае траектория будет пересекать квадратную грань ячейки пространственной сетки (VI.12), то есть вносить перечисленные выше вклады в 4 узла — в вершины этой грани. То есть формулы (VI.23) и (VI.24) изменятся с учетом линейного взвешивания уже по двум координатам.

Для программирования эти формулы также надо выписать для каждой из трех координат. Ввиду их громоздкости, а также не использования их в этой работе, здесь их также приводить не будем.

В результате описанного трассирования из источников на границе мы найдем на сетке значения необходимых для расчета полей моментов функций распределения каждой компоненты плазмы и плотность полного тока [кХ] в текущей итерации к . Затем в результате решения краевой задачи для системы из уравнения Гаусса (1.4) и уравнения Ампера (1.44) находим в следующей итерации магнитное поле В^^^ [кХ] . Далее из соответствующей специфике задачи системы уравнений находим электрическое поле Е^^ [кХ] в следующей итерации. Для контроля сходимости итерационного процесса вычисляется относительное изменение магнитного поля в ходе текущей итерации

2\\ В(к+1) — В(к) II

4К+1) = -, (™5)