автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Численное моделирование фарлей-бунемановской неустойчивости в ионосфере Земли

МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ имени М. В. Ломоносова

Факультет вычислительной математики и кибернетики

На правах рукописи

Ковалёв Дмитрий Владимирович

2 7 АВГ 2009

Численное моделирование фарлей-бунемановской неустойчивости в ионосфере Земли

Специальность 05.13.18 - математическое моделирование, численные методы и комплексы программ

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

Москва, 2009

003475854

Работа выполнена на кафедре автоматизации научных исследований факультета вычислительной математики и кибернетики Московского государственного университета имени М. В.Ломоносова.

Научный руководитель: кандидат физико-математических наук,

доцент Смирнов Александр Павлович

Официальные оппоненты: доктор физико-математических наук,

профессор Еленин Георгий Георгиевич

Защита состоится 23 сентября 2009 г. в 15 час. 30 мин. на заседании диссертационного совета Д 501.001.43 при Московском государственном университете имени М.В.Ломоносова по адресу: 119991, Москва, Ленинские горы, МГУ им. М. В. Ломоносова, факультет вычислительной математики и кибернетики, ауд. 685.

С диссертацией можно ознакомиться в библиотеке факультета вычислительной математики и кибернетики МГУ им. М.В.Ломоносова.

Автореферат разослан 1В 2009 г.

Ученый секретарь

доктор физико-математических наук, профессор Истомин Яков Николаевич

Ведущая организация: Институт прикладной математики

им. М. В. Келдыша РАН

диссертационного совета, профессор

Общая характеристика работы

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

Ионосферные аномалии интенсивно исследуются с помощью запусков ракет, лабораторных экспериментов и радиопередающих комплексов, таких как Jicamarca, EISCAT, SPEAR, HAARP, Сура. Но, несмотря на полувековую историю экспериментальных и теоретических исследований нелинейных процессов в плазме атмосферы Земли, на данный момент не существует единой теории, описывающей нелинейную динамику плазмы. Усложнение новых теоретических моделей приводит к тому, что влияние различных внешних факторов на нелинейное поведение плазмы часто не может быть оценено без компьютерного моделирования.

На текущий момент в мире численное моделирование фарлей-бунемановской неустойчивости выполняется в Центре космической физики Бостонского университета США с помощью программного кода, основанного на методе частиц. Моделирование проводится на вычислительном комплексе IBM Blue Gene, установленном в Бостонском университете.

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

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

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

Многомерная структура уравнений требует использования новейших высокопроизводительных многопроцессорных вычислительных комплексов. Отдельная глава диссертационной работы посвящена оптимизации программного кода, созданного для моделирования фарлей-бунемановской неустойчивости, для выполнения на современных массивно-параллельных вычислительных системах с распределенной памятью, таких как IBM Blue Gene, установленный на факультете ВМК МГУ имени М.В.Ломоносова, а также СКИФ МГУ "Чебышёв".

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

Цель диссертации.

1. Разработка модели фарлей-бунемановской неустойчивости.

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

3. Оптимизация разработанных алгоритмов и программ для использования на современных многопроцессорных вычислительных комплексах.

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

Научная новизна и практическая ценность.

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

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

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

Апробация работы.

Основные результаты работы докладывались на:

• XXIV генеральной ассамблее Международного союза геодезии и геофизики (Италия, Перуджа, Июль 2007);

• североамериканской встрече Международного союза радионаук, организованной национальными комитетами Канады и США (Канада, Оттава, Июль 2007);

• научном семинаре Центра космической физики Бостонского университета, (США, Бостон, Ноябрь 2007);

• 12 международном симпозиуме по экваториальной аэрономии (Греция, Крит, Май 2008);

• XVI Международной научной конференции студентов, аспирантов и молодых ученых "Ломоносов-2009" (Россия, Москва, Апрель 2009);

• научном семинаре кафедры автоматизации научных исследований факультета вычислительной математики и кибернетики МГУ имени М.В.Ломоносова под руководством зав. кафедры чл.-корр. РАН Д.П.Костомарова (Россия, Москва, 2006-2009).

Публикации. Основные результаты диссертации опубликованы в работах [1-9].

Структура и объем диссертации. Диссертация состоит из введения, трёх глав, заключения, приложения, содержащего рисунки, и списка литературы. Текст изложен на 147 страницах, диссертация содержит 47 рисунков. Список литературы включает 92 наименования.

Содержание работы

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

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

В §1.1 излагается модель неустойчивости. Фарлей-бунемановская неустойчивость - это низкочастотная неустойчивость в плазме, движимая достаточно сильным квазистационарным электрическим полем Е0, перпендикулярным магнитному полю Во'. Эта неустойчивость возникает в слабоионизированной Е-области ионосферы Земли, где электроны замагничены, в то время как ионы из-за частых столкновений с нейтральными частицами практически не подвержены влиянию Во. Ниже Е-области (примерно 95 км) неустойчивость не развивается из-за столкиовительного затухания, а на высотах более 120 км ионы становятся замагниченными. В Е-области из-за частых столкновений незамагниченные ионы сильно подвержены диффузии и не могут разгоняться, поэтому при превышении постоянным электрическим полем порогового значения Еыг — (Ю - 20) мВ/м, что обычно происходит в экваториальных и приполярных электроструях (электроструя - холловский ток в плазме, характерный для Е-области), функция распределения скоростей электронов смещается относительно ионной функции распределения на величину Ц) ~ с|Е0 х Во|/Во- Кроме низкочастотных электростатических колебаний возникают возмущения плотностей ионов и электронов с Щ <С /г±, здесь

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

Условия развития неустойчивости рассматриваются в пункте 1.1.1. Фарлей-бунемановская неустойчивость возникает в верхнем D- и нижнем Е-слое ионосферы Земли на высотах примерно от 80 до 120 км при наличии замагниченных электронов и незамагниченных ионов, что может быть записано в виде следующих условий:

l'en и vm ~ это средние частоты столкновений электронов и ионов с нейтральными частицами, a £îe,t = еВо/те{ - гирочастоты электронов и ионов, массы которых те,< соответственно. В Е-области ионосферы соотношение частот столкновений ионов и электронов с нейтральными частицами близко к постоянному и равно vm с? 10 Vin, а масса ионов приблизительно равна тридцати протонным массам.

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

жидкостную теорию, можно получить уравнение:

где электронные потоки заданы следующим образом

Ге|| = ntV4 = - —i- (V||(Tene) - епе УВФ),

теиеп

re± = neVe± = - ~2 [v± (Тепе) - епе (УхФ - Е0)] +

7Tlelle

GTI

+ neV0 + —£-егх У±Ф,

Triiili

Vex ~ это среднемассовая скорость электронов, Ф - возмущенный электрический потенциал, 0еп = зависимость частоты столкновений от температуры, Те - температура электронов, 7* - температура ионов (она же температура электронов в начале моделирования), Ео - постоянное электрическое поле, перпендикулярное геомагнитному полю Во, е2 = Во/Во - единичный вектор вдоль магнитного поля, Vo = сЕо х ег/Во - дрейфовая скорость, V_l = -f Cyg^ - оператор по координатам, перпендикулярным В0, х и у - координаты вдоль Vo и Ео соответственно, V|| = - оператор по координате z. Здесь и далее температура выражается в Джоулях, при задании температуры в Кельвинах её необходимо домпожать на константу Больцма-на.

Развитие неустойчивости исследуется в прямоугольной области для двумерной модели и в прямоугольном параллелепипеде для трехмерной модели. При решении уравнения для электронной плотности, как и остальных уравнений из модели фарлей-бунемановской неустойчивости, на границе области моделирования используются периодические граничные условия. Начальные условия для плотности электронов пе> как и для плотности ионов щ, которая будет описана далее, задаются в виде периодической функции таким образом, чтобы начальное значение электрического шля УФ было на порядок меньше движущего поля Ео- Это условие выполняется, если разность между начальными значениями плотности электронов и ионов мала.

В пункте 1.1.3 приводится уравнение изменения температуры электронов, которое описывается следующим нелинейным диффузионно-конвекционным уравнением:

п1'г~ = \wenV? - 5ЛЖ ~ Тп), (2)

где ^ = -щ + КхЛ + Уеущ + Уег^ ~ лагранжева производная по времени, = тетп/{те + тп) и тпе - приведенная масса, параметр 6е отвечает за относительное количество энергии, потерянной при одном неупругом соударении (величина 5е на два порядка больше соответствующей величины для упругого соударения 2тпе/тпп), пгп - средняя масса нейтральной частицы (в Е-области примерно равна ионной массе), Уе - среднемассовая скорость электронов.

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

Кинетическое уравнение для ионов описано в пункте 1.1.4, которое с учетом незамагниченности ионов (1) и важного эффекта затухания Ландау, записывается с помощью столкновительного члена вида БГК (Bhatnagax-Gгoss-Krook):

0/ , , т/ , е(Е0 — УФ) д/

_ + (у.7)/ + —--^ =

здесь V - это скорость ионов, /(г, V, £) - функция распределения ионов,

Причем невозмущенная функция ионов является максвелловской /о(г, V, Ь) = щ (з^г) ' ехР (- •

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

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

Наконец, уравнение Пуассона для электрического потенциала приведено в пункте 1.1.5. Оно записывается следующим образом:

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

где Шр,- - плазменная частота, а ф - параметр, отвечающий за высоту в Е-области ионосферы.

В §1.2 приводится описание методов моделирования фарлей-бунемановской неустойчивости. В пункте 1.2.1 рассмотрены методы решения уравнений для электронной температуры и плотности. Основная проблема при решении этих уравнений заключается в том, что при моделировании неустойчивости необходимо, во-первых, проводить моделирование как можно в большей области, во-вторых, как можно лучше аппроксимировать плотности и температуру, то есть увеличивать количество точек на характеристическую единицу длины. Из физических соображений для разрешения дебаевского радиуса на единицу характеристической длины должно приходиться как минимум 5-10 точек, в то время как для моделирования процесса в целом в плоскости (х, у) желательно рассматривать область длиной по каждому направлению как минимум в 20

г0У2Ф = е(пе - п,).

характеристических длин. Многие испробованные разностные схемы для конвекции (Лакс-Вендрофф, Лакс-Фридрихс, Мак-Кормак, различные комбинации этих методов и некоторые другие схемы) не дают достаточной аппроксимации на сетках с большим шагом по пространству (когда на единицу длины приходится 5-10 точек) и оказываются неустойчивыми при больших значениях скоростей переноса. Для проведения расчетов и соблюдения условий развития устойчивости приходится сильно уменьшать шаг по времени, что приводит к увеличению времени выполнения алгоритма.

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

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

Дугласа в трехмерном варианте может быть представлена в следующем виде:

пт*1*~пт = охАхпт+1/3 + (1 - <гх)Ахпт + АуПт + А2пт < "">+г/\-"т+'/3 = ауАу(пт+2^ - пт) = агАг(пт+1 - пт)

пт- значение искомой функции на шаге по времени с номером то, Аху2 -дискретные аналоги операторов дифференцирования, Ы - шаг по времени, - параметры схемы.

В пункте 1.2.2 предлагается решать уравнение Пуассона с помощью быстрого дискретного преобразования Фурье.

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

9/ д/ д/ п

а/ е ЭФ д} е / ЭФ\ дf ,г г.

дЬ тп{ дх дьх пц\ ду) дьу '

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

Схематический алгоритм решения уравнений, описывающих неустойчивость, на одном шаге по времени приводится в пункте 1.2.4:

1. Считая известными плотности электронов и ионов, решить уравнение

Пуассона.

2. Решить ионное уравнение.

3. Решить N раз уравнение для электронов совместно с уравнением Пуассона с шагом по времени ¡ц/'И.

Необходимость более частого вычисления решения электронного уравнения связана с тем, что скорость электронов на порядок больше скорости ионов. И так как ионное уравнение наиболее вычислительно сложно, то его стоит решать как можно меньшее количество раз, с максимально возможным шагом по времени. Условие Куранта для электронов из-за большей скорости более жесткое, чем для ионов, поэтому приходится производить N расчетов {И « 20 — 40 в двумерном моделировании, N и 30 — 90 в трехмерном моделировании в зависимости от параметров) с меньшим шагом по времени для того, чтобы сохранить устойчивость разностных схем. Уравнение Пуассона пересчитывается, так как при движении электронов меняется электрическое поле.

В §1.3 рассмотрены результаты тестирования, изложенных в §1.2 методах решения уравнений. Пункты 1.3.1, 1.3.2, 1.3.3 содержат описание результатов тестирования и сравнение с другими методами решения уравнения для электронной температуры и плотности, уравнения Пуассона и кинетического уравнения для ионов соответственно. В пункте 1.3.4 исследована зависимость численного решения двумерной модели неустойчивости от параметров пространственно-временной сетки. Как оказалось, основные принципы развития неустойчивости не изменяются при разумном изменении количества точек в области моделирования и шагов по времени. При этом наблюдается соответствие результатов ранее проведенным методом частиц численным расчетам фарлей-бунемановской неустойчивости.

Выводы по первой главе сформулированы в §1.5.

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

В §2.1 изучается эффективность распараллеливания алгоритмов

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

В §2.2 проводится анализ результатов двумерного моделирования фарлей-бунемановской неустойчивости. В пункте 2.2.1 исследуется влияние изменения электронной массы, которое обычно используется в методе частиц, на результаты моделирования. Проведенные численные расчеты с реальной и увеличенной электронной массой показали качественное соответствие результатов, однако количественные характеристики развившейся неустойчивости отличаются. Увеличение массы электронов приводит к уменьшению значений турбулентного электрического поля и к увеличению длины волны. Сделан вывод, что метод частиц, использующий модифицированные массы, может использоваться для качественных предсказаний и допускает погрешности в количественных предсказаниях. В пункте 2.2.2 описываются результаты двумерного моделирования в предположении изотермичности электронов, на основе четырех расчетов для различных параметров ионосферы анализируются основные особенности развития фарлей-бунемановской неустойчивости. Использование подхода, альтернативного методу частиц, позволило провести расчеты со значением внешнего поля, близким к пороговому. В соответствии с теоретическими ожиданиями уменьшение внешнего поля приводит к значительному увеличению времени развития неустойчивости, уменьшению амплитуд колебаний плазмы и увеличению длин волн в квазистационарном состоянии. В пункте 2.2.3 с помощью сравнения результатов трех серий моделирования с изотермичными и неизотермичными электронами исследуется влияние электронных тепловых эффектов на развитие фарлей-бунемановской неустойчивости. Сделан вывод, что предположение неизотермичности электронов позволяет добиться результатов более близких к наблюдениям в ионосфере Земли. Электронные тепловые эффекты влияют на неустойчивость следующим образом: уменьшается амплитуда колебаний плазмы и величина возмущенного электрического

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

Выводы по второй главе сформулированы в §2.3.

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

В §3.1 описываются основные приемы распараллеливания алгоритмов моделирования неустойчивости с помощью совместного использования технологий ОрепМР и MPI для выполнения на современных многопроцессорных вычислительных комплексах с распределенной памятью. В пункте 3.1.1 рассматривается оптимизация алгоритмов решения ионного уравнения, в пункте 3.1.2 - электронного уравнения, а в пункте 3.1.3 - уравнения Пуассона.

В §3.2 проводится анализ времени выполнения оптимизированных алгоритмов моделирования фарлей-бунемановской неустойчивости в зависимости от количества используемых MPI узлов и ОрепМР нитей. Тестирование алгоритмов проводится на вычислительных комплексах СКИФ МГУ и IBM Blue Gene/P. Показано, что распараллеленный программный комплекс имеет хорошую масштабируемость и может эффективно выполняться на вычислительных системах IBM Blue Gene/P и СКИФ МГУ "Чебышёв". Пункт 3.2.1 посвящен анализу времени выполнения алгоритмов решения ионного уравнения, пункт 3.2.2 - электронного уравнения и пункт 3.2.3 - уравнения Пуассона. В пункте 3.2.4 исследуется время выполнения программного комплекса в целом на вычислительных комплексах СКИФ МГУ и IBM Blue Gene/Р. Анализ работы оптимизированных алгоритмов для двумерной задачи проводится

в пункте 3.2.5.

В §3.3 описываются результаты двумерного моделирования на большой сетке (более 2000 точек вдоль каждого из направлений в геометрическом пространстве) для изотермичных и неизотермичных электронов, проводится сравнения с расчётами, проведенными на меньших сетках. Установлено, что среднеквадратическое значение турбулентного электрического поля и среднее значение амплитуд колебаний плазмы более стационарны по времени в моделировании на большой сетке, чем в расчетах на малых сетках. Проведен спектральный анализ колебаний плазмы. Результаты моделирования на большей сетке подтверждают ранее сделанные выводы о том, что электронные тепловые эффекты уменьшают амплитуду колебаний плазмы и замедляют развитие неустойчивости.

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

Выводы по третьей главе сформулированы в §3.5.

В заключении изложены основные результаты диссертации.

Основные результаты работы

• Разработана гибридная модель неустойчивости, возникающей в плазме Е-области ионосферы Земли (фарлей-бунемановской неустойчивости), состоящая из кинетического уравнения для функции распределения ионов, зависящего от шести независимых переменных: трех в геометрическом пространстве, двух в пространстве скоростей и от времени, а также из нелинейных гидродинамических уравнений для электронной

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

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

На основе предложенного метода создан программный комплекс для моделирования фарлей-бунемановской неустойчивости, который с помощью технологий MPI и ОрепМР оптимизирован для выполнения на современных массивно-параллельных вычислительных системах. Показано, что программный комплекс имеет хорошую масштабируемость и может эффективно выполняться на вычислительных системах IBM Blue Gene/P и СКИФ МГУ "Чебышёв".

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

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

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

Публикации

Основные результаты диссертации опубликованы в работах:

1. Kovalev D.V., Smirnov А.Р., Dimant Y.S. Modeling of the Farley-Buneman instability in the E-region ionosphere: a new hybrid approach // Ann. Geophys. - 2008. - Vol.26. - N9. - P. 2853-2870.

2. Ковалёв Д.В. Моделирование фарлей-бунемановской неустойчивости с использованием четырехмерного кинетического уравнения Ц Мат. моделирование - 2008. - Т.20. - №12. - С. 89-104.

3. Ковалёв Д.В., Смирнов А.П., Димант Я.С. О влиянии изменения электронной массы в численных расчетах фарлей-бунемановской неустойчивости // Вестник МГУ Сер. 15 Вычисл. матем. и киберн. — 2009. - Т.ЗЗ. - №1. - С. 19-26.

4. Ковалёв Д.В., Смирнов А.П., Димант Я.С. Исследование кинетических эффектов, возникающих при моделировании фарлей-бунемановской неустойчивости // Физ. плазмы — 2009. — Т.35. — №5. — С. 465-471.

5. Ковалёв Д.В., Смирнов А.П., Димант Я.С. Моделирование нелинейного развития фарлей-бунемановской неустойчивости с учетом электронных тепловых эффектов // Физ. плазмы — 2009. — Т.35. — №7. — С. 657-664.

6. Smirnov А.Р., Kovalev D.V., Dimant Y.S. Simulations of Farley-Buneman instability: new hybrid approach // Proceedings of IUGG XXIV General Assembly, Perugia, Italy, 2007. - JMS023, 1376 -http://www.iugg2007perugia.it/webbook/pdf/JM .pdf

7. Dimant Y.S., Kovalev D.V., Smirnov A.P. Numerical simulation of the Farley-Buneman instability: new hybrid approach // Proceedings of URSI 2007 CNC/USNC North American Radio Science Meeting, Ottawa, Canada, 2007.

8. Kovalev D.V., Smirnov A.P., Dimant Y.S. Hybrid-model simulations of Farley-Buneman instability with electron thermal effects // Book .of abstracts. 12th International symposium on equatorial aeronomy Creete, Greece, 2008. — http://iseal2.physics.uoc.gr/files/ISEA12-Book_of_abstracts.pdf

9. Ковалёв Д.В. О методах трёхмерного моделирования фарлей-бунемановской неустойчивости на современных многопроцессорных вычислительных комплексах // Сборник тезисов XVI Международной научной конференции студентов, аспирантов и молодых ученых "Ломоносов-2009". Секция "Вычислительная математика и кибернетика", Москва, Россия, 2009 — С. 35.

Подписано в печать: 07.08.2009

Заказ № 2351 Тираж - 100 экз. Печать трафаретная. Типография «11-й ФОРМАТ» ИНН 7726330900 115230, Москва, Варшавское ш., 36 (499) 788-78-56 www.autoreferat.ru

Введение

Глава 1. Постановка задачи и методы её решения

§ 1.1. Математическая модель неустойчивости.

1.1.1. Условия развития неустойчивости

1.1.2. Гидродинамическая модель электронов.

1.1.3. Уравнение изменения температуры электронов

1.1.4. Кинетическая модель ионов.

1.1.5. Электростатический потенциал.

1.1.6. Двужидкостная теория и пороговые значения для развития неустойчивости

§ 1.2. Методы моделирования.

1.2.1. Уравнения для электронной плотности и температуры

1.2.2. Уравнение Пуассона.

1.2.3. Кинетическое уравнение для функции распределения ионов.

1.2.4. Полный алгоритм решения уравнений неустойчивости на одном шаге по времени.

§ 1.3. Тестирование алгоритмов.

1.3.1. Уравнение для электронной плотности и температуры, двумерный вариант

1.3.2. Уравнение для электронной плотности и температуры, трехмерный вариант.

1.3.3. Уравнение для ионов

1.3.4. Зависимость численного решения двумерной модели от выбора параметров сетки

§ 1.4. Выводы.

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

§ 2.1. Эффективность распараллеливания алгоритмов решения двухмерной модели.

§ 2.2. Анализ результатов двумерного моделирования.

2.2.1. Влияние изменения электронной массы на развитие неустойчивости

2.2.2. Моделирование с изотермическими электронами

2.2.3. Моделирование с неизотермическими электронами

2.2.4. Кинетические эффекты.

2.2.5. Обсуждение применимости результатов моделирования

§ 2.3. Выводы.

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

§ 3.1. Основные приемы распараллеливания программного кода

3.1.1. Алгоритмы параллельного решения кинетического уравнения для ионов

3.1.2. Алгоритмы параллельного решения уравнения для электронов

3.1.3. Алгоритмы параллельного решения уравнения

Пуассона.

§ 3.2. Оценка эффективности параллельной программы

3.2.1. Кинетическое уравнение для ионов

3.2.2. Уравнение для плотности электронов

3.2.3. Уравнение Пуассона.

3.2.4. Производительность программного комплекса для трехмерной задачи.

3.2.5. Производительность программного комплекса для двумерной задачи

§ 3.3. Анализ результатов двумерного моделирования.

§ 3.4. Анализ результатов трехмерного моделирования

§ 3.5. Выводы.

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

Фарлей-бунемановская неустойчивость - это низкочастотная неустойчивость в плазме, движимая достаточно сильным квазистационарным электрическим полем Ео, перпендикулярным геомагнитному полю Во [1—3]. Неустойчивость возникает в слабоионизированной Е-области ионосферы Земли, где электроны замагничены, в то время как ионы из-за частых столкновений с нейтральными частицами практически не подвержены влиянию В0. Ниже Е-области (примерно 95 км) неустойчивость не развивается из-за столкновительного затухания, а на высотах более 120 км ионы становятся замагниченными. В Е-области из-за частых столкновений незамагниченные ионы сильно подвержены диффузии и не могут разгоняться, поэтому при превышении постоянным электрическим полем порогового значения Ец1Г — (Ю — 20) мВ/м, что обычно происходит в экваториальных и приполярных электроструях (электроструя - холловский ток в плазме, характерный для Е-области), функция распределения скоростей электронов смещается относительно ионной функции распределения на величину Т/о — с|Ео х Во|/#о- Кроме низкочастотных электростатических колебаний возникают возмущения плотностей ионов и электронов с /сц здесь и к± - компоненты турбулентного волнового вектора параллельные и перпендикулярные Во соответственно. Средние амплитуды колебаний плотностей обычно не превосходят нескольких процентов. Типичные длины волн колебаний, которые были зафиксированы радарами, лежат в метровом диапазоне. Помимо земной ионосферы фарлей-бунемановская неустойчивость может развиваться в атмосферах других планет, на звёздах [4,5], а также в других разновидностях плазмы, где электроны замагничены, а ионы незамагничены.

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

Волновые возмущения, создаваемые неустойчивостью, фиксируются радарами и приводят к помехам в радиоэфире. Именно поэтому возмущения плазмы в Е-области ионосферы были открыты вскоре после изобретения радара в 40ых годах XX века, а детальные исследования начались уже в 50х [6]. С помощью радаров удалось выявить два типа основных сигналов, исходящих от экваториальных электроструй: один с узким, второй с широким спектром [7,8]. Считается, что эти сигналы возникают из-за волн фарлей-бунемановской неустойчивости или длинноволновых градиентно-токовых колебаний, которые в результате нелинейного взаимодействия преобразуются в более короткие волны и детектируются радарами [9,10]. При дальнейшей изучении радарных наблюдений были найдены и классифицированы и некоторые другие типы сигналов [11].

Кроме радаров для изучения Е-области ионосферы используются ракетные запуски. Это достаточно дорогой тип эксперимента, однако он необходим для калибровки и подтверждения теорий и радарных наблюдений. Ранние эксперименты [12,13] подтвердили существование фарлеп-бунемановских волн и соответствие их основных свойств (распределение по высоте, уровни турбулентности электронов и значения фазовой скорости) теоретическим ожиданиям. В большинстве случаев колебания плазмы составляли 2-3% и практически всегда меньше 10%.

Лабораторные эксперименты по изучению фарлей-бунемановской неустойчивости проводились в цилиндрических установках с радиально направленными электростатическими полями, в которых плазма вращалась за счет Е хВ дрейфа. Полученные результаты подтверждали основные выводы линейной теории: наличие порога неустойчивости и преобладание в движении Е х В дрейфа плазмы. Из более поздних экспериментов [14,15] следовало, что в коротковолновом диапазоне не содержится достаточно энергии для аномального нагрева воли и то, что фазовая скорость волн лежит между скоростью звука и скоростью, предсказанной линейной теорией.

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

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

- параметр, характеризующий высоту в ионосфере, будет описан далее, щп - частота столкновения ионов и нейтральных частиц, У^ = - дрейфовая скорость, а С"1 = ТЛт'Тк ' СК0Р0СТЬ звука. Из соотношения для мнимой части циклической частоты в (1) следует, что плазма становится неустойчивой и рост волн начинается при Уа > С8{ф + 1). Более детальное описание можно найти в [2,10].

Теория Бунемана также предсказывает, что темпы роста неустойчивости пропорциональны значению волнового числа, в то время, как кинетический подход утверждает, что максимальное значение темпов роста неустойчивости достигается при конечном значении волнового числа. Последнее более достоверно, иначе бы неустойчивость приводила к взрывоподобным эффектам в ионосфере. Для того, чтобы показать, что кинетическое описание, используемое в модели Фарлея, предсказывает максимальные темпы роста для волн с конечной длиной в [17] была использована кинетическая модель для ионов и гидродинамическая модель для электронов. В [18] на основе сравнения результатов [17] и |19[ было показано, что доминирующие длины волн в радарных экспериментах совпадают с теми длинами, для которых был предсказан рост в [17]. С помощью расширения линейной кинетической теории в [20] было продемонстрировано, что максимальные темпы роста наблюдаются для волн, у которых значение компоненты волнового вектора, параллельной магнитному полю Во, мало. Авторы [21] измерили параллельную магнитному полю компоненту волнового вектора для наиболее выраженных волн в Е-области и убедились, что эти компоненты действительно малы. Тем не менее, до сих пор нельзя четко сказать является ли это линейным или нелинейным эффектом [22].

Исследования, основанные на линейной теории, продолжаются до сих пор. Теория расширяется за счет усложнения столкновительных операторов, что приводит к некоторым изменениям в дисперсионном соотношении (см., например, пункт 1.1.6 или [11,16]). Но применение линейной теории ограничено. В частности, с помощью неё можно получить корректные пороговые значения неустойчивости, описать начальные стадии развития неустойчивости и вывести некоторые закономерности, такие как зависимость нерегулярностей от электрического поля. Линейная теория сама по себе не может описать процесс нелинейного насыщения, амплитуды и спектральные характеристики развившейся турбулентности. Для объяснения наблюдений и для количественных предсказаний необходима нелинейная теория, которая должна ответить на множество важных вопросов, таких как: какие основные факторы приводят к нелинейному насыщению неустойчивости, какой уровень колебаний плотностей и турбулентного электрического поля, каковы преобладающие длины волн, фазовые скорости и спектральные характеристики в турбулентной плазме в зависимости от величины электрического поля Е0 и других параметров ионосферы.

Ответы на вышеперечисленные вопросы важны не только для объяснения возмущений в ионосфере. Такие факторы как величина колебаний турбулентного поля и средние значения /сц критичны для описания сильного нагрева электронов, наблюдаемого в приполярных электроструях при наличии довольно сильного постоянного внешнего электрического поля [23-25]. Простые измерения показывают, что этот эффект не описывается с помощью нагрева из-за наличия внешнего электрического поля Ец, необходим дополнительный разогрев за счет турбулентного поля [26]. Конечные значения /сц, соответствующие довольно большим значениям турбулентного электрического поля, могут иметь ключевое значение для этого нагрева [27]. В ранних моделях аномального нагрева электронов [28, 29] получались довольно разумные количественные оценки, однако возникали проблемы физической состоятельности модели. Более состоятельная модель аномального нагрева описана в [30,31]. Эта модель основана па эвристическом описании развившейся фарлей-бунемановской неустойчивости, поэтому для её улучшения необходимо более детальное изучение механизмов нелинейного насыщения неустойчивости.

Интересен также вопрос о величине фазовой скорости. Из наблюдений следуют противоречивые результаты (см. ссылки на статьи, описывающие радарные эксперименты, и [32,33]). В соответствии с общепринятой двужидкостной теорией фазовая скорость линейно неустойчивых волн должна быть пропорциональна косинусу угла между направлениями дрейфовой скорости Уц и волнового вектора к, достигая минимального значения, равного скорости звука для ионов С3, на границе зоны линейной устойчивости. Радарные наблюдения в экваториальном регионе показывают, что независимо от направления волнового вектора фазовая скорость наблюдаемых волн часто близка к ионной скорости звука [34]. В то же время, в приполярных электроструях, особенно при больших значениях внешнего электрического поляЕо, радарные наблюдения могут показывать значительные отличия между фазовой скоростью и скоростью звука [35], ближе к тому, что предсказывает линейная теория. Из ракетных экспериментов в экваториальных и полярных регионах [13, 19, 36] следует, что фазовая скорость наиболее ярко выраженных волн может быть сравнима со скоростью, предсказываемой общепринятой линейной теорией. Стоит отметить, что и на радарные, и на ракетные измерения влияют погрешности инструментов, сильная изменчивость и непредсказуемость ионосферы. Возможно, что результаты наблюдений не противоречат друг другу, а соответствуют различным условиям в ионосфере и характеризуют различные группы волн.

Создание общей нелинейной теории фарлей-бунемановской неустойчивости должно помочь разрешить эти противоречия.

Для описания насыщения и нелинейных характеристик неустойчивости использовались различные приемы. Так в [37,38] были предложены одномерные модели нелинейного поведения фарлей-бупемановских волн. Ни одна из работ не учитывала нелинейный дрейф электронов, который имеет ключевое значение в нелинейной теории. Автор [39] указал на этот недостаток и предположил, что короткие волны нелинейно возмущают орбиты электронов, вызывая "аномальную" диффузию, которая приводит к нелинейной коррекции фазовой скорости волн и насыщению самой неустойчивости. В [29] было отмечено, что нелинейный анализ, проведённый [37] и [39] некорректен, поэтому была продолжена работа над новыми теориями. Результаты первых компьютерных моделирований неустойчивости привели к появлению работ, в которых нелинейное насыщение волн описывалось с помощью небольшого количества сильно связанных мод [40-42]. Авторам этих работ удалось довольно точно описать основные особенности наиболее ярко выраженных возмущений в Е-области, которые регистрируются радарами. Тем не менее, до сих пор не понятно, когда моды с малым значением компоненты вдоль магнитного поля играют важную роль, а когда нет.

Первое компьютерное моделирование фарлей-бунемановской неустойчивости на ЭВМ описано в [43]. Это был численный эксперимент, основанный на гидродинамической модели с коэффициентом вязкости, зависящим от длины волны, что в какой-то мере позволяло описывать важный кинетический эффект затухания Ландау. Однако, из-за численных проблем применимость результатов была ограничена. Следующее моделирование неустойчивости описано в [44]. Авторы использовали метод частиц в ячейке для нахождения наиболее корректного механизма нагрева электронов в приполярных электроструях. На основании моделирования авторы пришли к выводу, что наиболее корректный механизм предложен в работах [28,39,45]. Работа [44] была расширена и дополнена в [46], также были проведены сравнения результатов моделирования и экспериментальных наблюдений. Все вышеперечисленные моделирования имели существенный недостаток, так как рассматривали волны в плоскости, параллельной магнитному полю Во и внешнему электрическому полю Ео, полностью игнорируя важный и большой по значению нелинейный член, возникающий из-за дрейфа электронов в плоскости, перпендикулярной Во. Из-за этого амплитуды и характеристики волн при насыщении неустойчивости значительно отличались от значений, получаемых в численных экспериментах с учетом нелинейного дрейфа, и их подобие результатам экспериментальных наблюдений было сомнительно.

В [47] были описаны результаты моделирования методом частица в ячейке, но уже в плоскости, перпендикулярной магнитному полю. Автор наблюдал формирование и движение фарлей-бунемановских волн, хотя и не смог достичь их насыщения из-за недостатка вычислительной мощности. Основной результат, который он наблюдал - это поворот волн так, что (к • Ео) < 0, где к - волновой вектор. По мнению автора насыщение неустойчивости было вызвано этим поворотом волн [48]. Проведенные моделирования не показывали нагрева волн из-за кинетического поведения электронов. Хотя выбор метода моделирования в [47] привёл к появлению ряда серьёзных проблем (см. пункт 2.2.1 и [11]), но эта работа положила начало сложным нелинейным моделированиям фарлей-бунемановской неустойчивости.

Дальнейшее изучение фарлей-бунемановской неустойчивости с помощыо компьютерного моделирования, помимо численных экспериментов, представленных в данной диссертационной работе, проводилось одним автором, сначала на основе гибридного подхода, включавшего метод частиц в ячейке и гидродинамические уравнения в предположении квазинейтральности, так и без гидродинамических уравнений. Моделирование на основе гибридной модели в плоскости, перпендикулярной магнитному полю, описано в [49,50]. Большинство полученных результатов было подтверждено и более поздними численными экспериментами. Основные выводы из моделирования были следующие: нелинейная динамика электронов приводит к насыщению фарлей-бунемаповских волн с амплитудами, которые сравнимы с амплитудами, наблюдаемыми в экспериментах; фазовые скорости меньше значения, предсказываемого линейной теорией, но выше скорости звука; фазовые скорости имесют слабую зависимость от высоты в ионосфере; рост волн наблюдается в направлении отличном от направления Ео х Во дрейфа электронов. Также в [49-51] отмечается тот факт, что фарлсй-бунемановская неустойчивость приводит к нелинейному сдвигу потоков плазмы в ионосфере Земли, что оказывает существенное влияние на глобальную динамику плазмы в Е-области ионосферы во время периодов сильной электромагнитной активности.

В [11] на основании моделирования с помощью метода частиц в ячейке для электронов и ионов были проведены исследования значимости кинетического описания ионов. Было показано, что эффекты ионной температуры играют важную роль в формировании и насыщении волн в Е-области ионосферы, особенно набольших высотах (выше 110 км) и при больших значениях внешнего электрического поля, которое часто встречается на высоких широтах. Эти эффекты влияют на порог неустойчивости, темпы роста и нелинейную динамику [16]. Хотя поворот волн во время развития неустойчивости, наблюдался еще в [47,49,50], он не связывался с термическими эффектами. Более того, моделирование [47] из-за сильной модификации основных параметров, в результате которой были получены соотношения нехарактерные для Е-области [11], больше описывают тепловую неустойчивость, а не фарлей-бунемановскую. Это также означает, что эффект насыщения неустойчивости на основе поворота волн, описанный в [48], имеет малое значения для турбулентности в Е-области.

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

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

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

§3.5. Выводы

• Предложены методы распараллеливания алгоритмов двумерного и трехмерного моделирования фарлей-бунемановской неустойчивости с помощью совместного использования технологий параллельного программирования ОрепМР и MPI.

• Показано, что распараллеленный программный комплекс имеет хорошую масштабируемость и может эффективно выполняться на вычислительных системах IBM Blue Gene/P и СКИФ МГУ "Чебы-шёв".

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

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

Заключение

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

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

• На основе предложенного метода создан программный комплекс для моделирования фарлей-бунемановской неустойчивости, который с помощью технологий MPI и ОрепМР оптимизирован для выполнения на современных массивно-параллельпых вычислительных системах. Показано, что программный комплекс имеет хорошую масштабируемость и может эффективно выполняться на вычислительных системах IBM Blue Gene/P и

СКИФ МГУ "Чебышёв".

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

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

1. D.T. Farley. A plasma instability resulting in field-aligned irregularities in the ionosphere // J. Geophys. Res. — 1963. — Vol. 68. — Pp. 6083-6097.

2. O. Buneinan. Excitation of field aligned sound waves by electron streams // Phys. Rev. Lett. 1963. - Vol. 10. - Pp. 285 -288.

3. D.V. Kovalev, A.P. Smirnov, Y.S. Dimant. Modeling of the Farley-Buneman instability in the E-region ionosphere: и new hybrid approach // Ann. Geophys. 2008. - Vol. 26, no. 9. -- Pp. 2853-2870.

4. V.A. Liperovsky, C.-V. Meister, E.V. Liperovskaya, K.V. Popov, S.A. Senchenkov. On the generation of modified low-frequency Farley-Buneman waves in the solar atmosphere // Astronomische Nachrichten. 2000. - Vol. 321, no. 2. - Pp. 129-136.

5. J.M. Fontenla,, W.K. Peterson, J. Harder. Chromospheric heating by the Farlcy-Bunoman instability /'/ Astronomy and Astrophysics. — 2008. —-Vol. 480, no. 3. - Pp. 839-846.

6. K.L. Bowles. Doppler shifted radio echoes from aurora // J. Geophys. Res. 1954. - Vol. 59, no. 4. - Pp. 553-555.

7. R. Cohen, K.L. Bowles. Secondary irregularities in the equatorial elec-trojet // J. Geophys. Res. 1967. - Vol. 72, no. 3. - Pp. 885-894.

8. B.B. Balsley. Some characteristics of non-two-steam irregularities in theequatorial electrojet // J. Geophys. Res. — 1969. — Vol. 74. no. 9. — Pp. 2333-2347.

9. K. Maeda, T. Tsuda, H. Maeda. Theoretical interpretation of the equatorial sporadic E layers // Phys. Rev. Lett. — 1963. — Vol. 11, no. 9. — Pp. 406-409.

10. R.N. Sudan. J. Akinrimisi, D.T. Farley. Generation of small-scale irregularities in the equatorial electrojet // J. Geophys. Res. — 1973. — Vol. 78, no. 1. Pp. 240-248.

11. M.M. Oppenheim, Y.S. Dimant. Ion thermal effects on E-region instabilities: 2d kinetic simulations // J. Atmos. Terr. Phys. 2004. — Vol. 66, no. 17. - Pp. 1655-1668.

12. M.C. Kelley, F. Mozer. Electric field and plasma density oscillations due to the high-frequency Hall current two-stream instability in the auroral E region // J. Geophys. Res. 1973. - Vol. 13. - Pp. 2214-2221.

13. B. Kustom, N. D'Angelo, R. Merlino. A laboratory investigation of the high-frequency Farley-Buneman instability // J. Geophys. Res. — 1985. Vol. 90, no. A2. - Pp. 1698-1704.

14. Y.S. Dimant. M.M. Oppenheim. Ion thermal effects on E-region instabilities: linear theory // J. Atmos. Terr. Phys. ~ 2004. — Vol. 66, no. 17. Pp. 1639-1654.

15. M.J. Schmidt, S.P. Gary. Density gradients and the Farley-Buneman instability // J. Geophys. Res. 1973. - Vol. 78, no. 34. - Pp. 82618265.

16. M.C. Kelley. The Earth's Ionosphere. — San Diego: Academic, 1989.

17. P.F. Pfaff, M.C. Kelley, E. Kudeki, et. a.l. Electric field and plasma density measurements in the strongly driven daytime equatorial electrojet. 2. Two-stream waves // J.Geophys. Res. 1987. -- Vol. 92, no. A12.- Pp. 13597—13612.

18. S.L. Ossakow, K. Papadopoulos, J. Orens, T. Coffey. Parallel propagation effects on the type 1 electrojet instability // J. Geophys. Res. — 1975. Vol. 80, no. 1. - Pp. 141-148.

19. E. Kudeki, D.T. Farley. Aspect sensitivity of equatorial electrojet irregularities and theoretical implications // J. Geophys. Res. — 1989. — Vol. 94, no. Al. Pp. 426-434.

20. B.J. Jackel. D.R. Moorcroft, K. Schlegel. Characteristics of very large aspect angle E-region coherent echoes at 933 MHz // Ann. Geophys. — 1997. Vol. 15, no. 1. - Pp. 54 -62.

21. K. Schlegel, J.-P. St.-Maurice. Anomalous heating of the polar E region by unstable plasma waves, 1, Observations //' J. Geophys. Res. — 1981.- Vol. 86, no. A3. Pp. 1447-1452.

22. J.-P. St.-Maurice, K. Schlegel, P.M. Banks. Anomalous heating of the polar E region by unstable plasma waves. II Theory //J. Geophys. Res. - 1981. - Vol. 86, no. A3. - Pp. 1453-1462.

23. J.-P. St.-Maurice, R. Laher. Are observed broadband plasma wave amplitudes large enough to explain the enhanced electron temperatures of the high-latitude E region? // J. Geophys. Res. 1985. Vol. 90, no. A3. - Pp. 2843-2850.

24. T.R. Robinson. Towards a self-consistent non-linear theory of radar auroral backscatter /'/' J. Atmos. Terr. Phys. — 1986. — Vol. 48, no. 5. — Pp. 417-422.

25. J.-P. St.-Maurice. A unified theory of anomalous resistivity and Joule heating effects in the presence of ionospheric E region irregularities // J. Geophys. Res. 1987. - Vol. 92, no. A5. - Pp. 4533-4542.

26. Y.S. Dimant, G.M. Milikh. Model of anomalous electron heating in the E region: 1. Basic theory // J. Geophys. Res. 2003. - Vol. 108, no. A9. - P. 1350.

27. G.M. Milikh, Y.S. Dimant. Model of anomalous electron heating in the E region: 2. Detailed numerical modeling // J. Geophys. Rus. 2003. Vol. 108, no. A9. - P. 1351.

28. C. Haldoupis. K. Schl^gel, G.C. Hussey, ot al. Radar observation of kinetic effects at meter scales for Farley-Buneman plasma waves //J. Geophys. Res. 2002. - Vol. 107, no. A10. - P. 1272.

29. E.E. Timofeev. M.K. Vallinkoski, P. Pollari, et al. Flow angle dependence of 1-m ionospheric plasma wave turbulence for near-threshold radar echo electric fields // J. Geophys. Res. 2002. - Vol. 107, no. A10. -P. 1286.

30. D.T. Farley. Theory of equatorial electrojet plasma waves: new developments and current status /'/' J. Atrnos. Terr. Phys. — 1985. — Vol. 47. no. 8-10. Pp. 729-744.

31. J.C. Foster. D. Tetenbaum. Phase velocity studies of 34-cm E-region irregularities observed at Millstone Hill // </. Atmos. Terr. Phys. — 1992. Vol. 54, no. 6. - Pp. 759 -768.

32. K. Rinneit. Plasma waves observed in the auroral E-region ROSE campaign // J. Atmos. Terr. Phys. — 1992. — Vol. 54, no. 6. — Pp. 683692.

33. A. Rogistcr. Nonlinear theoiy of cross field instability with application to the equatorial electrojet // J. Geophys. Res. — 1972. — Vol. 77, no. 16. Pp. 2975-2981.

34. J. Weinstock. A. Sleeper. Nonlinear saturation of type I irregularities in the equatoiial electrojet /'/ J. Geophys. Res. 1972. - Vol. 77, no. 19. - Pp. 2975-2981.

35. R.N. Sudan. Theory nonlinear theory of type I irregularities in the equatorial electrojet // Geophys. Res. Lett. — 1983. — Vol. 10, no. 10. — Pp. 983-986.

36. J.D. Sahr, D.T. Farley. Three-wave coupling in the auroral E-region /'/ Ann. Geophys. 1995. - Vol. 13, no. 1. - Pp. 38-44.

37. N.F. Otani, M.M. Oppenheim. A saturation mechanism for the Farley-Buneman instability // Geophys. Res. Lett. — 1998. — Vol. 25, no. 11.- Pp. 1833-1836.

38. Y.S. Dimant. Nonlinearly saturated dynamical state of a three-wave mode-coupled dissipative system with linear instability // Phys. Rev. Lett. 2000. - Vol. 84, no. 4. - Pp. 622-625.

39. A.L. Newman, E. Ott. Nonlinear simulations of type I irregularities in the equatorial electrojet // J. Geophys. Res. — 1981. — Vol. 86, no. A8.- Pp. 6879-6891.

40. S. Machida, C.K. Goertz. Computer simulation of the Farlcy-Bunenian instability and anomalous electron heating in the auroral ionosphere // J. Geophys. Res. 1988. - Vol. 93, no. A9. - Pp. 9993-10000.

41. F. Prirndahl. Polar ionospheric E region plasma wave stabilization and electron heating by wave-induced enhancement of the electron collision frequency // Phys. Scr. 1986. - Vol. 33, no. 2. - Pp. 187-191.

42. K. Schlegel, H. Thiemann. Particle-in-cell plasma simulations of the modified two-stream instability /'/' Ann. Geophys. — 1994. — Vol. 12, no. 10-11. Pp. 1091-1100.

43. P. Janhunen. Perpendicular particle simulation of the E region Farley-Buneman instability // ,J. Geophys. Res. — 1994. — Vol. 99, no. A6. — Pp. 11461-11473.

44. P. Janhunen. implications of flow angle stabilization un coherent Eregion spectra // J. Geophys. Res. — 1994. — Vol. 99, no. A7. — Pp. 13203-13208.

45. M.M. Oppenheim, N.F. Otani. Spectral characteristics of the Failey-Buneman instability: simulations versus observations // J. Geophys. Res. 1996. - Vol. 101, no. All. - Pp. 24573-24582.

46. M.M. Oppenheim, N.F. Otani. G. Ronchi. Saturation of the Farley-Buneman instability via nonlinear electron E x В drifts // J. Geophys. Res. 1996. - Vol. 101, no. A8. - Pp. 17273-17286.

47. M.M. Oppenheim. Evidence and effects of a wave-driven nonlinear current in the equatorial electrojet // Ann. Geophys. — 1997. — Vol. 15, no. 7. • Pp. 899 907.

48. M.M. Oppenheim, Y. Dimant, L.P. Dyrud. Large-scale simulations of 2D fully kinetic Farley-Buneman turbulence // Ann. Geophys. — 2008. Vol. 26, no. 3. - Pp. 543-553.

49. Д.В. Ковалёв. Моделирование фарлей-бунемановской неустойчивости с использованием четырёхмерного кинетического уравнения //' Мат. моделирование. 2008. - Т. 20, № 12. - С. 89-104.

50. Д.В. Ковалёв, А.П. Смирнов, Я.С. Димант. О влиянии изменения электронной массы в численных расчетах фарлей-бунемановской неустойчивости // Вестник МГУ\ сер. Вычислительная математика и кибернетика. — 2009. — Т. 33, № 1. — С. 19-26.

51. Д.В. Ковалёв, А.П. Смирнов, Я.С. Димант. Исследование кинетических эффектов, возникающих при моделировании фарлей -бунемановской неустойчивости // Физ. плазмы. — 2009. — Т. 35, № 5. С. 465 -471.

52. Д.В. Ковалёв, А.П. Смирнов, Я.С. Димант. Моделирование нелинейного развития фарлей-бунемановской неустойчивости с учётом электронных тепловых эффектов // Фаз. плаз,мы. — 2009. — Т. 35, № 7. С. 657-664.

53. D.V. Kovalev, А.P. Smirnov, Y.S. Dirnant. Simulations of Farley-Buneman instability: new hybrid approach // Proceedings of IUGG XXIV General Assembly. — Perugia, Italy: 2007.

54. Y.S. Dimant, D.V. Kovalev, A.P. Smirnov. Numerical simulation of the Farley-Buneman instability: new hybrid approach //' Proceedings of UR-SI 2007 CNC/USNC North American Radio Science Meeting. Ottawa, Canada: 2007.

55. D.V. Kovalev, A.P. Smirnov, Y.S. Dimant. Hybrid-model simulations of Farley-Buneman instability with electron thermal effects // Book of abstracts. 12th International symposium on equatorial aeronomy — Creete, Greece: 2008.

56. R.W. Schunk, A.F. Nagy. Ionospheres-Physics, Plasma Physics and Chemistry. — New York: Cambridge University Press.

57. A.B. Гуревич, А.Б. Шварцбург. Нелинейная теория распространения радиоволн в ионосфере. — М.: Наука, 1973.

58. А.А. Самарский. Теория разностных схем. — М.: Наука, 1989.

59. P.L. Bhatnagar, Е.Р. Gross, М. Krook. A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component, systems // Phys. Rev. — 1954. — Vol. 94, no. 3. — Pp. 511525.

60. Y.S. Dimant. R.N. Sudan. Kinetic theory of low-frequency cross-field instability in a weakly ionized plasma. I // Phys. Plasmas. — 1995. — Vol. 2, no. 4. Pp. 1157—1168.

61. Y.S. Dimant, R.N. Sudan. Kinetic theory of the Fariey-Bimoman instability in the E region of the ionosphere // J. Geoph/ys. Res. — 1997. — Vol. 100, no. A8. Pp. 14605-14624.

62. M. Rosenberg, V.W. Chow. Farlcy-Bunernan instability in a dusty plasma /'/ Planet. Space Set. ~ 1998. Vol. 46, no. 1. - Pp. 103-108.

63. B.A. Finlayson. Numerical methods for problems with moving fronts. —■ Seattle: Ravenna Park Publishing, Inc., 1992.

64. R. Liska, В WendrofF. Composite schemes for conservation laws // SIAM J. on Num. Analysis. 1998. - Vol. 35, no. 6. - Pp. 2250 2271.

65. Г.И. Марчук. Методы вычислительной математики. — М.: Наука, 1977.

66. А.В. Кажихов. А.Е. Мамонтов. Об одном классе выпуклых фуикций и точных классах корректности задачи Коши для уравнения переноса в пространствах Орлича // Сибирский матем. журнал. — 1998. Т. 39, № 4. - С. 831-850.

67. R.J. DiPerna, P.L. Lions. Ordinary differential equations, transport theory and Sobolev spaces // Invent, math. — 1989. — Vol. 98, no. 3. — Pp. 511-547.

68. А.А. Самарский, А.В. Гулин. Численные методы математической физики. — М.: Научный мир, 2000.74| W. Press, S. Teulkovsky, W. Vetterling, et al. Numerical recipes in C. Art of scientific computing. — Cambridge: Cambridge University Press.

69. M. Frigo. S.G. Johnson. The FFTW web page. 2005.http://www.fftw.org.

70. M. Frigo, S.G. Johnson. The design and implementation of FFTW3. // Proceedings of the IEEE. 2005.

71. N.A. Kuz'micheva, A.P. Smirnov. A numerical-analytic method for solving Landau's two-dimensional kinetic equation in self-similar variables // Сотр. Maths Math. Phys. 1994. - Vol. 34, no. 6. - Pp. 775784.

72. E. Sonnendrucker, J. Roche, P. Bertrand, A. Ghizzo. The serni-lagrangian method for the numerical resolution of the vlasov equation // J. Comput. Phys. 1999. - Vol. 149, no. 2. - Pp. 201-220.

73. F. Filbet, E. Sonnendrucker. Comparison of eulerian viasov solvurs // Computer Phys. Comm. 2003. - Vol. 150, no. 3. - Pp. 247-266.

74. А.А. Самарский, А.В. Гулин. Численные методы. — М.: Наука, 1989.

75. R.F. Pfaff, J. Sa.hr, J.F. Providakes, et al. The E-region rocket/radar instability study (ERRIS): scientific objectives and campaign overview // ,J. Atmos. Terr. Phys. 1992. - Vol. 54, no. 6. - Pp. 779 -808.

76. К. Lee, C.F. Kennel, J.M. Kindel. High-frequency Hall current instability // Radio Set. 1971. -■ Vol. G, no. 2. -■ Pp. 209-213.

77. T.F. Morse. Kinetic model equations for a gas mixture // Phys. Fluids. 1964. - Vol. 7, no. 12. - Pp. 2012-2013.

78. P. Stubbe. A new collisional relaxation model for small deviations from equilibrium /7 J. Pi. Phys. 1987. - Vol. 38, no. 1. - Pp. 95-116.

79. P. Stubbe. Theory of electrostatic waves in an E region plasma. 1. General formulation // J. Geophys. Res. — 1989. — Vol. 94, no. A5. — Pp. 5303-5315.

80. L.M. Kagan, M.C. Kelley. A thermal mechanism for generation of small-scale irregularities in the ionospheric E region // J. Geophys. Res. — 2000. Vol. 105, no. A3. - Pp. 5291-5302.

81. K. Schlegel, C. Haldoupis. Observation of the modified two-stream plasma instability in the midlatitude E region ionosphere // J. Geophys. Res. 1994. - Vol. 99, no. A4. - Pp. 6219-6226.

82. C. Haldoupis, K. Schlegel, D.T. Farley. An explanation for type 1 radar echoes from the midlatitude E-region ionosphere E-region ionosphere // Geophys. Res. Lett. 1996. - Vol. 23, no. 1. - Pp. 97-100.

83. C. Haldoupis, D.T. Farley. K. Schlegel. Type-1 echoes from the mid-latitude E-region ionosphere // Ann. Geophys. — 1997. — Vol. 15, no. 7. Pp. 908-917.

84. Суперкомпьютер IBM Blue Gene/'P на факультете BMK МГУ.http://hpc.cmc.msu.iu/bgp.

85. Суперкомпьютер СКИФ МГУ. http://paraHel.ru/cluster/b'kifinsu.html.

86. В.Э. Витковский, М.П. Федорук. Вычислительная производительность параллельного алгоритма прогонки на кластерных суперкомпьютерах с распределенной памятью // Выч. методы и программирование. 2008. - Т. 9. - С. 305-310.