Содержание к диссертации
Введение
I. Законы сохранения и полностью консервативные разностные схемы для нелинейного уравнения типа ландау и фоккера-планка
1.1. -Общие закономерности при конструировании консервативных разностных схем
1.2- Формы записи интеграла столкновений Ландау- Фоккера-Планка и полностью консервативная разностная схема
II. Релаксация к равновесию начального распределения газа частиц с дальнодеиствующими степенными потенциалами взаимодействия
47 2.1-Асимптотическое решение уравнения типа Ландау в высокоэнергетичной части распределения 49
2.2- Полностью консервативная разностная схема для изотропного уравнения типа Ландау
54 2.3- Релаксация начального распределения к равновесию для частиц со степенными потенциалами взаимодействия 60
III. Квазистационарные функции распределения для уравнения типа ландау (фоккера - планка) при наличии источников
3.1-Асимптотическое решение уравнения типа Ландау-Фоккера-Планка с источниками, локализованными в высокоэнергетичной области
3.2- Эволюция функции распределения при наличии внешней накачки 78
3.3-Сравнение с экспериментальными результатами по облучению полупроводниковой плёнки 87
IV. Система уравнений ландау-фоккера-планка. релаксация двухтемпературной плазмы 95
4.1- Релаксация двухтемпературной плазмы. Асимптотические оценки
4.2-Полностью консервативная разностная схема для системы изотропных уравнений Фоккера-Планка 105
4.3- Численный расчёт процесса релаксации ионов и электронов для начальных условий далёких от равновесия
4.4-Выводы
V. STRONG Двумерное уравнение ландау - фоккера - планка. релаксация импульса
STRONG 5.1-Полностью консервативная разностная схема для двумерного уравнения Фоккера-Планка
132 5.2- Уравнение Фоккера-Планка для изотропных потенциалов Розенблюта-Трубникова 141 5.3-Релаксация моноэнергетического анизотропного пучка ионов 146
5.4- Сравнение с результатами дискретного статистического моделирования 153
VI. Численные расчеты задач динамики столкновительнои плазмы, находящейся во внешнем электромагнитном поле 160
6.1- Распад и заполнение плазмы в открытой магнитной ловушке 160
6.2- Нагрев и ускорение электронов при воздействии статического и высокочастотного электромагнитных полей 169
6.3- Электронная функция распределения при алъфвеновском нагреве в токомаке и токи увлечения 188
6.4- Ускорение электронов алъфвеновской волной с переменной фазовой скоростью и поведение хоров во время суббури 196
6.5- Аномальное высыпание электронов в авроралъной зоне Земли
Заключение
Литература
- Формы записи интеграла столкновений Ландау- Фоккера-Планка и полностью консервативная разностная схема
- Полностью консервативная разностная схема для изотропного уравнения типа Ландау
- Эволюция функции распределения при наличии внешней накачки
- Численный расчёт процесса релаксации ионов и электронов для начальных условий далёких от равновесия
Введение к работе
Кинетическими уравнениями описывается динамика систем, состоящих из большого числа слабо взаимодействующих частиц. Типичными примерами служат разреженные нейтральные газы и плазма. Наиболее известное кинетическое уравнение - это нелинейное кинетическое уравнение Больцмана [1-3], которое описывает систему многих частиц, взаимодействующих по законам классической механики и является основным уравнением в моделях динамики разреженного газа. В общем виде уравнение для функции распределения частиц /, зависящей от пространственной координаты г, скорости v и времени t > О может быть представлено следующим образом
EL+y2L + I.EL = (?L\ +s.
dt дт mdv \dt J c Здесь m - масса частиц, F -внешняя сила, S - источники (стоки) частиц и
энергии. В правой части уравнения стоит так называемый оператор (интеграл)
столкновений
(9 = J[l f]=/dw/i# иа{щ ^ vww - /(v)/(w)i (1)
Дифференциальное сечение рассеяния сг(и, ц) задаётся как функция модуля относительной скорости и —\ v — w|>0h косинуса угла рассеяния /х = cosQ Є
[-1,1]-
Кроме других классических примеров, таких как уравнение Ландау [4,5] в
физике плазмы, или как его ещё называют в литературе, уравнение Фоккера-
Планка [6], кинетические уравнения играют важную роль в моделировании
гранулированных газов, заряженных частиц в полупроводниках, переносе
нейтронов, динамике народонаселения и во многих других областях (см., например, [7-Ю]).
Самый общий подход к численному решению полного кинетического уравнения основан на его расщеплении по физическим параметрам. Решение на одном шаге по времени At состоит из последовательности двух этапов. Вначале, на этапе столкновений. интегрируется пространственно однородное уравнение (внешняя сила и источники отсутствуют) для всех пространственных переменных с шагом по времени At:
^ = J[f,f], /(r,v,0) = /0(r,v).
Затем интегрируется уравнение переноса (левая часть полного уравнения) для того же шага At, используя функцию распределения, полученную на предыдущем шаге в качестве начального условия:
|[ = Л/,/], /(r,v,0) = /(r,v,Af).
По окончании этих двух этапов процесс может итерироваться для получения численного решения в последующие временные шаги. Первый этап действует только на скорость v, в то время как на втором этапе действие оказывается лишь на пространственную переменную г.
Хотя реализация этой идеи остаётся достаточно нелёгкой задачей, тем не менее, методологически она обладает определённым преимуществом из-за распараллеливания схемы расчёта [9,11-13].
Настоящая работа посвящена разработке численных методов решения пространственно-однородного нелинейного кинетического уравнения (системы уравнений) типа Ландау-Фоккера-Планка (ЛФП) и исследованию на их основе
ряда задач динамики разреженных газов и плазмы. Основное внимание уделяется конструированию численной схемы, которая обладает физическими особенностями исходного уравнения: законы сохранения, положительность функции распределения, Я-теорема, поскольку все эти свойства существенно зависят от моделирования именно столкновительного шага.
Нелинейный интеграл столкновений типа ЛФП является приближением интеграла столкновений Больцмана при учёте рассеяния на малые углы. Кинетическое уравнение ЛФП лежит в основе всех математических моделей, описывающих динамику столкновительной плазмы, и на протяжении десятилетий имеет самое широкое применение в чисто научных и прикладных задачах.
Это уравнение, широко используемое в приложениях, является предельным случаем уравнения Больцмана для экранированного кулоновского взаимодействия [4-6,14]. Модели кинетических процессов, обусловленных кулоновскими столкновениями, занимают значительное место в практических приложениях, связанных с высокотемпературной плазмой, как лабораторной так и магнитосферной, полупроводниковой плазмой, а также в плазмохимических задачах [15 - 21].
В кулоновских столкновениях рассеяние на малые углы вносит основной вклад в столкновительный член. Интеграл столкновений для заряженных частиц впервые был получен Ландау [4] из уравнения Больцмана с учетом эффекта малости передаваемого при кулоновских столкновениях импульса и эффекта экранирования заряда частицы вне сферы дебаевского радиуса другими частицами. Столкновительный член в уравнении Ландау явялется интегро-
дифференциальным оператором и имеет симметричную форму, схожую с формой интеграла столкновений Больцмана:
К*)„=ifdwUi1 (-ё)/(v)/(w)lv'wєл''-'<2)
где Г = 27ге4 L/m2, относительная скорость u = v — w, L - так называемый кулоновский логарифм и симметричное ядро Uij определено как
vr Грубым условием применимости уравнения (2) служит неравенство
е2п1/3 < Т,
означающее, что средняя энергия кулоновскго взаимодействия мала по сравнению со средней кинетической энергией (п - плотность числа частиц, Г - температура, выраженная в энергетических единицах).
Уравнение Ландау (2) было переоткрыто через 20 лет в работе [б] в форме нелинейного уравнения Фоккера-Планка. Другая форма записи того же самого уравнения оказалась очень полезной для численных расчётов, связанных с физикой горячей плазмы. Это уравнение часто называют в литературе уравнением Ландау-Фоккера-Планка (ЛФП), мы используем ту же терминологию. Уравнение ЛФП для n-компонентной плазмы имеет вид
1(ЁЬ) д lfadK | 1 д (f д"9а )} ap = i п (3)
Га \ dt J с dvi \ dvi 2 dvj \ dvidvj J J '
где функции
ha = У" К<*^1 + —) Г dw/^w, t) I v
V та J
wl"1
5 = Yl KaP І dwMw> b) I v ~ w I (4)
в >
- так называемые потенциалы Розенблюта-Трубникова [5,6]. Здесь Га = A-nZj/гПа2, Кар = (Za/Zp) Аар, та, Za - масса и заряд.
Симметричная форма уравнения, полученная Ландау ближе по характеру к оператору столкновений Больцмана. Однако, в этом операторе содержится симметричное ядро Щ (неотрицательная симметричная матрица), требующее при численной реализации большой объем компьютерной памяти даже для двумерного случая. В частности, возможно поэтому форма Фоккера-Планка была более употребительна в численных расчетах. Кроме того, разделение столкновительного интеграла на две части (трение и диффузия) иногда облегчает анализ конкретной математической модели, связанной с физической задачей.
Позже классическое уравнение Ландау было обобщено на случай произвольных потенциалов взаимодействия [22,23]. Подробнее об этой модели будет сказано в главе I.
В отсутствие источников и стоков частиц и энергии для уравнения (2-4) справедливы три закона сохранения: плотности числа частиц, импульса (средней скорости), энергии (температуры):
[ fdv, v = - [ /vdv, Т = -^- [ /(v - v)2dv.
Ум3 njR3 3kBnJ]R3
Для интеграла столкновений справедлива Я—теорема Больцмана. Если определить Я-функцию как
Я(/) = / fln(f)dv,
ін(і) _ , , .,„, lnU)dv 0
JR3 \dt J(
dt ,/КЗ \ VI / c
Я-функция монотонно убывает, достигая своего минимума на распределении Максвелла
m3/2 / т\ v — v |2\
Я—теорема Больцмана означает, что любая равновесная функция распределения, т.е., любая функция, для которой интеграл столкновений (dfa/dt)c — 0,
С математической точки зрения кинетическое уравнение ЛФП представляет собой весьма сложный инструмент для добывания необходимой физической информации и поддаётся аналитическому решению с трудом и при существенных упрощениях. Поэтому даже в линеаризованном и линейном случаях для задач, имеющих практическое значение, используются численные методы решения уравнения ЛФП.
Численные модели кинетических процессов, обусловленных кулоновскими столкновениями, занимают значительное место в приложениях. Существует обширная литература, посвященная этой теме (см., например, [16,17]) . Однако
следует отметить, разработка математических моделей на основе нелинейного кинетического уравнения ЛФП, как и развитие численных алгоритмов, адекватно отражающих основные физические законы, остаётся весьма сложной проблемой.
Необходимо добавить, что нелинейное кинетическое уравнение ЛФП содержит в себе многие основные особенности уравнений физической кинетики. Поэтому актуальность его исследования обусловлена как чисто математическим интересом в пониманиии качественных свойств решений нелинейных кинетических уравнений, так и важностью проблемы количественного описания, связанного с практическими приложениями этих уравнений в теории газов и плазмы.
Классической проблемой теории разностных схем является исследование асимптотических свойств схемы при неограниченном измельчении пространственно-временной сетки. Однако, значительное внимание уделяется аспектам теории, связанным со свойствами схемы, обеспечивающей заданную точность на реальных грубых сетках. Кроме того, практические задачи физики в большинстве своём очень сложные и не покрываются доказанными теоремами. При проведении практических расчётов приходится сочетать физическую интуицию с соображениями здравого смысла, поскольку зачастую трудно найти аналитические тесты для реальных задач. Численное моделирование всё более обращается к комплексным задачам и результаты его трактуются как численный эксперимент, при этом проверкой корректности расчёта служит сравнение с результатами физического эксперимента. Такой подход конечно же не может рассматриваться в качестве надёжной верификации модели. В связи с этим численное решение уравнения, максимально приближенное по своим основным
свойствам к решению точного уравнения, является важнейшим тестом для более сложных моделей. Именно это является одной из целей настоящей работы: решение некоторых базовых моделей и как результат создание своего рода набора тестовых задач (бенчмарк).
Важность корректного решения оператора ЛФП заключается, в частности, в том, что для пространственно неоднородного уравнения зачастую необходимо включать эффект столкновений. В этом случае используются упрощённые модели интегралов столкновений, поскольку полновесное описание столкновительных эффектов слишком сложно и дорого в вычислительном смысле. Тогда применяется, например, модель БГК [24] и проверяется её надёжность. Если рассматривается поглощение (затухание) волн, то этой моделью можно пользоваться в случае, когда надо модифицировать экстремальные эффекты пространственой неоднородности функции распределения в разреженной космической плазме и если специально аккуратно проверять используемые скорости релаксации. Для высоких же моментов функции распределения (выше второго), или реалистического описания хвостовой части распределения эта модель не годится. Надёжность этой модели можно сверить с полным оператором ЛФП (для расчета температуры) и упростить расчет пространственно-неоднородной задачи [25].
При решении физических задач естественно рассматривать разностную схему как дискретную математическую модель реального физического процесса. При этом необходимо потребовать, чтобы модель правильно отражала принципиальные стороны исследуемого процесса. В физике такими базовыми
сторонами являются в первую очередь фундаментальные законы сохранения -массы, энергии, импульса. Разностные схемы, обладающие аналогами законов сохранения, называются консервативными [26-28].
Принцип консервативности как один из основных принципов построения разностных схем был выдвинут в [29,30]. Интегро-интерполяционный метод конструирования консервативных схем широко применяется в настоящее время. Там же был указан пример расходимости неконсервативной схемы в классе разрывных коэффициентов. При численном расчёте практических задач предпочтение отдается консервативным разностным схемам.
Далее принцип консервативности был усилен и сформулирован принцип полной консервативности [31,32]. С физической точки зрения необходимость такого усиления связана с тем, что помимо законов сохранения, могут оказаться важными некоторые балансные соотношения (например, между различными видами энергии - внутренней и кинетической). Дисбаланс вносят фиктивные источники, причина появления которых состоит в несогласованности системы разностных уравнений. В [33] сформулировано следующее правило отбора: разностная схема должна одновременно аппроксимировать различные виды записи исходной системы дифференциальных уравнений, имеющие непосредственный физический смысл. В [34] полная консервативность для уравнений магнитной гидродинамики в двумерном случае была достигнута в результате применения вариационного метода.
Ясно, что сама проблема построения консервативных и полностью консервативных схем связана с тем, что эквивалентные формы записи
дифференциального уравнения (системы) приводят, вообще говоря, к неэквивалентным разностным схемам. Даже из приведенных выше уравнений (2) и (3) видно, что записать уравнение можно в разных формах - это будет по сути одно и то же уравнение. При переходе к дискретному случаю и построению разностных схем мы должны выбрать для аппроксимации конкретный вид уравнения. Как оказывается, это зачастую является атрибутивным условием при наложении определенных требований на модель. В частности, если мы требуем от численной модели таких же макроскопических зависимостей (сохранения инвариантов), что и от точной модели. Именно эта сторона понятия полной консервативности является главной для уравнения типа ЛФП.
Дело в том, что для нелинейных кинетических уравнений типа Больцмана и ЛФП возникает нетривиальная ситуация, когда для одного уравнения справедливы несколько законов сохранения. Разностные схемы для кинетического уравнения, допускающие аналогичные ему эквивалентные представления, будем называть полностью консервативными схемами.
В [11] обсуждается проблема построения консервативного вычислительного алгоритма для кинетического уравнения Больцмана. При решении и других проблем сейчас общепринято использовать модели, сохраняющие инварианты системы при численном расчёте. Симметрии дифференциальных уравнений математической физики являются их неотъемлимым свойством и, следовательно, должны учитываться при построении дискретных аналогов (см.,например, [35,36]).
Наиболее распространенный метод для численного моделирования уравнения
ЛФП - это метод конечных разностей. Существует другой подход к решению уравнений больцмановского типа [37-39], описывающих столкновительные и излучательные процессы в разреженном газе и плазме, ставящие в соответствие системе уравнений математической физики систему стохастических дифференциальных уравнений для скачкообразных марковских процессов. Система уравнений так называемого физико-вероятностного аналога решается численно (см.,например, [40-42]). Однако эффективность (быстродействие) этих методов остаётся не очень высокой.
Для описания дальнодействующих потенциалов взаимодействия применение методов расчёта типа Монте-Карло, в отличие от газов больцмановского типа, обладают существенными трудностями. Лишь недавно появился алгоритм, позволяющий эффективно моделировать статистически уравнение типа ЛФП [43,44]. Пример прямого статистического моделирования уравнения ЛФП будет рассмотрен ниже в диссертации. Методы типа Монте Карло обладают тем преимуществом, что позволяют рассчитывать трёхмерные задачи и в принципе могут объединяться естественно алгоритмически с PIC (particle in cell) методом расчёта бесстолкновительного уравнения Власова. Следует отметить, что методы типа Монте Карло хорошо работают, когда нужно знать лишь несколько первых моментов функции распределения. Для высоких моментов функции распределения и при приближении к стационарному состоянию эти методы дают сильные статистические флуктуации. В частности, хвосты распределения считаются с помощью этих методов плохо. Так что при выборе метода расчёта надо исходить из требований конкретной физической задачи.
Уравнения Ландау и Фоккера-Планка можно рассматривать как одну из двух стартовых точек для данной диссертации. Другой стартовой точкой являются известные работы А.Н.Тихонова, А.А.Самарского, Ю.П.Попова, А.П.Фаворского и др., по консервативным и полностью консервативным разностным схемам (ПКРС) для уравнений математической физики.
Большую роль в развитии численных методов и применении их к задачам физики плазмы сыграли работы Ю.Н. Днестровского, Д.П. Костомарова, J.Killeen, A.A.Mirin и других авторов. Первый шаг в построении ПКРС для кинетического уравнения ЛФП был сделан А.В. Бобылевым и В.А. Чуяновым [45]. Эта работа и исследования автора диссертации были обобщены и развиты в последующих работах ([46-50]).
В диссертации обсуждаются основные проблемы, возникающие при построении разностных схем для уравнения ЛФП (2,3). Фактически, в диссертации рассматриваются не только уравнения ЛФП для кулоновского взаимодействия, но и их обобщение на произвольные потенциалы - так называемые уравнения типа Ландау [22, 23]. Обосновывается важность различных форм записи уравнения ЛФП для последующего использования их при конструировании разностной схемы в численном алгоритме. Отмечается роль инвариантов, имеющих непосредственный физический смысл для математической модели. Разработан подход к построению конечно разностных схем в сочетании со сравнительно высокой точностью (второй порядок аппроксимации по пространству скоростей) дал возможность в ряде случаев получить асимптотические решения, проверить аналитические подходы и другие методы моделирования, решить численно
важные прикладные задачи. Основные результаты диссертации опубликованы в работах [22,44,51-91].
В отсутствие источников и стоков уравнение ЛФП описывает процесс релаксации начального распределения к равновесному максвелловскому распределению. Задача релаксации - классическая задача кинетической теории газов и плазмы, решение которой является стартовой точки любой столкновительной модели [92-95]. Для нелинейного уравнения ЛФП эта задача не может быть численно решена корректно без использования полностью консервативных разностных схем. Она также является естественным тестом для любой более сложной численной модели кулоновских столкновений и рассматривается подробно в диссертации. Во второй главе рассматривается релаксация к равновесию начального распределения модельного газа частиц со степенными потенциалами взаимодействия U — а /г'3,1 < (5 < 4 на основе изотропного кинетического уравнения типа Ландау.
В связи с разработкой и широким использованием мощных источников частиц и энергии в последние десятилетия постоянно возрастает интерес к неравновесным состояниям разнообразных физических систем. Источник и сток энергии (частиц) может обеспечиваться ионными пучками, мощным лазерным излучением, током эмиссии, потоками заряженных частиц, выделяемых при реакциях синтеза или деления, и т.п. [96-101]. В работах [102-107] найдены степенные стационарные решения кинетического уравнения Больцмана, описывающие распределение частиц с потоком от источника к стоку. Показано, что при кулоновском взаимодействии распределение с постоянным
потоком энергии является локальным. Найдены точные степенные решения, обращающие в нуль интеграл столкновений Больцмана. Эти решения аналогичны колмогоровским спектрам в инерционном интервале [108,109]. Однако это формальная, методическая общность. Найденные, в том числе и в данной работе, решения устанавливаются в результате прямого взаимодействия (столкновений) частиц. В третьей главе исследуются квазистационарные функции распределения для уравнения ЛФП при наличии источников и стоков для изотропной однокомпонентной плазмы. Проводится сравнение с экспериментальными результатами по облучению полупроводниковой плёнки.
В четвёртой главе подробно исследована задача о релаксации в двухтемпературной электрон-ионной плазме. Хорошо известны классические результаты для этой задачи: приближённая формула Ландау для температур электронов и ионов и формально уточняющая этот результат так называемая формула Спитцера [НО]. Возможно ли уточнить эти формулы путём численных расчётов? Ответ оказывается положительным, что свидетельствует о высоком качестве построенных в диссертации разностных схем. Задача о релаксации температуры электронов и ионов является тестом для более сложных задач.
Пятая глава посвящена численным методам для неизотропного уравнения ЛФП с аксиальной симметрией (кулоновское взаимодействие). Обсуждается важный вопрос о приближенном вычислении потенциалов Розенблюта. Насколько законно это приближение? Частичный ответ на этот вопрос даёт сравнение с результатами, полученными методом Монте-Карло для системы частиц с кулоновским взаимодействием. Этот метод, предложенный в 2000 году
в [43] является, по-видимому, самым быстрым из существующих методов решения трёхмерных задач для уравнения ЛФП. Для сравнения была выбрана классическая задача о релаксации продольной и поперечной температур (иногда её называют задачей о релаксации импульса).
В шестой главе описаны основные для данной диссертации приложения развитых в работе численных методов к конкретным задачам столкновительной плазмы, находящейся во внепшем электромагнитном поле. Рассматривается ускорение и нагрев частиц в магнитных ловушках: открытых (пробочные, или зеркальные ловушки и магнитосфера планет) и замкнутых (токамаках). Глава 6 является самой большой по объёму, причём основное место в ней уделено приложениям, а используемые численные методы уже были описаны в предыдущих главах.
В предыдущих главах уравнение ЛФП для функции распределения f(v,fi,t) рассматривались в полном пространстве скоростей 0<и<оо,—1
Поскольку в простом пробкотроне даже при больших пробочных отношениях нельзя добиться больших коэффициентов мощности, в разное время был предложен ряд усовершенствованных вариантов простого пробкотрона: центробежные ловушки, многопробочные и другие соленоидальные системы. Открытая магнитная ловушка является полезным инструментом для изучения общефизических свойств плазмы. Кроме того, на их основе можно
создать высокопоточный генератор "термоядерных" (14 МэВ) нейтронов для материаловедческих и других исследований [115-120].
В диссертации задача рассматривается в постановке [16,18,113,114], когда для магнитного поля берётся приближение магнитной ямы. Тогда при аксиальной симметрии задачи магнитное поле входит лишь в граничные условия для функции распределения в пространстве скоростей. Характерным параметром задачи является пробочное отношение R — Вт/Во, где Вт и Во величина магнитного поля на конце и в центре ловушки. В качестве источника частиц выбирается моноэнергетический пучок, а стоки моделируются уходом через пробки частиц с большой продольной скоростью. Здесь важно отметить важность используемых ПКРС, поскольку при длительных расчётах не накапливаются ошибки, искажающие решение. Эта простая модель может быть полезна как базовая и в задачах, связанных с магнитосферной плазмой, где есть пробочные конфигурации, и в применении к токамакам (например, для расчёта влияния захваченных частиц).
Поскольку статическое электрическое поле и МГД волны могут генерироваться и сосуществовать одновременно как в космической, так и в лабораторной плазме, их совместной действие остаётся открытой проблемой и объектом изучения, благодаря широкому приложению . Задача ускорения потока частиц и появления убегания электронов за счет воздействия МГД волны и прямого электрического поля, или за счет воздействия нескольких пакетов МГД волн изучалась многими авторами (см., например, [121-132]). Понимание структуры функции распределения - важный вопрос как с точки зрения объяснения
наблюдаемых явлений, так и понимания основных плазменных процессов. В данной работе нагрев и ускорение электронов при воздействии статического и высокочастотного электрических полей рассматривается во втором разделе главы 6. Вводится основная для этой главы математическая модель: анизотропное (осесимметричное) уравнение ЛФП для электронной функции распределения f(v,n,t) (и =| v |, /х = cos0) с учётом постоянного внешнего электрического поля Е и квазилинейной диффузии
ft/= '[/,/] + Df + Ef. Операторы Df и Ef моделируют поглощение ВЧ поля и влияние постоянного электрического поля, соответственно. Принимая за выделенное направление электрическое поле Е: Е\\ — Ez, в цилиндрических координатах имеем
/ = 7V, 7 = #, Ес = ^, 1^ = 1.5/7.
Здесь Ес - так называемое поле Драйсера, vc - соответствующая ему критическая скорость, Vth - тепловая скорость, te - время электрон-электронных соударений. Величина 7 ^С 1, например, для параметров солнечных факелов (flares) она равна 7 = 0.1-0.2 [133-135].
Уравнение ЛФП является в диссертации базовым уравнением в моделях, связанных с взаимодействием волна-частица, которые используют квазилинейную теорию (формально это означает добавление к интегралу столкновений диффузионного оператора). Диффузионный оператор и уравнение ЛФП объединяется в алгоритме расчета естественным образом. Зачастую оператор столкновений играет роль регуляризирующей добавки. Механизм поглощения
волн за счёт электронного затухания Ландау рассматривается в рамках стандартной квазилинейной теории взаимодействия волна-частица [136-138]. Проведена серия вычислений для максвелловских начальных условий при различных параметрах задачи. Детально изучается структура функции распределения, ток, индуцированный волнами, качественные и количественные аспекты явления убегания, влияние квазилинейной диффузии на убегание электронов, а также возможности использования в рамках единого алгоритма двух систем координат: сферической и цилиндрической. Рассмотренная задача, которая носит в значительной степени методический характер, далее применяется к конкретной физической задаче о нагреве плазмы альфвеновскими волнами и создании токов увлечения в токамаке. Эта задача имеет долгую историю в литературе [139-145]. Исследование генерации тока увлечения в диссертации проведено на базе дрейфово-кинетического уравнения с интегралом ЛФП. Эта работа была первоначально сделана для строящегося в СФТИ токамака. В дальнейшем продолжение этой работы связано с токамаком TCA/BR. Эффекты влияния хвостовых электронов важны для при экспериментальных измерениях дисперсии альфвеновских волн, которые, в свою очередь, определяют так называемую альфвеновскую диагностику эффективной массы атомов или q-profile в больших токамаках.
Известно, что высыпание электронов во время суббури связано с взаимодействием волна-частица в магнитосферной экваториальной плоскости [146,147]. Эти волны могут генерироваться в магнитосфере Земли благодаря мазер-эффекту [148]. Динамика процесса высыпания обычно оценивается при
условии, что фазовая скорость свистов постоянна во времени. В следующем разделе изучается ускорение электронов гидродинамическими волнами с переменной фазовой скоростью, при этом типичные параметры динамики хоров взяты из наблюдаемых данных [149,150].
В последнем разделе исследуется аномальное высыпание электронов около земной авроральной зоны, индуцированное волнами [151-160]. Математическая модель, используемая нами, базируется на модели достаточно проста. Однако наблюдаемые данные, взятые из спутниковых измерений, могут быть объяснены с помощью результатов численного моделирования и имеют хорошее качественное согласие. Результаты моделирования могут быть использованы для более сложной физической модели.
Основные результаты, полученные в диссертации состоят в следующем:
Для нелинейного многокомпонентного одномерного и двумерного в пространстве скоростей уравнения типа Ландау и Фоккера-Планка впервые построен класс полностью консервативных разностных схем, наиболее адекватно отражающих симметрию точных уравнений и выполняющих основные законы сохранения.
На базе построенных полностью консервативных разностных схем разработаны эффективные вычислительные алгоритмы, позволяющие вести расчеты длительное время без накопления ошибок.
На основе разработанных математических моделей и комплексов программ
впервые и систематически исследованы:
- задача нелинейной релаксации двухтемпературной электрон-ионной плазмы
в сильно неравновесном случае;
- релаксация функции распределения и формирование високоенергетичного
хвоста для газов частиц с дальне-действующими потенциалами взаимодействия;
формирование квазистационарных неравновесных распределений электронов плазмы, подверженной облучению пучками быстрых ионов или электромагнитного излучения;
нагрев и ускорение ионов альфвеновскими волнами, нейтральными пучками и генерация токов увлечения в магнитных ловушках;
ускорение и высыпание электронов в авроральной зоне Земли, а также формирование протяженных плато в распределении частиц для магнитосферной плазмы.
Теоретическая и практическая ценность результатов диссертации заключается в разработке эффективных численных методов решения пространственно однородного нелинейного кинетического уравнения ЛФП и изучении общих закономерностей нелинейной динамики слабо столкновительной плазмы.
Исследованы конкретные прикладные физические задачи. Полученные результаты, разработанные математические модели и численные методы могут быть инкорпорированы в более сложные физические модели, а также служить для них надёжными тестами.
Работа выполнялась в рамках научных планов Института прикладной математики им. М.В.Келдыша РАН, проектов Государственного фонда фундаментальных исследований Российской Федерации.
Создано достаточно обширное программное хозяйство, отлаженное и
оттестированное, которое может быть рекомендовано в научных учреждениях, проводящих исследования в области разреженного газа и плазмы. Научные положения диссертации и разработанные на их основе методики, алгоритмы и программные комплексы использовались для совместных исследований лабораторной, космической и полупроводниковой плазмы в следующих научных организациях: РНЦ "Курчатовский институт", ИЯФ им. Г.И.Будкера СО РАН, ННЦ ХФТИ НАН Украины, Сухумский физико-технический институт (Сухуми, СССР), Институт физики плазмы Макса Планка (Гархинг, Германия), Университеты Рио-де-Жанейро, Кампинаса, Сан Пауло (Бразилия), Институт физики А. Вольта ( Павия, Италия).
Всего по теме диссертации опубликовано около 70 научных работ. Результаты исследований, приведенных в диссертационной работе, были представлены и обсуждались на Всесоюзных, Всероссийских и Международных конференциях, в том числе:
о XXV International Congress on Rarefied Gas Dynamics, St.Petersburg, July, 2006 о Международный конгресс по физике плазмы - ICCP06, Киев, май, 2006 о Сессия Совета по нелинейной динамике РАН, Москва, декабрь, 2004,2005 о International Conference on Conservation Laws and Kinetic Theory, Shanghai, China, July, 2005
о Международная Конференция MCC-04, Трансформация волн, когерентные структуры и турбулентность, Москва, ноябрь, 2004
о International Conference on Plasma Physics and Controlled Fusion, Alushta, Crimea, Ukraine, September, 2004
о Межгосударственное рабочее совещание "Плазменная электроника и новые методы ускорения", Харьков, Украина, 2000, 2003
о International Conference on Plasma Physics and Controlled Fusion, St.Petersburg, July, 2004
о Всесоюзные и Всероссийские конференции по физике плазмы, Звенигород, 1986, 1988, 1989, 2002, 2003
о International Topical Conference on Plasma Physics, Faro, Portugal, September, 2001
о International Topical Conference on Frontiers in Plasma Physics, ICTP, Trieste, Italy, 1997, 2000, 2006
о VIII Latin American Workshop on Plasma Physics, Tandil, Argentina, 1998
о International Conference on Computational Physics, Granada, Spain, 1998
о II Panamerican Workshop of Computational and Applied Mathematics, Gramado, Brazil, 1997
о V Latin American Workshop on Non- Linear Phenomena, Canela, RS, Brazil, 1997
о International Workshop on Astrophysics and Space Plasma, Guaruja, Brazil, June, 1995.
о VI Latin American Workshop and International Conference on Plasma Physics and Controlled Fusion, Foz do Iguacu, Brazil, 1994
о II Symposium on Plasma Dynamics. Theory and Applications. Trieste, Italy, July, 1992
о XVII European Conference on Plasma Physics and Controlled Fusion, Amster-
dam, Netherlands, 1990
о Всесоюзное рабочее совещание по открытым ловушкам, Сухуми, СССР, 1990
о XVI European Conference on Plasma Physics and Controlled Fusion, Venice, Italy, 1989
о International Symposium on Matematical Models, Analytical and Numerical Methods in Transport Theory, Minsk, USSR, 1986
Диссертация состоит из введения, шести глав, заключения, списка литературы.
Формы записи интеграла столкновений Ландау- Фоккера-Планка и полностью консервативная разностная схема
Из кинетического уравнения, описывающего эволюцию функции распределения частиц, вытекают как математические следствия, основные законы сохранения - числа частиц и их энергии. Разностная схема, вообще говоря, таким свойством не обладает. Это приводит к тому, что дискретные аналоги указанных законов в процессе расчета нарушаются (схема неконсервативна), в результате чего решение может существенно исказиться не только количественно, но и качественно. Дисбалансныен члены,порождающие искажения, формально имеют тот же порядок, что и аппроксимация схемы. Однако, они зависят от производных решения и поэтому могут значительными на реальных сетках, если решение изменяется сильно. Особое значение эти вопросы приобретают также для задач, где присутствуют частицы, резко отличающиеся по массе, что порождает процессы с разными масштабами времени.
Таким образом, целью является построение полностью консервативных схем, для которых из разностного уравнения для сеточной функции распределения вытекали бы как алгебраические следствия соотношения, выражающие выполнение законов сохранения частиц и энергии для дискретной модели.
Следует заметить,что для нелинейного кинетического уравнения Ландау-Фоккера-Планка имеется нетривиальная ситуация, когда для одного уравнения справедливы несколько, по крайней мере два, закона сохранения (плотность числа частиц и энергии).
Если разностная схема обладает только приближенным аналогом законов сохранения, то это легко ведет к накоплению ошибок, что особенно опасно при расчете нестационарных задач. Схемы, которые обладают более чем одним законом сохранения мы называем полностью консервативными. Общего метода построения полность консервативных схем не существует. Идеальной является ситуация, когда можно обеспечить выполнение законов сохранения лишь аппроксимацией на пространственной сетке. Таким образом, не ухудшая порядок аппроксимации по пространству скоростей, мы имеем свободу при выборе порядка аппроксимации по времени. Далее мы дадим общую идею подхода к конструированию ПКРС.
Важной особенностью рассматриваемых здесь уравнений является то, что они - интегро-дифференциальные. Поэтому, прежде всего, удобно на простейших примерах выделить два различных механизма консервативности, которые присутствуют в уравнениях типа ЛФП одновременно.
Разностная схема удовлетворяет численному аналогу закона сохранения плотности (числа частиц) (10). Действительно, умножая обе части разностного уравнения (9) на /ц+і/2 после суммирования по всем г — 1, ...М + 1 получим равенство к fc-i An Atf-/: = / ЛІ+І/2 Т f-f Г г=1 fc ах OX Л/2 М+1/2 которое равно нулю если предположить, что выполнены дискретные граничные условия: поток частиц через границы равен нулю. То есть мы получим разностный аналог закона сохранения.
В диссертации подобный оператор диффузионного типа в двумерном случае входит составной частью в модель нагрева и ускорения плазмы электромагнитными волнами. Нагрев моделируется в приближении квазилинейной диффузии на основе диффузионного оператора в двумерном пространстве скоростей.
2. Уравнение типа Колмогорова
К этому типу уравнения, в частности, относится односкоростное кинетическое уравнение (для нейтронов). = JdyK(x,y)[f(y)-f(x)]t K(x,y) = K(ytx), 0 х 1. (11) Закон сохранения массы (числа частиц) dn d Г1 , „ ч -r = I dxf(x, t) = 0. dt dt Jo Если известно точное аналитическое выражение к(х) = / dyK(x,y), Jo то это уравнение можно записать в виде = J dyK(x, y)f(y) - k(x)f(x). (12)
Естественно было бы строить разностную схему второго порядка точности по х, заменяя интеграл по формуле трапеции, а второе слагаемое - его точным выражением в уравнении (12). Однако, как нетрудно видеть, в результате получим неконсервативную схему, которая не допускает записи в форме аналогичной (11). Здесь опасность неконсервативности заключается в том, что теряется равновесное решение (/ = const), и поэтому для достаточно больших временных интервалов схема всегда дает неверные результаты (фиктивный источник или сток частиц): м Е І=І м An Л Л"-/Г1. _ #0. hi / J hm+l/2K{xii Zm)Jm kiji i+l/2 = _, tli+l/2 г г- т m=l 1=1
Для построения консервативной схемы второго порядка точности на основе уравнения нужно вместо точного выражения для функции к (х), брать приближенное выражение, полученное аппроксимацией интеграла по формуле трапеции. Та же схема получится, очевидно, автоматически, если использовать формулу трапеции, взяв за основу менее точное исходное уравнение (11).
Полностью консервативная разностная схема для изотропного уравнения типа Ландау
Для разностного уравнения в первой точке і = 2 из условия J3/2 = 0 следует т т «2 = -jrA3, /ЗІ = -jT-Ba, 72 = 0. Это означает, что функции в начальных точках связаны между собой соотношением в любой момент времени. Далее применяется стандартный метод прогонки [26,27]. При аппроксимации граничного условия точной задачи (второго рода), мы имеем определённую свободу. Если не обращать внимания на выполнение законов сохранения, а следить, скажем, за высокой точностью аппроксимации, мы можем при расчете получить постоянно накапливающуюся погрешность. Проверим, что значит формально полученная аппроксимация для функции распределения. Выражение для аппроксимации функции в первой точке можно получить исходя из конечного решения f(v,t) — Дг при t — 00. Предположим, что в дискретном случае имеем функцию максвелловского типа. Заметим, что
Поскольку V2 — vi + h — h (для равномерной сетки), то пренебрегая членами высших порядков по сравнению с О(h2) получим искомую аппроксимацию из разложения экспоненты. Заметим, что и в дальнейшем при конструировании разностных схем для других кинетических уравнений, формальное требование полной консервативности приводит к правильной и даже более естественной аппроксимации точных граничных условий для функции распределения.
Поскольку, как правило, выбраны схемы неявными по времени,то шаг по времени т, как правило, определяется лишь требуемой точностью расчета и характером распределения: т-1 max /fc).
Дивергентные свойства схемы не меняются, если использовать другую аппроксимацию по времени, то есть можно перейти к явной схеме или схеме порядка 0(т2 + h2). Нелинейные разностные уравненияы решаются на каждом шаге по времени методом итераций (аналогично квазилинейному уравнению теплопроводности), а на каждой итерации методом прогонки. Прогонка корректна при условии hi Ei/NiVi, hi = Vi — Vi-\.
Проверка порядка аппроксимации выполнялась обычным образом для численного решения нелинейных задач: при измельчении шага разностной сетки проверялась скорость сходимости к решению.
Схема немонотонна и может быть сделана монотонной, если ограничиться первым порядком аппроксимации по пространству. Это обстоятельство может быть полезным при решениях задач с очень большими интервалами по скорости не требующих высокой точности результата.
При решении кинетического уравнения типа Ландау (Фоккера-Планка) желательно использовать неравномерные сетки, как по времени, так и по скорости. Строго говоря, чтобы учесть экспоненциально быстрое спадание функций при больших скоростях, надо шаг по скорости в области больших энергий уменьшать, хотя каждый раз надо исходить из особенностей поставленной задачи. Приблизительная оценка для шага по скорости, которая в конкретных задачах может варьироваться, следующая
Для степенных потенциалов взаимодействия U — a/rs, где 1 s 4,симметричное ядро Q(v, w) может быть представлено как a(v,w) - [(n+A)vw-(v2+w2)], b(v,w) = [(n+A)vw+(v2+w2)}, n — (s-A)/s. Отрицательные значения n соответствуют мягким потенциалам (1 s 4). Для заряженных частиц s = \(п = — 3).
Для иллюстрации работы разностной схемы начальное распределение, имеющее вид моноэнергетического пучка f(v,0) — 6(v — l)/t 2. Такой выбор функции распределения соответствует нормировке (16).
Разностная аппроксимация этого распределения имеет вид /Ко) = { 2 / (vi+i — Vi-i) если vt = 1, О otherwise.
То есть, функция отлична от нуля в одной точке сетки, в которой она имеет немалое значение. Подобная аппроксимация моделирует разностную дельта-функцию и дает равенство единице как плотности, так и энергии частиц в разностном случае. Ограничимся представлением проведенных расчетов только с данным типом начального распределения.
Было проведено сравнение результата расчёта задачи Коши по данной схеме для уравнения в форме Ландау (для заряженных частиц) с аналогичным расчетом по уравнению Фоккера-Планка. Небольшая численная разница обусловлена, главным образом, тем, что аппроксимируется разные уравнения и разные граничные условия в нуле. Вид функции распределения при нулевой скорости для уравнения ФП более пологий в сравнении с расчётом для уравнения Ландау. Проведем оценку точности сохранения законов для построенной схемы, а затем перейдем к анализу физической задачи. Ошибка в сохранении плотности зависит от величины граничного значения vmax = L, которое можно оценить по распределению Максвелла. Если граничное условие f(L,t) — 0 близко по величине к машинному нулю, то число частиц сохраняется с машинной точностью (случайная ошибка).
Эволюция функции распределения при наличии внешней накачки
Для численного решения задачи мы заменяем бесконечный интервал в пространстве скоростей на конечный отрезок [0,г;тоах], который выбирается так, чтобы учесть высокоэнергетичные частицы и граничное условие для функции распределения будет по-прежнему f(vmax,t) = 0. Источники и стоки берутся, главным образом, в виде 5-функции: S± /± S(v — v±)/v2. Если интенсивности равны 1+ = 1-- v-2/v+2, тогда энергия, которая приходит извне равна нулю и мы имеем дело с аналогом потока частиц. Однако, плотность частиц в системе убывает (если источник находится при больших скоростях, чем сток), т.е., в данной ситуации реализуется аналог постоянного потока энергии в импульсном пространстве, но с несохранением плотности частиц в системе. Так как мы имеем дело с заряженными частицами, то при уменьшении плотности электронов в некоторой области из-за возникшего электрического поля туда начнут поступать тепловые электроны из соседней области. В рамках рассматриваемой модели, однородной по пространству, это можно учесть введением еще одного источника с такой интенсивностью, чтобы не происходило уменьшение числа частиц, а значит не возникало электрическое поле. Таким образом, можно сформулировать согласованную модель с двумя источниками с интенсивностями и одним стоком с интенсивностью , в рамках которой не будет происходить изменения энергии и плотности частиц. Для этого интенсивности источников и стоков должны удовлетворять двум соотношениям: Ith — I- + 1+ — 0 и hhv]h — I-v2_ + I+v+ = О Откуда получаются интенсивности источников, выраженные через интенсивность стока /_ Х- / v+-vth
Кроме того, мы рассматривали источники (стоки), распределенные экспоненциально S± 7±ехр[—b(v — v±) ]. (51) для того, чтобы исследовать зависимость неравновесных функций распределения от формы источника. Следует заметить, что, строго говоря, введение отрицательных источников (стоков), не зависящих от функции распределения может быть весьма проблематично. В таком случае всегда можно создать такие начальные условия, что решение станет отрицательным в области стока после некоторого промежутка времени. Поэтому зачастую стоки моделируются членами, отвечающими поглощению частиц (пропорциональными искомой функции распределения).
В дискретном случае функция 5(v — vi) отлична от нуля лишь при v — v\ как и в предыдущей главе. Начальное распределение выбирается максвелловского типа или типа 5-функции. Стоит отметить, что результаты моделирования практически не зависят от вида начальной функции распределения, разве что в самый начальный момент.
Сначала рассмотрим эволюцию функции распределения для различных показателей /?і,2,з = 1)2,4 под действием лишь источников, находящихся в высокоэнергетичной части распределения. Формирование неравновесного распределения может быть разделено на три части. В течение первого короткого периода система помнит начальные условия. Этот период не сильно отличается для разных показателей /3 и равен приблизительно t 1 как и в задаче о релаксации без учета источников. В течение второй стадии происходит формирование основной части распределения. Продолжительность этой стадии существенным образом зависит от местоположения источника г + и не зависит от его интенсивности. Функция распределения принимает форму плато или имеет слабо спадающую зависимость между источником и холодной областью, в зависимости от величины интенсивности источника. Установление стационарного распределения заканчивается формированием хвоста распределения. Эта стадия существенно зависит от показателя
Численный расчёт процесса релаксации ионов и электронов для начальных условий далёких от равновесия
Далее, воспользовавшись соотношением v\ = vl+1 — 2vk+i/2hk+i, из уравнения (85) получим изменение энергии частиц на шаге по времени
Если An 0, то U% _i — О» Ц = 0 и остается проверить в последнем уравнении чему равно второе слагаемое:
По существу, мы повторяем в дискретном случае (следуем) вывод законов сохранения из уравнения - интегрирование по частям, равенство нулю подстановок при v — 0 и v = оо, учет симметрии относительно перестановки индексов по сортам частиц. Таким образом, получаем V2 АЕ = — У2 mav2KAna и па пропорционально f _v На практике VK = L выбирается достаточно далеким, чтобы /_j « / было близким к машинному нулю для самого далекого хвоста распределения. Тогда при счете ошибки сохранения числа частиц и энергии при явной разностной схеме пропорциональны случайным ошибкам и ошибкам округления. Если же схема неявная, то нужно проводить итерирование и тогда сохранение энергии будет зависеть от точности итераций, (см. главу 2). Равенство нулю функции распределения при VK — L делает однозначным выбор начального значения момента Р% — 0. Ниже приведенные оценочные расчеты точности выполнения равенств в зависимости от значения L (в двумерном случае). Постановкой специальных граничных условий на хвосте VK = L в дискретном случае можно добиться точного выполнения законов сохранения (записав, например, особым образом уравнение в граничной точке К), однако на практике простое условие / = 0, аппроксимирующее lim oo/" = 0, вполне хорошее.
Неявная схема для многокомпонентной плазмы особенно важна из-за долгого времени расчета и различных масштабов характерных времен электронов и ионов, когда не надо обращать внимание на устойчивость и соотношения шагов по скорости и времени.
Перепишем уравнение разностное в форме, используемой в численных расчетах. Введем следующие обозначения к к i=l i=l Щ = — V Y, Мб«Д-і + 0-bQ vkhk), (88) «fc-l/2ftfc V Gak = V H Мб«Д - O.5Q0NZvkhk+l), (89) Vk+l/2flk+l где pP определено в (30). Тогда разностное перепишем в стандартной форме ft = 4г [Щ+г!ш - (Щ + Gak) Ла + GUfti) , к = 2,3,-К-1
Надо добавить граничное /а(г л-) = 0 и начальное условие /Q = fa(vk,0). Схема построена на симметричном шаблоне и имеет порядок аппроксимайии 0(т + h2) для равномерной сетки. Ее дивергентные свойства не меняются, если использовать другую аппроксимацию производной по времени, поскольку вся забота была как раз,чтобы не использовать временную производную. То есть все сохраняется,если перейти к явной схеме или ко второму порядку аппроксимации по времени 0(т2 + h2). Нелинейные разностные уравнения решаются на каждом шаге по времени методом итераций (аналогично квазилинейному уравнению теплопроводности), а на каждой итерации уравнения решаются прогонкой. Прогонка корректна при выполнении следующих условий на шаг по времени -і и по скорости Ш mavk ka0QpN0k Последнее условие является более сложным и отражает немонотонность схемы и быстрое спадание функций распределения на хвостах. С ростом скорости и, следовательно, ростом номера /с шаг hk должен уменьшаться. Эта оценка (можно использовать а, (5 = е, г) может быть менее точной, но более простой: К игТ к Leme 37? Ve
Поскольку при численном расчете нужно помнить об экономии, точности и аппроксимации вычислений шаги по скоростям должны быть разные для разных функций распределения, меняясь плавно. Шаг по времени должен быть достаточно мелким, особенно начале расчёта, исходя из масштаба характерного времени наиболее быстрых частиц (электронов) и особенностей начальных распределений.
