WWW.LIB.KNIGI-X.RU
БЕСПЛАТНАЯ  ИНТЕРНЕТ  БИБЛИОТЕКА - Электронные матриалы
 


Pages:   || 2 |

«МЕТОДЫ ОЦЕНКИ И ПОВЫШЕНИЯ ТОЧНОСТИ РЕШЕНИЯ ЗАДАЧ ФИЗИКИ ПЛАЗМЫ МЕТОДОМ ЧАСТИЦ В ЯЧЕЙКАХ ...»

-- [ Страница 1 ] --

Федеральное государственное бюджетное учреждение наук

и

Институт вычислительной математики и математической геофизики

Сибирского отделения Российской академии наук

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

МЕСЯЦ Екатерина Александровна

МЕТОДЫ ОЦЕНКИ И ПОВЫШЕНИЯ ТОЧНОСТИ

РЕШЕНИЯ ЗАДАЧ ФИЗИКИ ПЛАЗМЫ

МЕТОДОМ ЧАСТИЦ В ЯЧЕЙКАХ

05.13.18 - математическое моделирование,

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

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

В.А. Вшивков Новосибирск - 2014 Оглавление Введение........................................... 5 1 Обзор методов численного решения системы уравнений Власова-Максвелла 10

1.1 Плазма, основные характеристики.......................... 10

1.2 Кинетическое описание бесстолкновительной плазмы............... 12

1.3 Методы решения уравнения Власова, основанные на восстановлении функции распределения f (t, x, v)

1.4 Методы частиц..................................... 15

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

.......................... 18 1.5.1 Уравнения движения модельных частиц, форма модельной частицы.. 18 1.5.2 Сеточные ядра................................. 22 1.5.3 Ядра модельных частиц........................... 23 1.5.4 Общая схема метода частиц в ячейках................... 25 1.5.5 Проблема численных шумов метода частиц в ячейках.......... 26 2 Алгоритм уменьшения счетных эффектов метода частиц в ячейках на примере моделирования распада разрыва плотности ионов в одномерной постановке........................................ 31

2.1 Постановка задачи о распаде разрыва плотности ионов............. 31 2.1.1 Исходные уравнения.............................. 31 2.1.2 Начальные и граничные условия...................... 32

2.2 Схема метода частиц в ячейках.........

–  –  –

Литература......................................... 98 Введение Метод частиц в ячейках, возникший еще в середине прошлого столетия, на сегодняшний день широко применяется при моделировании поведения плазмы. Его также используют и при расчете динамических процессов в других средах (жидкости, газы, сплошные среды и др.). Но именно в решении задач физики плазмы он получил наиболее широкое распространение. Как самостоятельный раздел вычислительная физика плазмы сформировалась одновременно с развитием вычислительной техники во второй половине XX века. Это привело к выделению математического моделирования в отдельное направление исследовательской работы. Ввиду больших возможностей диагностики моделируемых в численных экспериментах кинетических процессов и относительно невысоких экономических затрат на их проведение (по сравнению со строительством установок для проведения реальных экспериментов), численные методы играют все более важную роль в решении задач физики плазмы [16].

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

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

Метод частиц в ячейках относится к группе смешанных алгоритмов, вычисления в которых проводятся в два этапа:

на эйлеровой и на лагранжевой сетке попеременно.

Такие преимущества, как отсутствие аппроксимационной вязкости и возможность осуществления сквозного счета при резких изменениях поведения решения в областях больших градиентов, поставили методы частиц для задач с большими объемными деформациями и неустойчивостями на первое место. Вследствие этого методы частиц нашли широкое применение именно для моделирования таких сложных явлений, как неустойчивости в физике плазмы. Главным преимуществом метода частиц в ячейках по сравнению с эйлеровыми методами является простота реализации и экономичность, так как в методах частиц не вычисляется напрямую функция распределения f (t, x, u). Развитая теория, простота реализации данного метода (в том числе и на параллельных ЭВМ), наглядность получаемых распределений частиц, большое число решенных на его основе задач в разнообразных областях физики плазмы составляют сильные стороны этого метода.

Актуальность работы. Основным недостатком метода частиц в ячейках всегда являлось наличие в решении нефизических эффектов, так называемых «численных шумов». Эта проблема очень сложна, так как причин возникновения шумов много [18]. Влияние различных факторов на решение трудно разделить, особенно при моделировании таких неустойчивых сред, как плазма. Шумовые гармоники взаимодействуют друг с другом, накладываются на гармоники неустойчивых решений, что может приводить к развитию численных неустойчивостей.

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

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

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

• разработать методику выбора оптимальной формы ядра, которая обеспечивает по сравнению со стандартным PIC-ядром уменьшение самосилы,

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

Научная новизна работы заключается в следующем:

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

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

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

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

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

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

Представленные в диссертации исследования проводились в рамках Интеграционных проектов СО РАН №113, №40, Научной школы Годунова - 4292.2008.1 и по проектам, поддержанным Российским фондом фундаментальных исследований (№ 08-01-00615, 11-01-00249, 11-01-00178).

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

На защиту выносятся:

• алгоритм вычитания шумовой добавки,

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

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

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

Апробация работы. Основные положения диссертационной работы докладывались и обсуждались на семинарах «Математическое моделирование больших задач» под руководством д.ф.-м.н. Вшивкова В.А. (ИВМиМГ СО РАН), на семинаре «Математическое и архитектурное обеспечение параллельных вычислений» под руководством д.т.н. Малышкина В.И. (ИВМиМГ СО РАН, декабрь 2009), на семинаре «Задачи механики и математической физики» под руководством д.ф.-м.н. Медведева С.Б. (ИВТ СО РАН, июнь 2013), на объединенном семинаре ИВМиМГ СО РАН и кафедры вычислительной математики НГУ под руководством д.ф.-м.н. Ильина В.П. (ИВМиМГ СО РАН, октябрь 2013), на семинаре «Законы сохранения и инварианты» под руководством д.ф.-м.н. Медведева С.Б. (ИВТ СО РАН, февраль 2014), на Всероссийской конференции по математике и механике, посвященной 130летию ТГУ (2008, Томск), на Международной конференции, посвященной 100-летию со дня рождения С.Л. Соболева (2008, Новосибирск), на Конференции молодых ученых ИВМиМГ СО РАН (2009, 2011, Новосибирск), на XIII Всероссийской молодежной конференции-школе «Современные проблемы математического моделирования. Математическое моделирование, численные методы и комплексы программ» (2009, Дюрсо), на VI Всесибирском конгрессе женщин-математиков (2010, Красноярск), на XV Байкальской Всероссийской конференции «Информационные и математические технологии в науке и управлении» (2010, ИркутскБайкал), на Всероссийской конференции по вычислительной математике КВМ-2011 (Новосибирск), на Российско-монгольской конференции молодых ученых по математическому моделированию, вычислительно-информационным технологиям и управлению (2011, Иркутск (Россия) – Ханх (Монголия)), на XIX Всероссийской конференции «Теоретические основы и конструирование численных алгоритмов и решение задач математической физики» (2012, Дюрсо), на XIII Всероссийской конференции молодых ученых по математическому моделированию и информационным технологиям (2012, Новосибирск), на Международной конференции «Mathematical modeling and computational physics» (2013, Дубна).

Основные результаты опубликованы в 15 работах, из которых 2 в журналах, рекомендованных ВАК [1]- [15].

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

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

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

Глава 1 Обзор методов численного решения системы уравнений Власова-Максвелла

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

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

Состояние плазмы характеризуется большим числом параметров [19]: степень ионизации, плотность n, температура T и др.

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

Плотность плазмы для космических сред изменяется в диапазоне от 106 см 3 для межзвездного пространства до 1023 1024 см 3 и выше в звездах. Пределы изменения плотности лабораторной плазмы от 1013 1014 см в исследованиях по управляемому термоядерному синтезу в токамаках до 1023 1024 см в специальных мишенях для лазерного термоядерного синтеза.

Температура плазмы в зависимости от ее происхождения также изменяется в широких пределах. Плазма с T 105 K считается низкотемпературной, с T 106 K - высокотемпературной. Температура плазмы для УТС должна превышать 108 K. В случае неизотермической плазмы различают температуру электронов Te, ионов Ti и неионизованных атомов Ta, если таковые есть.

Различие масс ионов mi и электронов me ведет к тому, что в плазме существуют процессы разного временного масштаба. В простом случае однородной полностью ионизованной водородной плазмы (ne = ni ), можно выделить четыре характерных времени: e - время релаксации (время установления максвелловского распределения) функции распределения электронов в результате их столкновений между собой, i - аналогичное характерное время для ионов, ei - время релаксации относительного движения электронов и ионов, T - время выравнивания температур электронов и ионов в неизотермической плазме. Cамое «быстрое»

время - e, т.е. релаксация импульса и энергии электронов происходит быстрее, чем релакmi сация ионов при той же температуре. ei e, i. Самым «медленным» процессом me e

–  –  –

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

–  –  –

плазма считается бесстолкновительной.

В зависимости от интересующего нас диапазона параметров плазма описывается математическими моделями разной степени детализации: кинетические модели описание со столкновениями и без, магнитогидродинамическое приближение (МГД), модели термоядерной плазмы и другие модели. Кроме того для разных компонент плазмы могут использоваться разные подходы. Такие модели называют гибридными. Например электроны рассматриваются в приближении МГД, а ионы в кинетическом приближении [20].

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

1.2 Кинетическое описание бесстолкновительной плазмы Свою концепцию описания широкого круга кинетических процессов в полностью ионизованной бесстолкновительной плазме Власов предложил в 1938 году [22,23]. В основе модели Власова лежит кинетическое уравнение без столкновительного члена для функции распределения частиц по скоростям f (t, r, v)

–  –  –

E, B – напряженность электрического и магнитного полей, c – скорость света, j – плотность тока,

- плотность пространственного заряда.

Электрическое и магнитное поля, входящие в уравнение через силу Лоренца, определяются из уравнений Максвелла

–  –  –

В отсутствии магнитного поля система уравнений Максвелла сводится к уравнению Пуассона (1.4) для потенциала, где E =.

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

На данный момент существет целый ряд методов для численного решения системы уравнений Власова-Максвелла [17, 18]:

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

1.3 Методы решения уравнения Власова, основанные на восстановлении функции распределения f (t, x, v) Исторически методы частиц возникли и развивались как альтернатива сеточным алгоритмам. Данная группа алгоритмов, получившая общее название Эйлеровы алгоритмы (Eulerian solvers), объединяет в себе методы вычисления значений функции распределения f (t, x, v) для уравнений Власова-Пуассона или Власова-Максвелла на связанной с областью пространственной сетке. Сюда относятся следующие методы: конечно-разностные методы (finite-difference method), метод конечных объемов (finite volume method) или метод баланса потоков (flux balance method), спектральный метод (spectral method) или метод преобразований, полу-Лагранжевы методы (semi-Lagrangian methods), метод конечных элементов и различные их комбинации и модификации. В данной классификации приводятся разные, в том числе и английские, названия схожих методов, которые встречаются в литературе.

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

Достоинством конечно-разностных методов является развитая теория. Эти методы имеют большую доказательную базу. Некоторые из этих методов для уравнения Власова можно найти в работах [24–26]. Но данные методы не получили широкого распространения как из-за сложности вычислений при большой размерности задачи, так и по ряду других причин. Конечно-разностные методы приводят к рассеиванию или диссипации мелкомасштабных структур, возникающему вследствие схемной вязкости. Конечно-разностные схемы плохо работают на задачах, в которых имеются области больших градиентов. В таком случае приходится сильно измельчать или строить неравномерную сетку, что только усложняет вычисления.

Есть работы, посвященные решению системы уравнения Власова методом конечных элементов [27, 28], методом конечных объемов (метод баланса потоков) [29–32]. Но широкого распространения для задач физики плазмы эти методы не получили.

Метод преобразований или спектральный метод - это использование разложения функции распределения по ортонормированным системам базисных функций, например разложение в ряд Фурье или использование в качестве базисных функций полиномов Эрмита или Чебышева [33–35]. Данный метод также применяется для узкого класса задач, чаще всего одномерных.

С того момента, когда Ченг и Кнорр [36] предложили свою схему расщепления как эффективный метод интегрирования уравнений Власова-Пуассона, был разработан целый ряд методов решения уравнения Власова, основанный на этой схеме [37–40]. В работах [41–43] описан метод, получивший название полу-Лагранжев (semi-Lagrangian method). Вычисления в полу-Лагранжевых методах происходят путем переноса по характеристическим кривым значений функций c предыдущего шага и интерполяции на новом шаге. Sonnendrucker и др. в своих работах [41] не разделяют схемы с расщеплением и полу-Лагранжевы методы. В работе же [44] полу-Лагранжевыми называются лишь методы без расщепления. Общим в них является перенос значений по характеристикам и использование схем интерполяции. Эти методы интересны прежде всего вследствии низкого уровня шумов, присущему этим методам.

Они свободны от таких численных артефактов, как диффузия и диссипация, и могут быть полезны в случае необходимости воспроизведения тонких эффектов, когда например поведение частиц в хвосте функции распределения играет важную роль. Но в двух-, а особенно трехмерных расчетах эти методы очень затратны. Обзор и сравнение разных вариантов данных методов между собой, а также с методом конечных-разностей, методом баланса потоков и спектральным методом можно найти в работая [26, 45, 46].

Стоит упомянуть также метод водяного мешка [47], который является аналогом лагранжевого метода применительно к кинетическому уравнению Власова в фазовом пространстве. В нем рассчитывается эволюция выделенных контуров в фазовой плоскости, плотность заряда и напряженность электрического поля восстанавливаются по этим контурам (Water-bag method). В силу своей сложности данный метод используется как правило только в одномерной постановке.

1.4 Методы частиц

Методы частиц - это особая группа вычислительных алгоритмов, которые объединяет между собой способ дискретизации исходной математической задачи. В этом методе среда (в нашем случае плазма) представляется набором достаточно большого числа модельных частиц. Движение частиц подчиняется законам классической механики в самосогласованном электромагнитном поле, определяемом из уравнений Максвелла. Методы частиц широко применяются при решении задач не только физики плазмы, но и динамики жидкостей и газов, механики сплошной среды и т.д. Но самое широкое применение они получили именно в физике плазмы. Большим плюсом этих методов по сравнению с приведенными выше, является экономичность, т.к. в методах частиц не вычисляется напрямую функция распределения f (t, x, v). Особенно заметным это преимущество становится при переходе от одномерных постановок к задачам большей размерности.

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

В вычислительной физике плазмы методы частиц начали использоваться с конца 1950х годов. Первый одномерный электростатический вариант - модель плоских листов - был предложен О. Бунеманом в 1959 г [48] и развит позже в работах Дж. Доусона [49].

Метод листов относится к лагранжевым алгоритмам, в которых рассчитывается действие каждой частицы друг на друга. Сюда же можно отнести бессеточный метод частиц конечного размера - FSP (gridless Finite Size Particle) [50]. В нем вместо точечных частиц используются облака частиц, имеющие форму функции Гаусса, а вычисление сил, действующих на частицу, происходит в Фурье-пространстве без перехода к пространственной сетке.

Сложность чисто лагранжевых алгоритмов, их еще называют P P -алгоритмами [51] (англ.

particle-particle – частица-частица), равна O(N 2 ), где N - полное число частиц. Такие вычисления возможны только для сравнительно небольшого числа частиц. Например для сильно столкновительных систем, когда число частиц в Дебаевской сфере ND невелико.

Позже, под влиянием результатов группы Ф.Х. Харлоу из Лос-Аламосской лаборатории, разработавших в 1955 г. для задач газовой динамики метод частиц в ячейках [53–55], методы частиц в ячейках были разработаны и для плазмы.

Метод частиц в ячейках - это эйлерово-лагранжев алгоритм (P M -алгоритм, англ.

particle-mesh – частица-сетка). Силы в нем вычисляются гораздо быстрее, чем в лагранжевых алгоритмах, но менее точно. С помощью сеточных значений можно хорошо представить только поля, вариации которых имеют длину волны большую, чем шаг пространственной сетки. Сложность P M -алгоритма равна O(N log M ), где M - число узлов сетки.

Метод частиц в ячейках иногда в литературе называют методом крупных частиц, что ведет к смешению понятий. Метод крупных частиц [56] является модификацией схемы Харлоу и при решении задач физики плазмы применяется к системе уравнений магнитной гидродинамики (МГД), а не к кинетическим уравнениям. Для повышения экономичности схемы авторы данного метода фактически отказались от разбиения среды на частицы, «крупной»

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

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

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

В группу методов под общим названием PIC (particle-in-cell), помимо классического варианта метода частиц в ячейках, рассматриваемого в этой работе, относят также ряд его модификаций с учетом столкновеительных эффектов. Область применения этих методов находится между бесстолкновительными кинетическими моделями и моделями МГД для плазмы. Реализуется это несколькими путями. Например один сорт частиц (чаще это электроны) считается по модели МГД, а другой по кинетической модели, или же плазма описывается МГД-приближением двухжидкостной модели. В некоторых методах рассматривается уравнение Фоккера-Планка для описания столкновений и используется метод Монте-Карло для его решения. Описание разных вариаций этих методов можно найти в [57–61] Необходимо отметить также модификацию метода частиц в ячейках для замагниченной плазмы. Волны в горячей однородно-замагниченной плазме распространяются как малые возмущения на фоне однородного равновесного состояния. Для описания замагниченной плазмы используется гирокинетическое приближение, являющееся промежуточным между МГД и кинетическим подходом. В этом случае необходимо расчитывать низкочастотные процессы, что означает большие времена, при этом использовать малый временной шаг, т.к.

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

Cохраняя полноту кинетического приближения, гирокинетический подход позволяет устранить высокочастотные составляющие движения частиц при помощи осреднения по фазе циклотронного вращения, что дает возможность сохранить учет физических эффектов, вносимых циклотронным движением, но снимает ограничения по величине шага по времени, что существенно для низкочастотных процессов. Модификация метода частиц для этого подхода получила название f -метод [62–65].

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

Математическое обснование справедливости использования методов частиц для численного решения уравнения Власова можно найти например в работе Аресеньева А.А. [66].

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

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

В пространстве моделирования вводится сетка, называемая эйлеровой, в ячейках которой задаются сеточные величины (в нашем случае это F,, j). На эйлеровом этапе конечноразностными методами решаются уравнения Максвелла и находится новое значение сил на сетке. На лагранжевом этапе по полученному значению сил определяют новые положения частиц (узлов лагранжевой сетки) из уравнений движения частиц, которые являются характеристиками кинетического уравнения Власова. Иногда метод частиц сравнивают с методом конечных элементов, базисные функции в котором представляют частицы. Отличие состоит в том, что частицы двигаются по области.

Метод частиц в ячейках сочетает в себе ряд плюсов, которые ставят его в задачах физики плазмы на первое место. Данный метод довольно прост в реализации, менее затратен по сравнению с конечно-разностными методами вычисления функции распределения в трехмерном пространстве или с методами «частица-частица», когда рассчитывается непосредственное взаимодействие всех частиц со всеми, а так же физически нагляден. Существует ряд пакетов, построенных на этом методе и позволяющих проводить моделирование методом частиц для различных плазменных задач как в декартовой, так и в цилиндрической системе координат, в одно-, двух- и трехмерной постановках, с адаптивным изменением масс, с учетом столкновений и т.п., например KARAT [67], Mandor [68], CFHALL [69], VORPAL [70], DEMOCRITUS [71], MAGIC [72], такие коды The UCLA Plasma Simulation Group как OSIRIS [73], QuickPIC [74] и другие, ссылки на которые можно найти в [75].

1.5.1 Уравнения движения модельных частиц, форма модельной частицы

Рассмотрим классический метод частиц в ячейках применительно к системе уравнения Власова-Максвелла [51, 52, 76]. В уравнении Власова (1.1) опустим индекс сорта частиц и рассмотрим уравнение для функции распределения частиц f (t, r, v). Для точечных частиц функция распределения имеет вид (r rj (t))(v vj (t)).

f (t, r, v) = (1.8) j Здесь – дельта-функция Дирака, j, rj, vj – номер, координата и скорость частицы.

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

–  –  –

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

–  –  –

Индекс p здесь и далее - номер модельной частицы (particle). Такой подход позволяет учитывать коллективные процессы (дальнодействующие силы), устраняя при этом столкновительные эффекты, связанные с близким прохождением пары частиц друг относительно друга. Это приводит к сглаживанию сил на близких расстояниях. В англоязычной литературе такие частицы получили название «облака» (cloud-in-cell) из-за того, что они имеют конечный размер и в то же врем свободно проходят сквозь друг друга.

Свойства S :

1) Область определения S конечна.

2) Функция S частицы симметрична и положительна: S(x, y) = S(y, x) 0.

S(x, y) S(x, y) = 3).

x y

4) Условие нормировки: S(x, y)dx = 1.

Обычно берут Sv (v, vp ) = (v vp ), т.е. в этом случае модельная частица - это группа частиц с одинаковой скоростью. Функция Sr (r, rp ) = Sr (rrp ) задает форму, размер частицы и распределение плотности заряда внутри нее. Она называется ядром модельной частицы (shape function). В процессе перемещения форма частицы и ее ориентация относительно осей координат остается неизменной. Таким образом частицы имеют конечную лебегову меру в пространстве координат и сингулярную в пространстве скоростей.

Подставив (1.9), (1.10) в уравнение Власова, получаем

–  –  –

При интегрировании все члены равны 0, кроме интегрирования по тому же направлению, по которому взят момент. Например для x, используя интегрирование по частям, получаем

–  –  –

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

–  –  –

В формуле (1.20) под интегралом E и B стоят непрерывные функции. В методе же частиц в ячейках значения полей вычисляются чаще всего конечно-разностными методами из уравнений Максвелла (1.2-1.5). В результате вычислений имеем сеточные величины Eg и Bg. И для решения уравнений Максвелла нам также нужны сеточные значения плотности g и тока jg. Индекс g - номер узла сетки (grid).

В качестве сеточного значения плотности и тока берем усредненное по небольшому объему пространства вокруг узла сетки g - g (|g | - его лебегова мера)

–  –  –

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

–  –  –

Интерполяционные формулы могут быть выбраны произвольными. Для упрощения формул берут B-сплайн нулевого порядка и в трехмерном пространстве интерполяционные формулы выбирают в виде произведения одномерных S(x, y, z) = S(x)S(y)S(z). Тогда, если объемы интегрирования выбраны согласованно, получаем

–  –  –

1.5.3 Ядра модельных частиц Одним из основных моментов, характеризующих метод частиц в ячейках, является выбор формы ядра модельной частицы S и соответственно сеточного ядра R. Функция ядра характеризует геометрическую форму, размер модельной частицы и распределение плотности заряда внутри нее. В качестве ядра можно выбрать любую функцию с конечным носителем, удовлетворяющую свойствам функции ядра, которая проводит сглаживание сил и устраняет близкие взаимодействия(соударения), не играющие роли в разреженной плазме.

Ниже представлены некоторые одномерные ядра частиц [18, 78]. Первым и простейшим является NGP-ядро (nearest-grid-point) или ZSP (zero-size-particle), в котором частица отдает свой вклад в плотность в ближайший к ней узел сетки.

Ядро NGP SN GP (x) = (x). (1.29)

–  –  –

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

Параболическое ядро S(x)parab = RP IC (x). (1.33)

–  –  –

Важным моментом здесь является переход с Лагранжевой сетки на Эйлерову (1.23), (1.24) и обратно (1.27), (1.28), который осуществляется при помощи функции R (1.21). Отметим, что в общем случае не обязательно использовать одинаковую функцию при переходе от сетки к частицам и от частиц к сетке.

–  –  –

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

Совокупность проявлений вычислительных ошибок метода частиц принято называть «численным шумом» - термин, заимствованный из теории сигналов. Шумовые гармоники накладываются на гармоники решения и при больших временах бывает трудно судить о решении. Источников нефизических шумов в таких сложных алгоритмах, как метод частиц в ячейках, много: наличие пространственной сетки (на которой определяются значения электромагнитных полей, распределения плотности, скорости, тока), разбиение среды на модельные частицы (изменяется потенциал взаимодействия между частицами на масштабе сеточной ячейки), переходы с частиц на сетку и обратно, дискретный шаг по времени. Все эти погрешности приводят к нарушениям в численных моделях законов сохранения энергии, импульса, изменяются дисперсионные свойства, возрастает число столкновений, происходит стохастический саморазогрев ансамбля частиц, который ограничивает допустимый интервал, отслеживаемый в расчетах.

К вопросу количественной оценки уровня нефизических флуктуаций подходили с разных сторон [79,80]. Подробное описание ряда причин шумов и методов их устранения можно найти в работах [51, 52, 81, 82]. На настоящий момент единого подхода к решению этой проблемы нет. Для уменьшения численных шумов увеличивают число модельных частиц, меняют их форму, выбирают оптимальный по времени и пространству шаг, меняют начальное распределение частиц, используют различные модификации метода исходя из специфики поставленной задачи и др. Предлагаются также алгоритмы сглаживания, но они не устраняют сам шум и могут сглаживать физические эффекты, что нежелательно.

• Количество частиц При моделировании бесстолкновительной плазмы каждая модельная частица соответствует десяткам миллионов физических частиц. Обозначим штрихом величины, относящиеся к модельной плазме. Из условий сохранения заряда и массы следует, что e = e, m = m, где = / - отношение плотностей реальной и модельной плазмы. Тогда, как показано в [16, 18], следующие характеристики реальной плазмы сохраняются и в модельной не смотря на сравнительно малое число модельных частиц. Это собственная плазменная частота pe,i = pe,i, дебаевский радиус экранирования (в предположении равенства реальной и модельной тепловой скорости vT = vT ) D = D. Также сохраняется отношение заряда частицы к массе e/m = e /m, а значит динамика частиц не меняется.

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

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

L) и неравенство во многих задачах запас «бесстокновительности» достаточно велик ( L (L - размер области) для модельной плазмы можно сохранить. При выборе числа частиц (или же при выборе размера области) необходимо также учитывать, что число частиц в дебаевской сфере не должно быть маленьким ND 1.

• Адаптивное изменение количества частиц Следует отметить, что проблема шумов не теряет своей остроты при повышении числа модельных частиц в расчетной области даже при расчетах на суперЭВМ, так как с изменением плотности вещества могут появиться ячейки с небольшим количеством частиц, и общий уровень шумов в расчете будет определяться именно этими ячейками. С другой стороны, в наиболее заполненных ячейках при этом может быть частиц в десятки раз больше необходимого уровня, что создает избыточную нагрузку на вычислительную систему, таким образом, повышая время расчета. Для преодоления этих проблемм создаются методы частиц в ячейках с адаптивным изменением числа частиц [71, 83–88]. Кроме того создан метод частиц в ячейках с частицами переменных размеров [89].

• Форма модельной частицы Соответствие модельных частиц их реальным прототипам далеко от полного совпадения не только количественно.

Модельные частицы конечного размера принято интерпретировать как облака частиц. Но форма частиц и распределение признака в них в процессе моделирования остаются неизменными, что может противоречить реальной ситуации. Описанию свойств ядер и того, какие погрешности вносит такой способ дискретизации, посвящено большое число работ различных авторов. Большинство основных выводов и ссылки на первые работы можно найти в следующих монографиях и статьях [18, 51, 52, 77, 80, 81, 90].

В работе [91] приведено сравнение сохранения полной энергии для NGP и PIC ядер, в [92] для параболического и PIC ядер. С увеличением числа частиц в шаблоне и повышением уровня гладкости полная энергия системы сохраняется лучше, флуктуации плотности и поля уменьшаются. В работе [93] получены погрешности аппроксимации признаков с лагранжевой на эйлерову сетку для NGP и PIC ядер. Показано, что при использовании NGP ядра сеточная функция плотности g аппроксимирует исходную плотность (x) со вторым порядком по h, а отклонение g = |g (xg )| обратно пропорционально среднему числу частиц в ячейке.

При использовании PIC ядра g аппроксимирует исходную плотность (x) также со вторым порядком по h, а отклонение g = |g (xg )| обратно пропорционально квадрату среднего числа частиц в ячейке.

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

В [51] на примере показано, что использование схемы PIC/NGP вносит сравнительно неопасные колебания, в то время как NGP/PIC в некоторых может вести к экспоненциальному росту неустойчивости.

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

Для понижения уровня шума, связанного с начальным распределением частиц по скоростям, используется методика «тихого старта» (или «спокойного старта») [81, 95]. Суть ее состоит в представлении Максвелловского распределения с использованием меньшего, чем общее число частиц J, набора скоростей N. Таким образом плазма представляется в начальный момент времени как N пучков. Для этого алгоритма время вычислений с низким уровнем шумов ограничено временем нарастания пучковой неустойчивости.

• Шаг по времени и пространству Неправильно выбранный шаг по времени и пространству также вносит свои ошибки.

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

c h, pe 2. Шаг по пространству: h C0 D, где C0 - константа первого порядка, зависящая от формы ядра частицы. Для PIC C0.

Все эти условия приводят к затратности вычислений, т.к. зачастую нас интересуют процессы масштаба ионных времен pi, а шаг по времени должен составлять доли электронного периода 2pi, что в силу mi me приводит к большому числу шагов по времени. Вследствии этого в расчетах иногда используют уменьшенное отношение масс mi /me (например 100). Таким образом шаг по времени определяют самые мелкомасштабные процессы. Для того, чтобы избежать этих ограничений, разрабатываются неявные схемы метода частиц в ячейках [96–99].

В работе [100] исследована зависимость времени столкновений и разогрева модельной плазмы (рост кинетической энергии), вызванного дискретностью, от формы ядра частицы и соотношений параметров h и D.

Еще одна погрешность, связанная с сеткой - возникновение так называемой «самосилы» (self-force). Авторы [90,101] приводят анализ причин появления этой нефизической силы.

Более подробно этому вопросу посвящена третья глава данной диссертации.

• Адаптивные и неструктурированные сетки Для уменьшения требований на размер шага сетки созданы методы частиц в ячейках на адаптивных [102–104] и на неструктурированных [105] сетках для задач со сложной геометрией. Эти методы обладают своими трудностями, вызванными в том числе изменением размера модельных частиц из-за изменения размера ячеек сетки. Кроме того эти методы значительно усложняют вычислительный алгоритм метода частиц.

• Система уравнений Максвелла и вычисление тока

–  –  –

то закон сохранения заряда d + divh j = 0 (1.36) dt в разностном виде может не выполняться. Для исправления этой ошибки есть два пути.

Первый - корректировка E с учетом выполнения уравнения Пуассона [18, 106, 107], что довольно затратно, т.к. требуется на каждом шаге решать эллиптическое уравнение. Иногда не требуют полной сходимости при решении уравнения Пуассона, а делают только несколько приближений для понижения уровня ошибки [108, 109]. Второй - вычисление j по частицам с учетом уравнения (1.36) [78, 110–112]. Этот способ в общем случае не сохраняет импульс системы.

Кроме всего упомянутого для уменьшения численных шумов используются такие техники, как фильтры Фурье для обрезания шумовых гармоник [95], разрабатывают модификации алгоритма, основанные на особенностях поставленной задачи, например вышеупомянутый f -метод [63] или метод с периодическим восстановлением функции распределения, предложенный в работе [113]. Все эти методы, используемые по отдельности или в сочетании друг с другом, имеют свои плюсы, так же как и свою область применимости. Но все они ведут к значительному усложнению алгоритма.

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

Глава 2

–  –  –

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

2.1.1 Исходные уравнения

Для задачи о распаде разрыва плотности ионов в бесстолкновительной полностью ионизованной неизотермической плазме, когда температура электронов значительно превышает Ti ), выбрана комбинированная модель [114,115], в которой рассчитемпературу ионов (Te тывается только движение ионной компоненты плазмы (без магнитной индукции), а плотность электронов описывается распределением Больцмана

–  –  –

Здесь v – скорость ионов, f (t, x, u) – функция распределения ионов, n0 – невозмущенная плотность плазмы, Te – температура электронов, e – заряд электрона, E(t, x) – напряженность электрического поля (t, x) – потенциал электрического поля, n(t, x) = f (t, x, u)du – плотность ионов, v(t, x) = 1n(t, x) (u)uf (t, x, u)du – средняя скорость ионов.

2.1.2 Начальные и граничные условия

–  –  –

Здесь C - значение плотности в области слева от разрыва, b - коэффициент сглаживания, точка разрыва находится в центре области x0 = L/2. Для простоты далее везде сохранено название распад разрыва.

Начальные условия для потенциала

–  –  –

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

–  –  –

Значения i и ni задаются в узлах сетки, а Ei+1/2 в центрах ячеек.

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

–  –  –

Эйлеров этап. На эйлеровом этапе решается уравнение Пуассона для потенциала с найденными значениями плотности nk+1 и вычисляется напряженность электрического поля.

i Решение нелинейного уравнения

–  –  –

В работе [116] показана монотонная сходимость последовательности k+1,s к решению уравнения (2.18), если k+1,1 k+1. А в работе [117] показано, что k+1,1 k+1 при любом начальном приближении k+1,0.

В конечных разностях уравнение (2.19) будет иметь вид

–  –  –

За границами области берем E1/2 = E3/2, En+1/2 = En1/2.

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

В силу простоты реализации в качестве такой схемы для уравнения Власова была выбрана схема Лакса-Вендроффа [118]. Порядок аппроксимации этой схемы второй. Уравнение Пуассона решалось как в п.2.1.3.

На [0, L] [umin, umax ] вводится сетка с шагом h по x и hu по u. Шаг по времени время изменяется от 0 до tm ).

Индексом i будем обозначать узлы сетки по x (i = 1,..., n), j – узлы сетки по u (j = 1,..., km).

Обозначим F = f, G = uf, H = Ef. Тогда уравнение Власова имеет вид

–  –  –

Для того, чтобы показать, что метод частиц в ячейках при соответствующих параметрах способен давать хороший результат, проведем сравнение метода частиц с решением по схеме Лакса-Вендроффа с начальной температурой v = 0.05 (2.7) (разброс частиц по скоростям или дисперсия функции распределения).

2.3.1 Сравнение метода частиц и конечно-разностного метода

На рисунках 2.2 и 2.3 приведены результаты расчетов по методу частиц в ячейках и по схеме Лакса-Вендроффа со следующими параметрами:

общие: v = 0.05, tm = 0.5, параметры метода частиц: шаг по времени = 0.001, на рисунке 2.2 сетка по n = 200, число частиц m = 4000000, на рисунке 2.3 сетка по n = 100, число частиц m = 5000.

параметры схемы Лакса-Вендроффа: [0, 4] [0.8, 0.8], сетка 1600 2000.

Средняя скорость в начальный момент нулевая.

Рис. 2.2: Плотностьи средняя скорость ионов в моменты времени t = 0 и t = 0.5, полученные по схеме Лакса-Вендроффа и по методу частиц с параметрами n = 200, m = 4000000.

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

Рис. 2.3: Плотность и средняя скорость ионов в момент времени t = 0.5, полученные по схеме Лакса-Вендроффа и по методу частиц с параметрами n = 100, m = 5000.

2.4 Алгоритм уменьшения счетных эффектов метода частиц в ячейках (алгоритм вычитания шумовой добавки) Нефизические колебания в скоростях приводят к тому, что средние скорости частиц, вычисленные по формуле (2.12), являются негладкими функциями. И хотя амплитуда наблюдаемых колебаний с увеличением числа частиц уменьшается, эти колебания приводят к аналогичным колебаниям электрического поля E.

Основная идея предлагаемого алгоритма уменьшения шума заключается в следующем:

1. Предположим, что вычисленная на k шаге средняя скорость v k гладкая функция (в начальный момент времени v|t=0 = 0, а значит это выполнено).

2. Сделаем шаг с отключеным электрическим полем (E = 0). Частицы, имея разные скорости, переместятся в новое положение. Если снова посчитать среднюю скорость, то она будет негладкой, хотя на частицы не действуют никакие силы.

3. Обозначим шумовую добавку V, она будет определена ниже.

4. Вычитая V из скорости каждой частицы в данной ячейке на шаге k, получим скорректированное значение скорости частиц.

5. Теперь, если вновь на шаге k посчитать среднюю скорость по формуле (2.12), она будет с шумом. Но тогда после прохождения шага, получаем среднюю скорость уже «без шума» (вернее, уровень шума в средней скорости будет меньше).

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

–  –  –

1. На k шаге средняя скорость частиц является гладкой (на первом шаге это обеспечивает способ задания начальных данных).

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

–  –  –

Назовем желаемой среднюю скорость v k+1, получаемую по какой-либо конечноразностной схеме для (2.22) (v k = v k ).

2.2. Сделаем шаг без действия электрического поля. Посчитаем среднюю скорость и обозначим ее v k+1. Она будет с шумом.

–  –  –

где m - масса одной частицы.

Обозначим координату частиц после сдвига, корректирующего положения частиц через xj = xj + di. Плотность после сдвига будет иметь вид

–  –  –

На рисунке 2.5 показана блок-схема цикла модифицированого метода частиц в ячейках.

Рис. 2.5: Блок-схема цикла модифицированого метода частиц в ячейках

2.6 Схемы, использованные в Алгоритме III В пункте 2.1. алгоритма в ходе расчетов использовались две схемы для уравнений (2.23), (2.24):

–  –  –

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

2.7.1 Зависимость уровня шума от количества частиц

–  –  –

Общие параметры: tm = 0.05, = 0.001.

Параметры метода частиц: общее число частиц m = 10000, число узлов сетки n = 100 (число частиц в ячейке lp = 100).

Параметры конечно-разностной схемы: число узлов сетки по x и u равно n = 1600 и km = 3000 соответственно.

–  –  –

Общие параметры: tm = 0.1, = 0.001.

Параметры метода частиц: число частиц m = 40000, число узлов сетки n = 100 (число частиц в ячейке lp = 400).

Параметры конечно-разностной схемы: число узлов сетки по x и u равно n = 1600 и km = 3000 соответственно.

В этих тестах для уравнений (2.23), (2.24) использовалась противопотоковая схема, время выбрано небольшим tm = 0.1. На рисунке 2.6-2.9 cиним изображены средняя скорость и плотность, полученные стандартным методом частиц в ячейках без вычитания шумовой добавки, зеленым - средняя скорость и плотность, полученные по методу частиц в ячейках с Алгоритмом III вычитания шумовой добавки. Из графиков видно, что решение, полученное по методу частиц и методу частиц с вычитанием шумовой добавки ведет себя так же, как и решение, полученное по схеме Лакса-Вендроффа для уравнения Власова. Видно также, что с вычитанием шумовой добавки графики более гладкие, особенно это заметно на тесте с меньшим числом частиц.

Сравнение с решением по схеме Лакса-Вендроффа проводится как с точным решением.

Но аналитически как должно вести себя решение точно известно только в областях A и B, Рис. 2.6: Тест 1: плотность и средняя скорость, вычисленные по схеме Лакса-Вендроффа, по обычному (PIC без вычитания) и модифицированному (PIC с вычитанием шума) методу частиц. Число частиц в ячейке lp = 100.

Рис. 2.7: Тест 1: плотность и средняя скорость, область A. Число частиц в ячейке lp = 100.

куда еще не дошло возмущение. Обозначим эти области A = [0, 1], B = [3, 4] и C = A B.

Средняя скорость в этих областях должна быть равна vA = vB = 0, плотность nA = 2, nB = 1.

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

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

–  –  –

Рис. 2.8: Тест 2:, вычисленные по схеме Лакса-Вендроффа, по обычному (PIC без вычитания) и модифицированному (PIC с вычитанием шума) методу частиц. Число частиц в ячейке lp = 400.

Рис. 2.9: Тест 2: плотность и средняя скорость, область A. Число частиц в ячейке lp = 400.

В таблице 2.1 индексом pic помечены величины, относящиеся к опо обычному (PIC без вычитания) методу частиц в ячейках, индексом new - относящиеся к модифицированному (PIC с вычитанием шума) методу (метод частиц в ячейках с Алгоритмом вычитания шумовой добавки III). Здесь, как и раньше, lp - число частиц в ячейке. Число узлов сетки n = 100, время t = 0.1, = 0.001. Из таблицы 2.1 видно, что для метода частиц с вычитанием шумовой добавки величина уровня шума на два-три порядка меньше, чем для обычного метода частиц как в плотности, так и в средней скорости.

Таблица 2.1: Зависимость уровня шума от количества частиц для метода частиц в ячейках (величины с индексом pic) и для метод частиц в ячейках с Алгоритмом III (индекс new).

–  –  –

2.7.2 Сравнение схем для уравнений на v, n.

Необходимо отметить, что на решение в центральной области влияет выбор схемы для нахождения v, n. На рисунках 2.11 и 2.10 приведены результаты расчетов с использованием стандартного метода частиц без вычитания шума, противопотоковой схемы (2.34)-(2.38) и схемы предиктор-корректор (2.39)-(2.42) для параметров tm = 0.3, = 0.001.

Параметры метода частиц: число частиц m = 10000, число узлов сетки n = 101.

Параметры конечно-разностной схемы: число узлов сетки по x и u равно im = 1600 и km = 3000 соответственно.

Рис. 2.10: Зависимость качества решения от выбранной схемы. а - плотность, б - средняя скорость. Красным цветом - решение, полученное по схеме Лакса-Вендроффа для уравнения Власова.

Уменьшить уровень шума в областях A и B позволяет алгоритм с использованием любой из рассмотренных схем. Но из этих графиков видно, что использование схемы предикторкорректор (рисунки 2.11 и 2.10, черным цветом) меньше влияет на поведение решения в области распространения возмущений, чем противопотоковая схема. Использование противопотоковой схемы подавляет передний пик в графиках плотности и средней скорости (рисунки 2.11 и 2.10, зеленым цветом). Это связано с тем, что противопотоковая схема имеет Рис. 2.11: а - плотность, б - средняя скорость для различных схем. Красным цветом - решение, полученное по схеме Лакса-Вендроффа для уравнения Власова.

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

–  –  –

В данной главе на примере задачи распада разрыва плотности ионов предложен алгоритм уменьшения счетных шумов метода частиц. Предложенный алгоритм был реализован в программе на языке Fortran для решения одномерной задачи распада разрыва плотности ионов. Проведенное исследование метода показало его перспективность. Алгоритм вычитания шумовой добавки позволяет существенно уменьшить разброс в средней скорости и плотности в методе частиц в ячейках для поставленной выше задачи. Качество решения зависит от выбора подходящей схемы для уравнений нахождения v, n Алгоритма вычитания шума III.

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

На проблему самовоздействия обращают мало внимания. В большом ансамбле частиц она теряется на фоне других численных шумов и выделить ее не всегда возможно. Механизм возникновения самосилы показан в работах [101, 119] на примере моделирования методом частиц в ячейках движения одной частицы в отсутствии внешних полей. Также авторами показана связь саморазогрева в ансамбле модельных частиц модельной плазмы с погрешностью аппроксимации силы, действующей на каждую частицу. В работе [120] приведено сравнение величины самосилы для разных ядер на примере одномерной задачи распада разрыва в дисперсионной среде. В работе [104] авторы рассматривают проблему самосилы для метода частиц в ячейках с адаптивной сеткой, и предлагают алгоритм, который позволяет уменьшить величину самосилы до любой требуемой степени точности. При этом авторы (как и большинство других) утверждают, что проблемы самосилы для методов с неадаптивной сеткой не существует, ссылаясь на работы Р. Хокни, Дж. Иствуда, Ч. Бэдсела и А. Ленгдона [51, 52] и утверждая, что если функции распределения заряда в узлы Эйлеровой сетки и интерполяции силы в местоположение частицы одинаковы, то самовоздействие равно нулю.

Но данное условие не является достаточным. В [51, 52] утверждается, что самосила будет равна нулю при выполнении следующих условий

1) функции для распределения заряда и интерполяции силы одинаковы, т.е. SF = R (S = R в обозначениях монографии Хокни и Иствуда [51]),

2) разностные аппроксимации производных правильно сцентрированы по пространству [51], т.е. разностные уравнения, связывающие сеточные функции E и симметричны в пространстве [52].

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

–  –  –

= 4, (3.1) где = e(ni ne ) в случае двухкомпонентной водородной плазмы, e – заряд электрона,

–  –  –

3.1 Самосила и VSP-ядро в одномерном случае Рассмотрим механизм возникновения самовоздействия в методе частиц и предложенное в [101] VSP-ядро, устраняющее самосилу в одномерном случае.

Пусть модельная частица находится в ячейке [xi1, xi ] и имеет координату xj = xi1 + h (h - шаг сетки, 0 1). По любой из моделей ее заряд q распределяется между несколькими соседними узлами, в зависимости от вида функции R(x). Для ядра PIC

–  –  –

где (x, y) - координаты частицы, (xi, yk ) - координаты узла сетки.

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

–  –  –

Требуется вычислить величину электрического поля, действующего на эту частицу в ячейке (i, k) для трех приведенных ядер.

Значение плотности, полученное по QCP1 и QCP2 – ядрам, можно выразить через плотности (полученные по ядру PIC) следующим образом

–  –  –

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

–  –  –

3.2.4 Потенциал поля одиночного заряда Для частицы с единичным зарядом, расположенной в узле (i,k), значение плотности i,k = 2, в остальных узлах сетки плотность равна нулю. Уравнение Пуассона на сетке для h него будет иметь вид

–  –  –

где Q(i, k) = 1 в узле(i, k), в остальных узлах Q = 0. Данное уравнение решалось комбинированным итерационным методом (subroutine poisson2 в [81]): последовательная верхняя релаксация в направлении x, прогонка в направлении y. Граничные условия определяются по формуле (3.24), вывод которой приводится ниже. Приведем полученные значения потенциала точечного заряда (так как потенциал определяется с точностью до константы, из всех i отнято для удобства 0 ) 0 = 0, 1 =, 2 = 4, 3 = 4, 5664, 4 = 4, 8584, 5 = 5, 3333.

3.2.5 Сравнение ядер PIC, QCP1, QCP2 На рисунке 3.5 изображены направление и величина вектора напряженности электрического поля, действующего на частицу, в зависимости от положения частицы в ячейке. Из рисунка видно, что на частицу с PIC-ядром действует сила, притягивающая ее к центру ячейки, частица с QCP1-ядром также притягивается к центру, частицу же с формой ядра QCP2 выносит на грань ячейки.

Рис. 3.5: Напряженность электрического поля (направление и модуль), действующего на частицу, в зависимости от положения частицы в ячейке для PIC, QCP1 и QCP2 ядер.

В таблице 3.1 представлено сравнение самосилы для PIC, QCP1 и QCP2 ядер в двумерном случае. В средней колонке указано число узлов с ненулевым значением плотности для данного ядра, в правой – максимум по ячейке модуля вектора напряженности электрического поля.

Таблица 3.1: Максимум модуля напряженности электрического поля (max|E|) по ячейке для разных ядер.

–  –  –

Из таблицы видно, что максимальная величина самосилы получается при использовании PIC ядра. QCP1 и QCP2 по сравнению с PIC показывают лучшие характеристики.

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

–  –  –

Формулы (3.13) и (3.14) отличаются только коэффициентом. Значит можно попытаться так скомбинировать QCP1 и QCP2 ядра, чтобы самосила у нового ядра была нулевой.

Возьмем комбинацию ядер со следующими коэффициентами

–  –  –

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

–  –  –

Для тестирования нового ядра (1 = 0.1728, 2 = 0.0772) рассмотрим задачу движения одной и двух модельных частиц в электростатическом поле. Движение частиц рассчитывалось с использованием противопотоковой схемы, а уравнение Пуассона решалось как и выше комбинированным итерационным методом: последовательная верхняя релаксация в направлении x, прогонка в направленнии y.

Граничные условия для потенциала были получены следующим образом. В цилиндрической системе координат уравнение Пуассона имеет вид

–  –  –

где r(i, k) - расстояние от граничного узла с индексом (i, k) до частицы. Отметим важный факт. Из уравнения (3.24) видно, что все решения подобны, т.е. при выборе другого шага по пространству происходит линейное увеличение r в s раз. Тогда

–  –  –

а значит разности потенциалов в двух соседних узлах, и, как следствие E не меняются при изменении h.

На рисунке 3.6 показаны результаты расчетов траектория движения одиночной модельной частицы с разной формой ядра, а также с PIC ядром и билинейной интерполяцией силы по методу Хокни. Внешние поля отсутствуют, начальная скорость частицы равна (ux, uy ) = (2, 0) и (ux, uy ) = (1, 1). Видно, что при использовании PIC-ядра (ширины h и 2h) с билинейной интерполяцией для силы на смещенной сетке траектория частицы сильно искажается. Частица летит в заданном направлении, но в пределах ячейки (ширина ячейки сетки в тестах h = 0, 2) частица совершает колебания. Такие же колебания, но с меньшей амплитудой наблюдаются для ядра QCP1 (для QCP2 картина аналогична). Для нового ядра (NEW) и для метода Хокни с PIC ядром траектория частицы - прямая линия, что и говорит об отсутствии влияния сетки на динамику частицы.

Рис. 3.6: Траектория движения одиночной частицы в отсутствии внешних сил. (ux, uy ) = (2, 0) и (ux, uy ) = (1, 1) Рассмотрим движение двух неподвижных в начальный момент времени частиц с зарядами разного знака. Начальное положение частиц (x1, y1 ) = (10.5, 9.57) и (x2, y2 ) = (10.5, 10.43), шаг сетки h = 0.2. Две разнозаряженные частицы в отсутствии внешних полей должны совершать колебания относительно центральной точки. На рисунке 3.7 представлена y-компонента координаты траектории частиц в зависимости от времени. Видно, что при использовании PIC-ядра (PIC, линия с точками) из-за ошибки в расчете силы частицы не могут вылететь из ячеек. Сила, действующая на частицу со стороны сетки оказывается больше, чем сила от соседней частицы. Их движение абсолютно нефизично. Если использовать метод Хокни (центральные разности для производной потенциала), то частицы начинают колебаться (пунктирная линия), но центр этих колебаний постепенно меняется. Это связано с возникающим небольшим отклонением по оси x, которое понемногу накапливается. Частицы же с новым ядром (NEW) совершают колебания, центр которых остается практически неизменным.

Рис. 3.7: y-компонента координаты траектории движения двух частиц с зарядом разного знака в зависимости от времени в отсутствие внешних полей.

–  –  –

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

Причина их возникновения может заключаться в следующем. В [18] приведено сравнение дисперсионного уравнения для плазмы с точечными частицами

–  –  –

в одномерном электростатическом случае продольных колебаний малой амплитуды. Из этих соотношений видно, что учет конечного размера частиц приводит в дисперсионном уравнении к замене 0e на 0e R(k). Поэтому дисперсионные свойства плазменных колебаний могут изменяться в зависимости от вида ядра. Исследовать аналогичным образом в более сложном случае влияние на дисперсионное уравнение частицы с ядром конечного размера невозможно. Но в случае, если в дисперсионное уравнение 0e R(k) входит в нечетной степени, смена знака этого выражения может приводить к развитию неустойчивости. Этот факт необходимо учитывать.

Образы Фурье разных ядер в одномерном случае имеют следующий вид:

–  –  –

Из формул и рисунка 3.8 видно, что образ Фурье ядра PIC всегда положителен, поэтому свойства устойчивости не меняются. Ядро VSP (как и NGP) имеет при некоторых k отрицательные значения, Rmin 0, 5. Это приводит к тому, что некоторые гармоники становятся неустойчивыми. Для параболического ядра (интерполяционный многочлен 2-го порядка) Rmin 0, 01, и слабая неустойчивость отдельных гармоник незаметна на фоне других счетных погрешностей. То же самое справедливо и в двумерном случае.

Поэтому было решено рассмотреть ядро, которое представляет собой комбинацию QCP1 и QCP1 с PIC ядром

–  –  –

3.4.2 Выбор оптимальных параметров Как было показано выше на примере одномерного случая, для того, чтобы ядро было устойчивым, необходимо, чтобы его образ Фурье был положительной функцией (R(k, m) 0) или по крайней мере не принимал больших отрицательных значений. Фурье-образ комбинированного ядра равен

–  –  –

Переобозначим x = kh, y = mh. На рисунке 3.9 показан вид функции G(x, y) для произвольных значений 1, 2. Необходимо найти область параметров 1, 2, в которой G(x, y) 0. Это требование сводится к тому, чтобы найти область параметров, в которой все минимумы функции G(x, y) неотрицательны.

Далее будем рассматривать G как функцию от (1, 2).

1. Рассмотрим плоскость x = 0 (для y = 0 аналогично). Обозначим t := x или y.

–  –  –

линейная функция, как и по 1, 2. Значит минимум ее нужно искать на границе допустимой по параметрам 1, 2 области (рисунок 3.10).

Для того, чтобы решить систему (3.31), необходимо найти минимум по (1, 2 ) максимума по d [0, 1] модуля функции F, то есть min(1,2 ) maxd |F |.

Для любых (1, 2 ) maxd |F | достигается при d = 0. Найдем точку на области параметров (1, 2 ), в которой |F | при d = 0 минимален на границе.

• 0 1 0.5, 2 = 0.

–  –  –

Таким образом минимум самосилы достигается при параметрах 1 = 0.25, 2 = 0.5, 0 = 1 1 2. Форма ядра с этими параметрами (назовем его QCP) приведена на рисунке 3.11.

–  –  –

Обозначим через max|E| максимум в ячейке модуля вектора напряженности электрического поля. Для QCP max|E| = 1, 1107/h, и следовательно величина самосилы для QCP ядра с указанными выше параметрами, по сравнению с PIC ядром, уменьшилась в 4 раза.

3.5 Саморазогрев модельной плазмы Свойства нового ядра были проверены на эффект саморазогрева модельной плазмы.

Хорошо известно [121], что в реальной плазме хаотические поля вызывают стохастический нагрев плазмы. А поскольку конечность шага по времени и пространственного шага, используемых в методе частиц, приводит к возникновению случайных полей (из-за ошибок аппроксимации) и к нарушению законов сохранения, то и в модельной плазме наблюдается постепенный рост полной энергии, и, соответственно, температуры [100], что и получило название «саморазогрев». Расчеты на саморазогрев проводились по двумерной модели аналогично численным экспериментам на саморазогрев в монографии [51].

Решается система, состоящая из уравнений движения частиц и уравнения Пуассона для потенциала электрического поля

–  –  –

Плазма состоит из ионов и электронов, магнитного поля нет (B = 0). Граничные условия по всем направлениям периодические. Уравнение Пуассона решается методом преобразования Фурье. Для уравнений движения частиц используется противопотоковая схема первого порядка. Отношение масс электронов и ионов выбрано me /mi = 100, плотность ne = ni = 1.

Размер области 2x2, cетка 32x32, шаг по времени = 0.025. Общее число частиц 3 · 106.

Частицы распределены в начальный момент времени равномерно по области, распределение по скоростям максвелловское (с дисперсией электронов ve = 0.002, ионов vi = 0.0002), средняя скорость нулевая v0 = 0.

–  –  –

На рисунке 3.12 приведена величина приращения полной энергии |(W (t)W (0))/W (0)| в зависимости от времени для описанной задачи с использованием различных ядер. Видно, что для параболического ядра, обладающего большей гладкостью, чем PIC-ядро, полная энергия растет медленнее. Использование QCP ядра позволило добиться результата, сравнимого с параболическим ядром. Но формулы для него более простые, по затратам этот метод сравним с PIC. Метод Хокни с PIC ядром на этой же задаче показал повышение темпа роста полной энергии. Таким образом можно заключить, что использование QCP-ядра позволяет лучше сохранять полную энергию системы, чем PIC-ядро.

3.6 Моделирование эволюции протопланетного диска с QCP- ядром

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

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

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

За основу была взята программа Вшивкова В.А. и Снытникова А.В, описанная в статье [122], в которой PIC-ядро заменено на QCP-ядро. Рассматривается движение только пылевой компоненты диска без газа.

Эволюция пылевой компоненты протопланетного диска, вращающегося в гравитационном поле некоторой звезды описывается системой уравнений звездной динамики, состоящей из уравнения Власова-Лиувилля и уравнения Пуассона:

–  –  –

где f (t, r, u) – функция распределения частиц пыли по скоростям, a – ускорение частиц, которое определяется формулой a =, - объемная плотность вещества, - гравитационный потенциал, - гравитационная постоянная. Частицы движутся в трех измерениях, но при изучении эволюции протопланетного диска предполагается наличие центрального тела с массой, на порядок большей массы диска. В этом случае движение может считаться двумерным, а сам диск бесконечно тонким, так что объемная плотность в правой части уравнения Пуассона равна 0. Плотность вещества в плоскости диска используется в граничных условиях. Уравнение Пуассона, тем не менее, является строго трехмерным. В задаче моделирования эволюции протопланетного диска область решения имеет вид кругового цилиндра (система координат цилиндрическая)

0 r RM, 0 2, 0 z ZM,

при этом диск считается бесконечно тонким и находится на нижней грани цилиндра (z = 0).

Такая модель позволяет изучать процессы структуризации протопланетного диска и звездных скоплений под действием сил гравитации. В отсутствие газа или если его плотность пренебрежимо мала (не влияет на динамику системы), численная модель позволяет изучить коллективные процессы самоорганизации во вращающемся диске отдельных частиц без парных столкновений. Известно, что холодный твердотельно вращающийся диск из частиц неустойчив относительно развития крупномасштабных джинсовских возмущений. Поэтому на временах, меньших определяемых инкрементом развития джинсовской неустойчивости, в неустойчивом диске не должны проявляться физически неустойчивые моды. Ясно, что на фоне физической неустойчивости сложно заметить какие-то другие неустойчивости, в том числе связанные с численными методами. Неустойчивые моды непосредственно от самого численного метода можно выявить на начальных этапах процесса, поэтому рассмотрим поведение прироста полной энергии для времен t = 4 и t = 8 (один и два оборота диска соответственно), для разной начальной дисперсии скоростей частиц. Параметры, при которых проводились расчеты: начальный радиус диска r1 = 1, радиус расчетной области Rm = 6, сетка 121x128, 400000 частиц.

На рисунке 3.14 приведен прирост полной энергии |W W 0|/|W 0| (W 0 = W (0)) в зависимости от времени при начальной дисперсии скоростей v = 0.8 и v = 4, а на рисунке 3.13 Рис. 3.13: Плотность пылевой компоненты диска в момент времени t = 8: (а ) v = 0.8 и (б) v = 4.

Рис. 3.14: Прирост полной энергии при начальной дисперсии: (а ) v = 0.8 и (б) v = 4.

плотность пылевой компоненты диска. Как видно из графика плотности, в случае v = 0.8 (рисунок 3.13 а) начинают образовываться сгустки вещества, при v = 4 (рисунок 3.13 б) диск разлетается. Характер решения не зависит от вида модельного ядра. В то же самое время из рисунка 3.14 видно, что почти везде энергия лучше сохраняется для QCP-ядра.

Резкий рост энергии для PIC-ядра на рисунке 3.14 а) объясняется недостаточно хорошим учетом энергии частиц, попавших в центр, поэтому сравнивать кривые на первом графике имеет смысл только до одного оборота диска. Использование QCP ядра при моделировании эволюции протопланетного диска на нескольких тестах показало лучшее или аналогичное PIC сохранение полной энергии на временах порядка одного-двух оборотов диска.

3.7 Выводы

Показано, что использование нового QCP-ядра, предложенного в данной главе, приводит к лучшему, по сравнению с PIC-ядром, сохранению импульса частиц и полной энергии, и, как следствие, к увеличению времени саморазогрева модельной плазмы. Тот факт, что плотность, вычисляемая по новому ядру, выражается через плотности по PIC-ядру (3.25), позволяет легко модифицировать программу с PIC-ядром под QCP-ядро. Кроме того, вычисление плотности таким образом требует меньше затрат по времени, чем вычисление по параболическому и более сложным ядрам, так как количество необходимых дополнительных (сверх PIC) операций пропорционально размеру сетки, а не количеству модельных частиц, что важно при большом количестве частиц. Таким образом использование QCP-ядра экономичнее и по вычислительным затратам, и по времени реализации, чем ядра более высоких, чем PIC, порядков. Использование QCP-ядра можно охарактеризовать и как дополнительное промежуточное сглаживание на сетке. Отличие от обычного произвольного сглаживание заключается в том, что свойства данного метода исследованы и показано, что выбранные параметры являются оптимальными.

Глава 4 Число модельных частиц и точность на примере задачи моделирования кинетической неустойчивости теплого электронного пучка малой плотности в плазме Пучковая неустойчивость является одной из наиболее распространенных неустойчивостей в плазме и представляет большой интерес для широкого круга приложений, таких как ускорение заряженных частиц внешними электромагнитными полями и нагрев лабораторной плазмы. В настоящее время известны разнообразные пучковые неустойчивости. Вызывается пучковая неустойчивость резонансным взаимодействием по меньшей мере двух компонент заряженных частиц: одной из компонент являются частицы пучка, другой - частицы среды, сквозь которую распространяется пучок. Неустойчивость состоит в экспоненциальном нарастании со временем флуктуаций плотности заряда и электромагнитных полей. Экспоненциальный рост имеет место на линейной стадии развития неустойчивости. Когда же нелинейные процессы ограничивают этот рост, говорят о стадии насыщения неустойчивости. По виду частиц, вызывающих неустойчивость, различают электрон-электронные, электрон-ионные и пирсовские неустойчивости [123].

В рамках работы, проводимой с ИЯФ СО РАН по совместным проектам, была создана трехмерная численная модель взаимодействия пучка электронов с плазмой, основанная на методе частиц и ориентированная на исследование устойчивости и нагрева плазмы электронным пучком. Несмотря на то что задача о релаксации электронного пучка в плазме является классической проблемой физики плазмы [124, 125] и существует целый ряд теоретических моделей, описывающих различные режимы пучково-плазменного взаимодействия, исследования в этой области остаются актуальными [126,127]. Максимально приближенная к условиям лабораторных экспериментов постановка задачи зачастую требует отказа от привычных для теории идеализаций, таких как слабое или сильное магнитное поле, гидродинамический или кинетический характер пучковой неустойчивости, приближение случайных фаз возбуждаемых в плазме турбулентных пульсаций. Кроме того, при длительной инжекции пучка эволюция пучково-плазменной системы может проходить через целую последовательность стадий, определяемых совершенно разными нелинейными процессами.

Работа была направлена на изучение сценариев возбуждения плазменной турбулентности электронным пучком.

При этом наибольший интерес представляют параметры пучка и плазмы, характерные для экспериментов по нагреву плазмы в открытой ловушке ГОЛ-3 (ИЯФ СО РАН) [128]. Для исследования влияния пучковых нелинейностей на поведение неустойчивости в условиях развитой турбулентности необходимо численное моделирование, которое, с одной стороны, способно на больших временах отслеживать эволюцию возбуждаемой пучком турбулентности, а с другой - позволяет обеспечить достаточно подробное описание кинетических эффектов, связанных с захватом пучка.

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

–  –  –

Модель высокотемпературной бесстолкновительной плазмы представляется кинетическим уравнением Власова и системой уравнений Максвелла, которые в безразмерной форме имеют следующий вид [129]:

–  –  –

Здесь индексами i и e помечены величины, относящиеся к ионам и электронам, соответственно; qe = 1, qi = me /mi ; fi,e (t, r, p) – функция распределения частиц; mi,e, pi,e, ri,e

– масса, импульс, положение иона или электрона; E, B – напряженности электрического и магнитного полей. Для перехода к безразмерному виду в качестве единиц используются следующие базовые величины:

• скорость света c = 3 · 1010 см/с;

• масса электрона me = 9, 1 · 1028 г;

• плотность плазмы n0 = 1014 см

• время t = pe, где плазменная электронная частота pe = 5, 6 · 1011 c1.

–  –  –

находится плазма, состоящая из электронов и ионов водорода, и пучок электронов. Задаются плотности пучка nb и электронов плазмы ne = 1 nb, температура электронов плазмы Te и пучка Tb ; температура ионов считается нулевой Ti = 0. Начальное распределение частиц по скоростям максвелловское с плотностью распределения

–  –  –

Для решения системы уравнений Максвелла используется метод Лэнгдона-Лазинского, описанный в работе [130], в котором поля определяются из разностных аналогов законов

Фарадея и Ампера:

–  –  –

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

–  –  –

m+1 m + divh jm+1/2 = 0, (4.9) то разностные аналоги уравнений (4.4) и (4.5) выполняются автоматически. Для того, чтобы удовлетворить этим условиям, токи вычисляются по схеме, приведенной в статье [112], магнитное поле в начальный момент нулевое B0 = 0.

Большой объем данных, связанных с трехмерностью задачи, требует использования распределенных вычислительных систем. Распараллеливание выполнено методом декомпозиции расчетной области по направлению y, перпендикулярному направлению движения электронного пучка x (рисунок 4.1). Используется смешанная эйлерово-лагранжева декомпозиция. Сетка, на которой решаются уравнения Максвелла, разделена на одинаковые подобласти по одной из координат. С каждой подобластью связана группа процессоров. Далее, модельные частицы каждой из подобластей разделяются между процессорами связанной с этой подобластью группы равномерно, вне зависимости от координат. Каждый из процессоров группы решает уравнения Максвелла во всей подобласти. Далее решаются уравнения движения модельних частиц. После этого происходит суммирование значений тока по всей области. Один из процессоров группы производит обмен граничными значениями тока и полей с соседними подобластями, и затем рассылает полученные граничние значения всем процессорам своей группы. Подробности параллельной реализации описаны в статье Снытникова А.В. [131]. В настоящей диссертационной работе они не приводятся, так как не являются частью представляемых результатов автора диссертации.

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

4.3 Вычисление инкремента неустойчивости

Знание линейного инкремента неустойчивости важно для определения типа наиболее неустойчивых возмущений. Величина инкремента существенно влияет на установившуюся мощность накачки энергии. Таким образом точное вычисление линейного инкремента является важной задачей. Для оценки точности получаемого решения проводится сравнение инкремента развивающейся неустойчивости с его аналитическими оценками для той же задачи в одномерной постановке [132]. Известно, что на начальной стадии развития неустойчивости энергия электрического поля (W ) нарастает по экспоненте, W e2t.

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

1 = ln(W ). (4.10) 2 t Расчеты, однако, показали, что данный способ вычисления инкремента неустойчивости не подходит для вычислений в кинетическом режиме. Сложность перехода от гидродинамического режима к кинетическому режиму состоит в следующем. Когда разброс по скоростям электронов пучка увеличивается, а плотность пучка становится все меньше, во взаимодействие с возбуждаемой волной вступает все меньшая часть электронов пучка, и время насыщения неустойчивости увеличивается. При этом энергия электрического поля волны уменьшается и при недостаточном количестве частиц становится по величине сравнима с численными шумами. Такие расчеты требуют большого количества частиц, что не всегда возможно обеспечить. Поэтому важно, во-первых, подобрать менее чувствительный к численным шумам способ вычисления инкремента неустойчивости, а во-вторых, ответить на вопрос, сколько частиц в ячейке можно считать достаточным для вычислений.

В качестве такой диагностики был выбран расчет инкремента неустойчивости по амплитуде отдельной неустойчивой моды E0,

–  –  –

На графиках, представленных ниже, приведено сравнение результатов расчетов инкремента неустойчивости по энергии электрического поля (4.12) и по амплитуде главной моды (4.13) с решением, описанным в [132] (для одномерной задачи). Расчеты проводились при параметрах, соответствующих трем разным режимам развития пучковой неустойчивости.

Счетные параметры: сетка по пространству 100 4 4, шаг по времени = 0.001. Обозначим lp число частиц в одной ячейке сетки и рассмотрим вопрос точности получаемого решения в разных режимах в зависимости от lp, а также вопрос ранней идентификации неустойчивого решения.

1) Гидродинамический режим (kv ).

Основные параметры: отношение плотности пучка к плотности плазмы nb /n0 = 2·103, разброс электронов пучка по скоростям v/v0 = 0.035, длина области L = 1.2566. При заданных параметрах величина инкремента неустойчивости должна быть = 0.0722 [132].

На рисунке 4.2 показан темп роста энергии электрического поля и значение инкремента 1, вычисленное по формуле (4.12). Здесь и далее на графиках с индексом а прямой пунктирной линией показан темп роста волны с аналитическим инкрементом, положение этой линии снесено по координате t, максимумы графиков совмещены по времени. На графиках с индексом б прямой пунктирной линией показано точное значение инкремента. На графиках видна одна из особенностей моделирования холодной плазмы методом частиц - саморазогрев модельной плазмы, который выражается в том, что при малых температурах плазмы в самом начале вычислений энергия возрастает до определенного значения [51], величина которого зависит от счетных параметров. Как видно из рисунка 4.2, уровень саморазогрева уменьшается с увеличением числа частиц. Когда амплитуда развивающейся неустойчивой волны становится больше уровня шумов, связанных с саморазогревом, экспоненциальная стадия роста волны становится видна на графиках (линейный рост кривой на рисунке 4.2 а). Далее Рис. 4.2: Временная зависимость логарифма энергии электрического поля (а) и производная логарифма энергии электрического поля (инкремент 1 ) (б) при разном числе частиц в ячейке. Гидродинамический режим.

Рис. 4.3: Временная зависимость логарифма амплитуды электрического поля (а) и производная логарифма амплитуды (инкремент 2 ) (б) при разном числе частиц в ячейке. Гидродинамический режим.

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

Второй способ вычисления инкремента (по амплитуде главной волны электрического поля, формула (4.11)) показан на рисунке 4.3. Здесь, как видно из графиков, не играют роли шумы от саморазогрева модельной плазмы, что позволяет практически сразу наблюдать экспоненциальный рост. Более того, как видно из представленного графика, для той же точности вычисления инкремента можно использовать меньшее количество частиц.

Таблица 4.1: Инкремент и ошибка в % в зависимости от количества частиц в ячейке.

Гидродинамический режим.

–  –  –

25 0.058 20 0.066 9 50 0.062 14 0.068 6 250 0.067 7 0.069 4 500 0.068 6 0.07 3 2500 0.07 3 0.071 2 5000 0.07 3 0.071 2 В таблице 4.1 приведено значение инкремента, вычисленного по энергии и по амплитуде, а также характерная величина ошибки в процентах. В гидродинамическом случае ошибка по энергии примерно в 2 раза больше, чем по амплитуде. Из приведенных графиков, а также из таблицы видно, что в гидродинамическом режиме достаточным можно считать число частиц в ячейке, равное 25-50, если использовать вычисление инкремента по амплитуде.

Ошибка в этом случае не превышает 10%.

2) Переходный режим (kv ) Основные параметры: отношение плотности пучка к плотности плазмы nb /n0 = 2·103, разброс электронов пучка по скоростям v/v0 = 0.14, длина области L = 1.1424. При этих параметрах величина инкремента неустойчивости должна быть = 0.0232.

Рис. 4.4: Временная зависимость логарифма энергии электрического поля (а) и производная логарифма энергии электрического поля (инкремент 1 ) (б) при разном числе частиц в ячейке. Переходный режим.

Рис. 4.5: Временная зависимость логарифма энергии электрического поля (а) и производная логарифма энергии электрического поля (инкремент 2 ) (б) при разном числе частиц в ячейке. Переходный режим.

Величина инкремента неустойчивости по сравнению с гидродинамическим режимом уменьшилась в три раза. Таким образом энергия нарастает медленнее, насыщение неустойчивости происходит позже по времени. Как и в гидродинамическом режиме, на графике логарифма энергии электрического поля (рисунок 4.4) виден саморазогрев модельной плазТаблица 4.2: Инкремент и ошибка в % в зависимости от количества частиц в ячейке. Переходный режим.

–  –  –

250 0.017 27 0.02 14 500 0.019 18 0.02 14 2500 0.02 14 0.021 9 5000 0.021 9 5 0.022 мы. На графике амплитуды основной моды (рисунок 4.5) этого эффекта нет. Из приведенных рисунков и таблицы 4.2 видно, что в переходном режиме необходимо использовать гораздо большее число частиц в ячейке по сравнению с гидродинамическим режимом. При увеличении числа частиц значение инкремента приближается к точному. При недостаточном числе частиц (lp = 50) график инкремента 2 (рисунок 4.5 б) не имеет области, близкой к горизонтальной прямой, которая была четко видна в гидродинамическом расчете и которая появляется в переходном режиме при большем числе частиц. При lp = 250 область экспоненциального роста видна лучше. Далее с ростом числа частиц кривая ln(E0 )/2 в интервале [0,150] приближается к прямой линии, а 2 (t) выходит на уровень, примерно соответствующий точному значению инкремента. Таким образом в данном случае для того, чтобы вычислить инкремент, нам необходимо использовать для расчетов как минимум 250 частиц в ячейке.

По ресурсам вычисления становятся более затратными, считать приходится на большие времена с большим числом частиц, чем в первом расчете. Из таблицы 4.2 видно, что ошибка по энергии также в 1.5-2 раза больше, чем по амплитуде. Для того, чтобы ошибка вычисления инкремента по амплитуде была в пределах 10%, число частиц в одной ячейке должно быть lp 2500. Чтобы достичь той же точности при вычислении инкремента по энергии электрического поля, потребуется 5000 частиц в одной ячейке.

3) Кинетический режими (kv ).

Основные параметры: отношение плотности пучка к плотности плазмы nb /n0 = 2·104, разброс электронов пучка по скоростям v/v0 = 0.14, длина области L = 1.1424. Аналитическое значение инкремента неустойчивости в этом случае = 0.0027.

Рис. 4.6: Временная зависимость логарифма энергии электрического поля (а) и производная логарифма энергии электрического поля (инкремент 2 ) (б) при разном числе частиц в ячейке. Кинетический режим.

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

Как и в переходном режиме, при недостаточном числе частиц график инкремента 2 при lp = 500 не имеет области, близкой к горизонтальной прямой, и, соответственно, выбрать точку для вычисления инкремента невозможно. С ростом числа частиц область экспоненциального роста становится шире. Для вычисления значения инкремента в кинетическом режиме с ошибкой в 25%, требуется около 5000 частиц в ячейке. Отметим, что, в отличие от гидродинамического режима, графики в переходном и кинетическом режиме приведены со сглаживанием.

4.5 Фазовые плоскости

Помимо инкремента, признаком развития неустойчивости является характерное поведение электронов пучка на фазовой плоскости (x, vx ). Для визуализации этого поведения отслеживалась эволюция выбранных частиц, начальное положение которых представлено на рисунке 4.7 (верхний ряд). В нижнем ряду показано положение электронов пучка в моменты времени, соответствующие насыщению неустойчивости. Для кинетического режима (рисунок 4.7 в) показана только узкая область по скорости 0.17 vx 0.2, но ширина начального распределения частиц такая же, как в переходном режиме (рисунок 4.7 б, t = 0).

Видно, во-первых, что в гидродинамическом режиме возбуждаемой волной захватывается весь пучок, а в переходном и кинетическом все более узкая его часть, остальные же электроны во взаимодействие с волной не вступают. Этим, в частности, и обусловлавливается то, что в кинетическом режиме требуется брать такое большое число частиц. Несмотря на малое (с точки зрения точности вычисления инкремента в кинетическом режиме) количество частиц, видно, что происходит захват именно тех частиц, скорость которых близка к расчетной фазовой скорости волны v = 0.9v0 = 0.18. Это соответствует аналитическим оценкам для одномерной задачи [133].

Рис. 4.7: Электроны пучка на фазовой плоскости (x, vx ) для а - гидродинамического (v = 0.007, lp = 50), б - переходного (v = 0.02, lp = 500) и в - кинетического (v = 0.02, lp = 5000) режимов.

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

Обозначим vx ширину по vx области захваченных частиц. Для переходного и кинетического режима она равна vx = C0, k где значение константы C0 для всех расчетов в переходном и кинетическом режимах одинаково и равно C0 14, что соответствует теоретическому значению 13.6 [133]. Следовательно, можно заключить, что построенная нами модель позволяет не только качественно, но и количественно правильно моделировать пучковую неустойчивость.

Рис. 4.8: Электроны пучка на фазовой плоскости (x, vx ) для а - переходного (v = 0.02, lp =

50) и б - кинетического (v = 0.02, lp = 500) режимов.

Для сравнения на рисунке 4.8 приведены фазовые портреты пучка в те же самые моменты времени для переходного и кинетического режимов при недостаточном числе частиц (lp = 50 и 500, соответственно), когда вычислить инкремент невозможно, так как невозможно выделить область экспоненциального роста. При недостаточном числе частиц область захваченных частиц формируется гораздо раньше и разрушается быстрее. Видно, что в этом случае картина не такая четкая. Помимо основной волны здесь возникает ряд волн сравнимой с ней амплитуды. Эти волны не просто маскируют решение, а разрушают его. Эти волны имеют чисто вычислительную природу. Опасность их заключается в том, что они не просто накладываются сверху как помеха, которую легко отделить от решения, но взаимодействуют друг с другом и с решением. А так как плазма это среда, в которой развиваются неустойчивости, со временем такие ошибки также могут расти экспоненциально, тем самым разрушая решение. В таком случае назвать полученные картины решением исходной задачи нельзя. Таким образом можно заключить, что невозможность вычислить инкремент по формуле (4.13) вызвана ошибками в моделировании явления, а не малой чувствительностью диагностики. Следовательно, совпадение рассчитаного по формуле (4.13) инкремента с теоретически предсказанным значением может служить критерием соответствия результатов моделирования физическому явлению, а также количественной мерой его точности.

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

4.6 Выводы

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

Рис. 4.9: Электроны пучка на фазовой плоскости (x, vx ) для гидродинамического (v = 0.007, lp = 500) и переходного (v = 0.02, lp = 5000) режимов.

Заключение

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

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

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

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

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

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

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

• На основе метода частиц в ячейках создан комплекс программ для моделирования взаимодействия электронного пучка с плазмой в трехмерной постановке при параметрах, соответствующих условиям лабораторных экспериментов на установке ГОЛ-3.

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

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

Литература [1] Вшивков, В.А. Форма ядра частицы и проблема «самовоздействия» в методе частиц-вячейках / В.А. Вшивков, Е.А. Месяц // Научный вестник НГТУ. – № 1 (42). – 2011. – С. 47–56.

[2] Месяц, Е.А. О выборе числа частиц в методе частиц-в-ячейках для моделирования задач физики плазмы / Е.А. Месяц, А.В. Снытников, К.В. Лотов // Вычислительные технологии. – Т. 18. – № 6. – 2013. – С. 83–97.

[3] Mesyats, E.A. A noise-reducing algorithm for particle-in-cell plasma simulation / E.A. Mesyats // Bulletin of the Novosibirsk Computing Center. Ser.: Numerical Analysis. – 2009. – № 14.

–P. 21–30.

[4] Месяц, Е.А. Разработка алгоритма, уменьшающего шум в методе частиц-в-ячейках / Е.А. Месяц // Всероссийская конференция по математике и механике, посвященная 130летию Томского государственного университета и 60-летию механико-математического факультета: cборник тезисов. 22 – 25 сентября 2008. – Томск: ТГУ. – 2008. – С. 138–139.

[5] Месяц, Е.А. Исследование шумовых свойств метода частиц-в-ячейках / Е.А. Месяц // Дифференциальные уравнения. Функциональные пространства. Теория приближений.

Международная конференция, посвященная 100-летию со дня рождения С.Л. Соболева.

Тезисы докладов. 5 – 12 октября 2008. – Новосибирск. – 2008. – С. 529.

[6] Месяц, Е.А. Об одном новом ядре для метода частиц-в-ячейках / Е.А. Месяц // Труды конференции молодых ученых ИВМиМГ СО РАН. 21 – 24 апреля 2009. – Новосибирск:

ИВМиМГ СО РАН. – 2009. – С. 94–105.

[7] Вшивков, В.А. Исследование численных свойств метода частиц-в-ячейках / В.А. Вшивков, Е.А. Месяц // XIII Всероссийская молодежная конференция-школа «Современные проблемы математического моделирования. Математическое моделирование, численные методы и комплексы программ». – Ростов-на-Дону: ЮФУ. – 2009. – С. 149–156.

[8] Месяц, Е.А. Минимизация погрешности в методе частиц-в-ячейках / Е.А. Месяц // VI

Всесибирский конгресс женщин-математиков (в день рождения Софьи Васильевны Ковалевской): Материалы Всероссийской конференции. 15 – 18 января 2010. – Красноярск:

РИЦ СибГТУ. – 2010. – C. 284–287.

[9] Месяц, Е.А. Форма ядра частицы и проблема «самовоздействия» в методе частиц-вячейках / Е.А. Месяц // Труды XV Байкальской Всероссийской конференции «Информационные и математические технологии в науке и управлении»; ч. 1. 1 – 9 июля 2010. – Иркутск: ИСЭМ СО РАН. – 2010. – C. 197–203.

[10] Месяц, Е.А. Новое QCP-ядро с пониженной самосилой для метода частиц-в-ячейках и его использование при моделировании эволюции протопланетного диска / Е.А. Месяц // Труды конференции молодых ученых ИВМиМГ СО РАН. 12 – 14 апреля 2011. – Новосибирск: ИВМиМГ СО РАН. – 2011. – С. 62–74.

[11] Месяц, Е.А. Новое ядро с пониженной самосилой в методе частиц-в-ячейках и его использование при моделировании эволюции протопланетного диска [Электронный ресурс] / Е.А. Месяц // Всероссийская конференция по вычислительной математике КВМ-2011.

Тезисы докладов. 29 июня – 1 июля 2011. – Режим доступа:

http://www.sbras.ru/ws/show_abstract.dhtml?ru+220+16170 [12] Месяц, Е.А. Применение нового QCP-ядра с пониженным самовоздействием для метода частиц-в-ячейках при моделировании эволюции протопланетного диска / Е.А. Месяц // Российско-монгольская конференция молодых ученых по математическому моделированию, вычислительно-информационным технологиям и управлению, тезисы докладов.

Иркутск (Россия) – Ханх (Монголия). 17 – 21 июня 2011. – Иркутск: ИДСТУ СО РАН.

– 2011. – С. 55.

[13] Месяц, Е.А. Трехмерная численная модель насыщения двухпотоковой неустойчивости электронного пучка в плазме / Е.А. Месяц, А.В. Снытников // Тезисы докладов XIX Всероссийской конференции «Теоретические основы и конструирование численных алгоритмов и решение задач математической физики», посвященной памяти К.И. Бабенко.

Дюрсо. 10 – 16 сентября 2012. – Издательство ИПМ им. М.В. Келдыша. – 2012. – С. 73–74.

[14] Месяц, Е.А. Трехмерная численная модель релаксации электронного пучка в плазме / Е.А. Месяц, А.В. Снытников // Материалы XIII Всероссийской конференции молодых ученых по математическому моделированию и информационным технологиям. 15 – 17 октября 2012. – Новосибирск: ИВТ СО РАН. – 2012. – С. 28.

[15] Mesyats, E.A. Particle-in-cell simulation of kinetic instability of an electron beam in plasma / E.A. Mesyats, A.V. Snytnikov // Book of abstracts of the international conference «Mathematical modeling and computational physics». July 8-12, 2013. – Dubna. – 2013. – P. 129–130.

[16] Сигов, Ю.С. Вычислительный эксперимент: мост между прошлым и будущим физики плазмы. Избранные труды / Ю.С. Сигов; cоставители Г.И. Змиевская, В.Д. Левченко. – М.: ФИЗМАТЛИТ, 2001. – 288 с.

[17] Иванов, М.Ф. Численное моделирование динамики газа и плазмы методами частиц / М.Ф. Иванов, В.А. Гальбурт. – М.: Издательство МФТИ, 2000. – 168 с.

[18] Березин, Ю.А. Методы частиц в динамике разреженной плазмы / Ю.А. Березин, В.А.

Вшивков. Новосибирск: Наука, 1980. – 95 с.

[19] Ахиезер, А.И. Электродинамика плазмы / А.И. Ахиезер, И.А. Ахиезер, Р.В. Половин, А.Г. Ситенко, К.Н. Степанов; под ред. А.И. Ахиезера. – Mосква: Наука, 1974. – 720 с.

[20] Kwok, Dixon T.K. A hybrid Boltzmann electrons and PIC ions model for simulating transient state of partially ionized plasma / Dixon T.K. Kwok // Journal of Computational Physics. – 2008. – Vol. 227. – Iss. 11. – P. 5758–5777.

[21] Калиткинa, Н.Н. Математические модели физики плазмы / Н.Н. Калиткинa, Д.П. Костомаров // Математическое моделирование. – 2006. –Т. 18. № 11. – С. 67–94.

[22] Власов, А.А. Теория многих частиц / А.А. Власов. – М.: Гостехиздат, 1950. – 348 с.

[23] Арцимович, Л.А. Физика плазмы для физиков / Л.А. Арцимович, Р.З. Сагдеев. – М.:

Атомиздат, 1979. – 313 с.

[24] Байерс, Дж. Конечно-разностые методы в плазме без столкновений / Дж. Байерс, Дж.

Киллин // Вычислительные методы в физике плазмы / под ред. Б. Ольдера, С. Фернбаха, М. Ротенберга. – М.: Мир, 1974. – C. 259–303.

[25] Телегин, В.И. Об одной разностной схеме для уравнения Власова / В.И. Телегин // Журнал вычислительной математики и математической физики. – 1976. – Т. 16. – № 5.

– C. 1191–1197.

[26] Filbet, F. Comparison of eulerian Vlasov solvers / F. Filbet, E. Sonnendrucker // Computer Physics Communications. – 2003. – V. 150. – Iss. 3. – P. 247–266.

[27] Zaki, S.I. A finite element code for the simulation of one-dimensional vlasov plasmas. I.

Theory / S.I Zaki, L.R.T Gardner, T.J.M Boyd // Journal of Computational Physics. –1988.

– Vol. 79. – Iss. 1. – P. 184–199.

[28] Zaki, S.I. A finite element code for the simulation of one-dimensional vlasov plasmas. II.

Applications / S.I. Zaki, L.R.T. Gardner, T.J.M. Boyd // Journal of Computational Physics.

–1988. – Vol. 79. – Iss. 1. – P. 200-208.

[29] Boris, J.P. Flux-corrected transport. III. Minimal-error FCT algorithms / J.P. Boris, D.L.

Book // Journal of Computational Physics. – 1976. – Vol. 20. – Iss. 4. – P. 397–431.

[30] Fijalkow, E. A numerical solution to the Vlasov equation / E. Fijalkow // Computer Physics Communications. – 1999. – Vol. 116. – Iss. 2–3. – P. 319-328.

[31] Elkina, N.V. A new conservative unsplit method for the solution of the Vlasov equation / N.V. Elkina, J. Buchner // Journal of Computational Physics. – 2006. – Vol. 213. – Iss. 2. – P. 862–875.

[32] Filbet, F. Conservative numerical schemes for the Vlasov equation / F. Filbet, E.

Sonnendrucker, P. Bertrand // Journal of Computational Physics. – 2001. – Vol. 172. –Iss. 1.

–P. 166–187.

[33] Армстронг, Т. Решение уравнения Власова методом преобразований / Т. Армстронг, Р.

Хардинг, Г. Кнорр, Д. Монтгомери // Вычислительные методы в физике плазмы / под ред. Б. Ольдера, С. Фернбаха, М. Ротенберга. – М.: Мир, 1974. – C. 39–95.

[34] Armtsrong, T.P. Numerical studies of the nonlinear Vlasov equation / T.P. Armtsrong // Physics of Fluids. – 1967. – Vol. 10,– № 6. – P. 1269–1280.

[35] Knorr, G. Plasma simulation with few particles / G. Knorr // Journal of Computational Physics. –1973. – Vol. 13. – Iss. 2. – P. 165–180.

[36] Cheng, C.Z. The integration of the Vlasov equation in configuration space / C.Z. Cheng, G.

Knorr // Journal of Computational Physics. – 1976. – Vol. 22. – Iss. 3. – P. 330–351.



Pages:   || 2 |
Похожие работы:

«УДК 629.113 ББК-39.33 Х-86 А.Л. Хохлов Рабочая программа по дисциплине: "Ресурсосбережение при проведении технического обслуживания и ремонта" Димитровград, ТИ – филиал ФГОУ ВПО Ульяновская ГСХА, 2007. – 14 с. Рецензент: доцент кафедры Технические дисциплины к.т.н., Власов С.Н....»

«Т.Н. Зорина ДОКУМЕНТАЛЬНОЕ КИНО – ЖАНР КИНОИСКУССТВА И ИСТОРИЯ СТРАНЫ В контексте истории и событий Родившись как техническое изобретение братьев Люмьер, кинематограф быстро стал новым средством общения. Именно документальная съемка реальной жизни сделала его источником достоверной информации, не только поз...»

«ФЕДЕРАЛЬНОЕ АГЕНТСТВО ПО ОБРАЗОВАНИЮ ГОСУДАРСТВЕННОЕ ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ ВОЛГОГРАДСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ КАМЫШИНСКИЙ ТЕХНОЛОГИЧЕСКИЙ ИНСТИТУТ (ФИЛИАЛ) ВОЛГОГРАДСКОГО ГОСУДАРСТВЕННОГО ТЕХНИЧЕСКОГО УНИВЕРСИТЕТА О. Н. Иосифова ПРАКТИКУ...»

«ПЫЛЕСОС Руководство по эксплуатации VCC43A1 Благодарим за приобретение Пылесоса компании Midea. Перед началом эксплуатации пылесоса внимательно прочитайте данное Руководств...»

«МИНИСТЕРСТВО ОБРАЗОВАНИЯ РОССИИСКОЙ ФЕДЕРАЦИИ НОВОСИБИРСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ ИНСТИТУТ ДОПОЛНИТЕЛЬНОГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ РЕГИОНАЛЬНЫЙ ЦЕНТР ПЕРЕПОДГОТОВКИ ВОЕННОСЛУЖАЩИХ КАФЕДРА УЧЕТА И БАНКОВСКОГО ДЕЛА ФИНАНСОВЫЙ...»

«Министерство образования и науки Украины НАЦИОНАЛЬНЫЙ ТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ "ХАРЬКОВСКИЙ ПОЛИТЕХНИЧЕСКИЙ ИНСТИТУТ" В.О. Соловьев, И.М. Фык, Е.П. Варавина НЕТРАДИЦИОННЫЕ ИСТОЧНИКИ УГЛЕВОДОРОДОВ: ПРОБЛЕМЫ ИХ ОСВОЕНИЯ Учебное пособие Харьков НТУ "ХПИ" УДК ББК Н Рекомендовано к изданию решен...»

«Артур Джозеф Квеку Исследование и разработка модели и алгоритмов MACи физического уровней сетей WiMAX 05.12.13 Системы, сети и устройства телекоммуникаций Диссертация на соискание ученой степени кандидата технических наук Научный руководитель кандидат техническ...»

«ИСТОЧНИК ВТОРИЧНОГО ЭЛЕКТРОПИТАНИЯ РЕЗЕРВИРОВАННЫЙ СКАТ-1200Т исп.12/20 РУКОВОДСТВО ПО ЭКСПЛУАТАЦИИ ФИАШ.425519. 014 РЭ МЛ11 Благодарим Вас за выбор нашего источника резервного питания, который обеспечит Вам надежную работу систем сигнализации и связи на Вашем объекте. Настоящее руководство по эксплуатации предназначено...»

«Технический паспорт и инструкция по монтажу Насосная группа со смесителем MK. Поколение 8 ООО "Майбес РУС" 109129, Москва, Ул. 8-ая Текстильщиков, 11/2 Тел.:+7-495-727-20-26 Техника быстрого монтажа www.meibes.ru 1. Назначение изделия 1.1. Насосные группы MK Поколения 8 предназначены для подачи теплоносителя, поступающего...»

«EPG-111 Универсальный эпоксидный грунт Описание продукта. Двухкомпонентный универсальный эпоксидный грунт без растворителей, для исполнения полимерныхпокрытий по бетону. Компонент А – низковязкая модифицированная эпоксидная смола на основе бисфенолов A/F. Компонент В – модифицированный полиамин.Свойства...»

«КОД ОКП 42 2860 "УТВЕРЖДАЮ" Технический директор ЗАО "Радио и Микроэлектроника" С.П. Порватов Подп. и дата "_"_2008 г. Инв. № дубл. Взам. инв.№ Подп. и дата Инв. № подл Счетчики электрической энергии однофазные статические многотарифные СОЭБ-2П-65 СО...»

«МИНОБРНАУКИ РОССИИ Федеральное государственное бюджетное образовательное учреждение высшего профессионального образования "Ухтинский государственный технический университет" (УГТУ) ВОДОСНАБЖЕНИЕ И ВОДООТВЕДЕНИЕ ЖИЛОГО ЗДАНИЯ Методические указания к...»

«Утверждено постановлением Правительства Кыргызской Республики от "" 2015 года № _ Технические требования к моделям контрольно-кассовых машин, разрешенных к использованию на территории Кыргызской Республики 1.Общие положения 1.1. Настоящие Технические требования разработаны в целях внедрения и применения контрольно-кас...»

«Источник бесперебойного питания РУКОВОДСТВО ПОЛЬЗОВАТЕЛЯ BRICK 600/800 WWW.POWERMAN.RU Содержание Введение 1. 3 Инструкции по технике безопасности 2. 4 Принцип работы ИБП 3. 4 Установка 4. 5 Работа с ИБП 5. 7 Сигналы ИБП 6. 8 Батарея 7. 9 Хра...»

«ТЕОРИЯ И ПРАКТИКА ОБЩЕСТВЕННОГО РАЗВИТИЯ (2013, № 2) УДК 378 Рыжова Гульнара Алимовна Ryzhova Gulnara Alimovna старший преподаватель Senior Lecturer of кафедры иностранных языков the...»

«МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ Федеральное государственное бюджетное образовательное учреждение высшего профессионального образования "МОСКОВСКИЙ АВТОМОБИЛЬНО-ДОРОЖНЫЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ (МАДИ)" Северо-Кавказский филиал УТВЕРЖДЕНО Директор С-Кф МАДИ Г.Г. Ягу...»

«Федеральное агентство по образованию Сибирская государственная автомобильно-дорожная академия (СибАДИ) Кафедра эксплуатации и ремонта автомобилей ЭКОНОМИЧЕСКАЯ ОЦЕНКА ДЕЯТЕЛЬНОСТИ ПО ТЕХНИЧЕСКОМУ ОБСЛУЖИВАНИЮ И РЕМОНТУ ПОДВИЖНОГО СОСТАВА Методич...»

«СБОРНИК ДОКЛАДОВ И КАТАЛОГ УЧАСТНИКОВ ШЕСТОЙ МЕЖДУНАРОДНОЙ КОНФЕРЕНЦИИ "МЕТАЛЛУРГИЯ-ИНТЕХЭКО-2013" СБОРНИК ДОКЛАДОВ И КАТАЛОГ УЧАСТНИКОВ ШЕСТОЙ МЕЖДУНАРОДНОЙ КОНФЕРЕНЦИИ "МЕТАЛЛУРГИЯ-ИНТЕХЭКО-2013" СОДЕРЖАНИЕ СБОРНИКА ДОКЛАДОВ И КАТАЛОГА КОНФЕРЕНЦИИ Раздел №1. Список компа...»

«Пятницкий Илья Алексеевич Фотоиндуцированная агрегация продуктов фотоокисления псоралена и механизмы их иммуносупрессорного действия Диссертация на соискание ученой степени кандидата медицинских наук 03.01.02 – Биофизика 14.03.09 – клиническая иммунология, аллергология Научны...»

«МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РФ Федеральное государственное бюджетное образовательное учреждение высшего профессионального образования "Национальный исследовательский Томский политехнический университет" Энергетический институт МАТЕРИАЛЫ ДОКЛАДОВ "Энергосбережение – теория и пр...»

«Мобильные дробильные установки для окончательного измельчения ДРОБИЛЬНЬІЕ УСТАНОВКИ С КОНИЧЕСКОЙ ДРОБИЛЬНЬІЕ УСТАНОВКИ С РОТОРНОЙ ДРОБИЛКОЙ MOBICONE ДРОБИЛКОЙ MOBIFOX Мобильные дробильные установки для окончательного измельчения:...»








 
2017 www.lib.knigi-x.ru - «Бесплатная электронная библиотека - электронные матриалы»

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