Содержание к диссертации
Введение
1 Обзор методов численного решения системы уравнений Власова-Максвелла 10
1.1 Плазма, основные характеристики 10
1.2 Кинетическое описание бесстолкновительной плазмы 12
1.3 Методы решения уравнения Власова, основанные на восстановлении функции распределения f(t,x,v) 13
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 Схема метода частиц в ячейках 33
2.3 Схема Лакса - Вендроффа для уравнения Власова 36
2.3.1 Сравнение метода частиц и конечно-разностного метода 37
2.4 Алгоритм уменьшения счетных эффектов метода частиц в ячейках (алгоритм
вычитания шумовой добавки) 38
2.4.1 Алгоритм вычитания шума I 39
2.4.2 Алгоритм вычитания шума II 39
2.4.3 Алгоритм вычитания шума III 2.5 Корректировка положения частиц 41
2.6 Схемы, использованные в Алгоритме III 43
2.7 Результаты расчетов 44
2.7.1 Зависимость уровня шума от количества частиц 44
2.7.2 Сравнение схем для уравнений на v, п 47
3 Форма ядра частицы и проблема самовоздействия в методе частиц в ячейках 50
3.1 Самосила и VSP-ядро в одномерном случае 51
3.2 Самосила в двумерном случае 52
3.2.1 Самосила в двумерном случае, PIC-ядро 57
3.2.2 Самосила в двумерном случае, QCPl-ядро 57
3.2.3 Самосила в двумерном случае, С;СР2-ядро 58
3.2.4 Потенциал поля одиночного заряда 58
3.2.5 Сравнение ядер PIC, QCP1, QCP2 59
3.3.1 Новое ядро, тесты 61
3.4.1 Фурье-образ функции ядра частицы 66
3.4.2 Выбор оптимальных параметров 68
3.5 Саморазогрев модельной плазмы 72
3.6 Моделирование эволюции протопланетного диска с QCP-ядром 74
4 Число модельных частиц и точность на примере задачи моделирования кинетической неустойчивости теплого электронного пучка малой плотно сти в плазме 78
4.1 Основные уравнения 80
4.2 Методы и алгоритмы решения 81
4.3 Вычисление инкремента неустойчивости 84
4.4 Результаты расчетов 85
4.5 Фазовые плоскости 91
Заключение 96
Литература
- Методы решения уравнения Власова, основанные на восстановлении функции распределения f(t,x,v)
- Схема Лакса - Вендроффа для уравнения Власова
- Самосила в двумерном случае, QCPl-ядро
- Вычисление инкремента неустойчивости
Введение к работе
Актуальность работы. Метод частиц в ячейках, возникший еще в середине прошлого столетия, на сегодняшний день широко применяется при моделировании поведения плазмы. Его также используют и при расчете динамических процессов в других средах (жидкости, газы, сплошные среды и др.), но именно в решении задач физики плазмы он получил наиболее широкое распространение. Как самостоятельный раздел вычислительная физика плазмы сформировалась одновременно с развитием вычислительной техники во второй половине XX века. Это привело к выделению математического моделирования в отдельное направление исследовательской работы. Ввиду больших возможностей диагностики моделируемых в численных экспериментах кинетических процессов и относительно невысоких экономических затрат на их проведение (по сравнению со строительством установок для проведения реальных экспериментов), численные методы играют все более важную роль в решении задач физики плазмы.
Существует ряд математических моделей описания плазмы: кинетическое описание, магнитогидродинамическое приближение, модели термоядерной плазмы и другие. В 1938 году Власов предложил свою концепцию описания широкого круга процессов в плазме. В основе модели Власова лежит кинетическое уравнение без столкновительного члена для функции распределения частиц по скоростям fit, х, v). Электрическое и магнитное поля, входящие в уравнение через силу Лоренца, определяются из уравнений Максвелла. Система уравнений замыкается формулами для плотности зарядов и токов, которые выражаются через функции распределения частиц. Система уравнений Власова-Максвелла послужила основой для большого числа работ по теории колебаний, устойчивости, коллективных процессов в плазме.
В диссертации рассматривается приближение бесстолкновительной плазмы и метод частиц в ячейках применительно к системам уравнений Власова-Пуассона и Власова-Максвелла. Для многих нелинейных процессов найти решение этой системы возможно только численно, а не аналитически. Существует ряд методов численного решения уравнений Власова-Максвелла для бесстолкновительдной плазмы: метод преобразований, метод водяного мешка, конечно-разностные методы, методы частиц и другие.
Метод частиц в ячейках принадлежит группе вычислительных алгоритмов, объединенных способом дискретизации, которая носит название «методы частиц». В отличие от конечно-разностных методов, объединенных аппроксимацией дифференциальных или интегральных операторов исходных уравнений, методы частиц объединяет концепция представления среды в виде множества модельных частиц, которые являются носителями набора ха-
4 рактеристик этой среды (масса га, заряд q, импульс p и т.п.). Множество точек, которое составляют координаты частиц, принято называть «лагран-жевой» сеткой, в противовес «эйлеровой» сетке, привязанной к области. В этой терминологии методы частиц можно разделить на лагранжевы методы и смешанные эйлерово-лагранжевы методы. Метод частиц в ячейках относится к группе смешанных алгоритмов, вычисления в которых проводятся в два этапа: на эйлеровой и на лагранжевой сетке попеременно.
Такие преимущества, как отсутствие аппроксимационной вязкости и возможность осуществления сквозного счета при резких изменениях поведения решения в областях больших градиентов, поставили методы частиц для задач с большими объемными деформациями и неустойчивостями на первое место. Вследствие этого методы частиц нашли широкое применение именно для моделирования таких сложных явлений, как неустойчивости в физике плазмы. Главным преимуществом метода частиц в ячейках по сравнению с эйлеровыми методами является простота реализации и экономичность, так как в методах частиц не вычисляется напрямую функция распределения. Простота реализации данного метода (в том числе и на параллельных ЭВМ), наглядность получаемых распределений частиц, большое число решенных на его основе задач составляют сильные стороны этого метода.
Основным недостатком метода частиц в ячейках всегда являлось наличие в решении нефизических эффектов, так называемых «численных шумов». Эта проблема очень сложна, так как причин возникновения шумов много. Влияние различных факторов на решение трудно разделить, особенно при моделировании неустойчивых сред. Шумовые гармоники взаимодействуют друг с другом, накладываются на гармоники неустойчивых решений, что может приводить к развитию численных неустойчивостей.
Самым простым методом уменьшения численных шумов является увеличение количества частиц до близкого к реальным значениям, что из-за ограниченности ресурсов ЭВМ не представляется возможным. Для уменьшения численных эффектов подбирают более гладкую форму ядра модельной частицы, вычисляют оптимальный шаг по времени и по пространству, используют различные модификации метода, исходя из специфики поставленной задачи. Используются также алгоритмы сглаживания, но они не устраняют сам шум и могут искажать физические эффекты. На данный момент единого подхода к решению проблемы численных шумов нет. Поэтому создание алгоритмов метода частиц с пониженным уровнем шума остается актуальной проблемой вычислительной математики и математического моделирования.
Цель диссертационной работы заключается в исследовании свойств метода частиц в ячейках и разработке общих алгоритмов уменьшения счет-
5 ных шумов данного метода. Для достижения этой цели было решено
создать алгоритм подавления шума в методе частиц в ячейках,
разработать методику выбора оптимальной формы ядра, которая обеспечивает по сравнению со стандартным PIC-ядром уменьшение самосилы,
создать комплекс программ для моделирования взаимодействия электронного пучка с плазмой в трехмерной геометрии и провести оценки точности получаемого решения.
Научная новизна работы заключается в следующем:
-
Разработан новый алгоритм уменьшения счетных эффектов метода частиц в ячейках, основанный на вычитании шумовой добавки, обеспечивающий подавление счетных флуктуаций при сохранении самого решения. Данный алгоритм реализован при решении задачи о распаде разрыва плотности ионов в неизотермической плазме.
-
Проведено исследование причины возникновения самосилы для метода частиц в ячейках на смещенных сетках. Для уменьшения этой ошибки разработан экономичный подход, основанный на использовании нового ядра модельной частицы с меньшей, по сравнению с PIC-ядром, самосилой.
-
Показано, что использование нового ядра приводит к увеличению времени саморазогрева модельной плазмы и, как следствие, к лучшему сохранению импульса и полной энергии.
-
На основе метода частиц в ячейках создан комплекс программ для моделирования взаимодействия электронного пучка с плазмой в трехмерной постановке при параметрах, соответствующих условиям лабораторных экспериментов на установке ГОЛ-3, который позволяет проводить вычисления как в гидродинамическом, так и в кинетическом режимах развития пучковой неустойчивости.
-
Проведены оценки величины погрешности решения, получаемого данным комплексом программ. Для различных режимов приводятся оценки достаточного количества частиц.
Научная и практическая ценность работы состоит в создании алгоритмов подавления шума в методе частиц в ячейках, а также в создании трехмерной программы моделирования взаимодействия пучка электронов малой плотности с плазмой, охватывающей три режима развития неустойчивости. Созданные алгоритмы могут быть использованы при решении различных за-
дач физики плазмы методом частиц для уменьшения уровня нефизических флуктуаций.
Представленные в диссертации исследования проводились в рамках Интеграционных проектов СО РАН №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 Всероссийской молодежной конференции-школе «Современные проблемы математического моделирования. Математическое моделирование, чис-
7 ленные методы и комплексы программ» (2009, Дюрсо), на VI Всесибирском конгрессе женщин-математиков (2010, Красноярск), на XV Байкальской Всероссийской конференции «Информационные и математические технологии в науке и управлении» (2010, Иркутск-Байкал), на Всероссийской конференции по вычислительной математике КВМ-2011 (Новосибирск), на Российско-монгольской конференции молодых ученых по математическому моделированию, вычислительно-информационным технологиям и управлению (2011, Иркутск (Россия) - Ханх (Монголия)), на XIX Всероссийской конференции «Теоретические основы и конструирование численных алгоритмов и решение задач математической физики» (2012, Дюрсо), на XIII Всероссийской конференции молодых ученых по математическому моделированию и информационным технологиям (2012, Новосибирск), на Международной конференции «Mathematical modeling and computational physics» (2013, Дубна).
Основные результаты опубликованы в 15 работах, из которых 2 в журналах, рекомендованных ВАК. Список публикаций приведен в конце автореферата.
Личный вклад соискателя заключается в обсуждении постановок задач, разработке численных алгоритмов и методов решения, создании и тестировании программ, проведении расчетов и анализе полученных результатов. Все выносимые на защиту результаты принадлежат лично автору. Представление результатов совместных исследований и разработок согласовано с соавторами.
Структура и объем работы. Диссертация состоит из введения, 4 глав и заключения. Список использованной литературы содержит 133 наименования (включая 15 публикаций автора). Текст диссертации содержит 110 страниц машинописного текста, включая 36 рисунков и 4 таблицы.
Методы решения уравнения Власова, основанные на восстановлении функции распределения f(t,x,v)
Исторически методы частиц возникли и развивались как альтернатива сеточным алгоритмам. Данная группа алгоритмов, получившая общее название Эйлеровы алгоритмы (Eulerian solvers), объединяет в себе методы вычисления значений функции распределения /(, x,v) для уравнений Власова-Пуассона или Власова-Максвелла на связанной с областью пространственной сетке. Сюда относятся следующие методы: конечно-разностные методы (fnite-diference method), метод конечных объемов (fnite volume method) или метод баланса потоков (fux 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). В силу своей сложности данный метод используется как правило только в одномерной постановке.
Методы частиц - это особая группа вычислительных алгоритмов, которые объединяет между собой способ дискретизации исходной математической задачи. В этом методе среда (в нашем случае плазма) представляется набором достаточно большого числа модельных частиц. Движение частиц подчиняется законам классической механики в самосогласованном электромагнитном поле, определяемом из уравнений Максвелла. Методы частиц широко применяются при решении задач не только физики плазмы, но и динамики жидкостей и газов, механики сплошной среды и т.д. Но самое широкое применение они получили именно в физике плазмы. Большим плюсом этих методов по сравнению с приведенными выше, является экономичность, т.к. в методах частиц не вычисляется напрямую функция распределения f(t,x,v). Особенно заметным это преимущество становится при переходе от одномерных постановок к задачам большей размерности.
Данные методы с самого начала были ориентированы на решение задач с сильными деформациями, сдвигами, неустойчивостями на поверхностях раздела сред, где конечно-разностные схемы оказались непригодны, поэтому они нашли широкое применение именно для таких сложных задач, как моделирование неустойчивостей в физике плазмы. Развитая теория, простота реализации этого метода (в том числе и на параллельных ЭВМ), физическая наглядность получаемых распределений частиц, большое число решенных на его основе задач в разнообразных областях физики плазмы составляют сильные стороны этого метода.
В вычислительной физике плазмы методы частиц начали использоваться с конца 1950-х годов. Первый одномерный электростатический вариант - модель плоских листов - был предложен О. Бунеманом в 1959 г [48] и развит позже в работах Дж. Доусона [49].
Схема Лакса - Вендроффа для уравнения Власова
Проблема точности метода частиц в ячейках особенно актуальна в задачах физики плазмы. Недостаточное число частиц как в гидродинамических кодах, так и при моделировании плазмы приводит к флуктуациям сеточных величин. Однако, если в первом случае это как правило всего лишь усиливает зашумленность результатов, то в случае плазмы, за счет присущих этому состоянию вещества дальнодействующих взаимодействий, флуктуации в плотности заряда могут вызывать развитие новых неустойчивостей, соответствующих физике модели, но не соответствующих реальному исследуемому физическому процессу.
Совокупность проявлений вычислительных ошибок метода частиц принято называть «численным шумом» - термин, заимствованный из теории сигналов. Шумовые гармоники накладываются на гармоники решения и при больших временах бывает трудно судить о решении. Источников нефизических шумов в таких сложных алгоритмах, как метод частиц в ячейках, много: наличие пространственной сетки (на которой определяются значения электромагнитных полей, распределения плотности, скорости, тока), разбиение среды на модельные частицы (изменяется потенциал взаимодействия между частицами на масштабе сеточной ячейки), переходы с частиц на сетку и обратно, дискретный шаг по времени. Все эти погрешности приводят к нарушениям в численных моделях законов сохранения энергии, импульса, изменяются дисперсионные свойства, возрастает число столкновений, происходит стохастический саморазогрев ансамбля частиц, который ограничивает допустимый интервал, отслеживаемый в расчетах.
К вопросу количественной оценки уровня нефизических флуктуаций подходили с разных сторон [79,80]. Подробное описание ряда причин шумов и методов их устранения можно найти в работах [51,52,81,82]. На настоящий момент единого подхода к решению этой проблемы нет. Для уменьшения численных шумов увеличивают число модельных частиц, меняют их форму, выбирают оптимальный по времени и пространству шаг, меняют начальное распределение частиц, используют различные модификации метода исходя из специфики поставленной задачи и др. Предлагаются также алгоритмы сглаживания, но они не устраняют сам шум и могут сглаживать физические эффекты, что нежелательно.
Количество частиц
При моделировании бесстолкновительной плазмы каждая модельная частица соответствует десяткам миллионов физических частиц. Обозначим штрихом величины, относящиеся к модельной плазме. Из условий сохранения заряда и массы следует, что е = ае, т! = am, где а = р/р - отношение плотностей реальной и модельной плазмы. Тогда, как показано в [16,18], следующие характеристики реальной плазмы сохраняются и в модельной не смотря на сравнительно малое число модельных частиц. Это собственная плазменная частота шре,г = »ег, дебаевский радиус экранирования (в предположении равенства реальной и модельной тепловой скорости v T = VT) XD = X D. Также сохраняется отношение заряда частицы к массе е/т = е /га , а значит динамика частиц не меняется.
В то же время для ряда величин имеем следующие соотношения. Температура Т" = аТ, длина свободного пробега А = Х/а, число частиц в дебаевской сфере N = N/a, частота столкновений v1 = и а. Поскольку а 1, то длина свободного пробега частицы в модельной плазме уменьшается, частота столкновений становится больше, уровень тепловых флуктуаций превышает реальный.
Увеличить количества частиц с целью как можно ближе подойти к реальным значениям приведенных параметров - самый простой способ уменьшения численных шумов. Но из-за ограниченности ресурсов ЭВМ это не всегда представляется возможным. Тем не менее во многих задачах запас «бесстокновительности» достаточно велик (А 3 L) и неравенство А 3 L (L - размер области) для модельной плазмы можно сохранить. При выборе числа частиц (или же при выборе размера области) необходимо также учитывать, что число частиц в дебаевской сфере не должно быть маленьким N D 1.
Адаптивное изменение количества частиц
Следует отметить, что проблема шумов не теряет своей остроты при повышении числа модельных частиц в расчетной области даже при расчетах на суперЭВМ, так как с изменением плотности вещества могут появиться ячейки с небольшим количеством частиц, и общий уровень шумов в расчете будет определяться именно этими ячейками. С другой стороны, в наиболее заполненных ячейках при этом может быть частиц в десятки раз больше необходимого уровня, что создает избыточную нагрузку на вычислительную систему, таким образом, повышая время расчета. Для преодоления этих проблемм создаются методы частиц в ячейках с адаптивным изменением числа частиц [71,83-88]. Кроме того создан метод частиц в ячейках с частицами переменных размеров [89].
Форма модельной частицы
Соответствие модельных частиц их реальным прототипам далеко от полного совпадения не только количественно. Модельные частицы конечного размера принято интерпретировать как облака частиц. Но форма частиц и распределение признака в них в процессе моделирования остаются неизменными, что может противоречить реальной ситуации. Описанию свойств ядер и того, какие погрешности вносит такой способ дискретизации, посвящено большое число работ различных авторов. Большинство основных выводов и ссылки на первые работы можно найти в следующих монографиях и статьях [18,51,52,77,80,81,90].
В работе [91] приведено сравнение сохранения полной энергии для NGP и PIC ядер, в [92] для параболического и PIC ядер. С увеличением числа частиц в шаблоне и повышением уровня гладкости полная энергия системы сохраняется лучше, флуктуации плотности и поля уменьшаются. В работе [93] получены погрешности аппроксимации признаков с лагранжевой на эйлерову сетку для NGP и PIC ядер. Показано, что при использовании NGP ядра сеточная функция плотности рд аппроксимирует исходную плотность р(х) со вторым порядком по h, а отклонение 8д = \рд — р(хд)\ обратно пропорционально среднему числу частиц в ячейке. При использовании PIC ядра рд аппроксимирует исходную плотность р(х) также со вторым порядком по fa,, а отклонение 8д = \рд — р(хд)\ обратно пропорционально квадрату среднего числа частиц в ячейке.
Одна из возможных модификаций метода частиц, связанная с ядрами, это использование разных ядер при раздаче заряда частицы и интерполяции сил, действующих на частицу. В [51] на примере показано, что использование схемы PIC/NGP вносит сравнительно неопасные колебания, в то время как NGP/PIC в некоторых может вести к экспоненциальному росту неустойчивости.
Самосила в двумерном случае, QCPl-ядро
Приведем результаты расчетов с использованием описанного Алгоритма III для вычитания шумовой добавки. Красным цветом на всех рисунках изображены графики средней скорости и плотности, которые были получены конечно-разностным методом для системы уравнений Влпсова-Пуассона (схема Лакса-Вендроффа).
Параметры метода частиц: число частиц т = 40000, число узлов сетки п = 100 (число частиц в ячейке 1р = 400). Параметры конечно-разностной схемы: число узлов сетки по ж и и равно п = 1600 и km = 3000 соответственно. В этих тестах для уравнений (2.23), (2.24) использовалась противопотоковая схема, время выбрано небольшим tm = 0.1. На рисунке 2.6-2.9 cиним изображены средняя скорость и плотность, полученные стандартным методом частиц в ячейках без вычитания шумовой добавки, зеленым - средняя скорость и плотность, полученные по методу частиц в ячейках с Алгоритмом III вычитания шумовой добавки. Из графиков видно, что решение, полученное по методу частиц и методу частиц с вычитанием шумовой добавки ведет себя так же, как и решение, полученное по схеме Лакса-Вендроффа для уравнения Власова. Видно также, что с вычитанием шумовой добавки графики более гладкие, особенно это заметно на тесте с меньшим числом частиц.
Сравнение с решением по схеме Лакса-Вендроффа проводится как с точным решением. Но аналитически как должно вести себя решение точно известно только в областях A и B, куда еще не дошло возмущение. Обозначим эти области А = [0,1], В = [3, 4] и С = A В. Средняя скорость в этих областях должна быть равна VA = VB = 0, плотность ПА = 2, пв = 1 . Поэтому можно легко оценить влияние алгоритма вычитания шума в данных областях.
Чтобы численно оценить полученные результаты, рассмотрим насколько отличаются от точного решения в областях А и В графики плотности и средней скорости для обычного и модифицированного метода частиц. Введем суммарное отклонение от точного решения по области С для средней скорости и плотности следующим образом
В таблице 2.1 индексом ріс помечены величины, относящиеся к опо обычному (PIC без вычитания) методу частиц в ячейках, индексом new - относящиеся к модифицированному (PIC с вычитанием шума) методу (метод частиц в ячейках с Алгоритмом вычитания шумовой добавки III). Здесь, как и раньше, 1р - число частиц в ячейке. Число узлов сетки п = 100, время t = 0.1, г = 0.001. Из таблицы 2.1 видно, что для метода частиц с вычитанием шумовой добавки величина уровня шума на два-три порядка меньше, чем для обычного метода частиц как в плотности, так и в средней скорости. Таблица 2.1: Зависимость уровня шума от количества частиц для метода частиц в ячейках (величины с индексом pic) и для метод частиц в ячейках с Алгоритмом III (индекс new).
Необходимо отметить, что на решение в центральной области влияет выбор схемы для нахождения v, п. На рисунках 2.11 и 2.10 приведены результаты расчетов с использованием стандартного метода частиц без вычитания шума, противопотоковой схемы (2.34)-(2.38) и схемы предиктор-корректор (2.39)-(2.42) для параметров tm = 0.3, г = 0.001. Параметры метода частиц: число частиц т = 10000, число узлов сетки п = 101. Параметры конечно-разностной схемы: число узлов сетки по x и u равно im = 1600 и km = 3000 соответственно.
Зависимость качества решения от выбранной схемы. а - плотность, б - средняя скорость. Красным цветом - решение, полученное по схеме Лакса-Вендроффа для уравнения Власова. Уменьшить уровень шума в областях A и B позволяет алгоритм с использованием любой из рассмотренных схем. Но из этих графиков видно, что использование схемы предиктор-корректор (рисунки 2.11 и 2.10, черным цветом) меньше влияет на поведение решения в области распространения возмущений, чем противопотоковая схема. Использование проти-вопотоковой схемы подавляет передний пик в графиках плотности и средней скорости (рисунки 2.11 и 2.10, зеленым цветом). Это связано с тем, что противопотоковая схема имеет
Рис. 2.11: а - плотность, б - средняя скорость для различных схем. Красным цветом - решение, полученное по схеме Лакса-Вендроффа для уравнения Власова. первый порядок точности, а в методе частиц в ячейках используется грубая пространственная сетка.
В данной главе на примере задачи распада разрыва плотности ионов предложен алгоритм уменьшения счетных шумов метода частиц. Предложенный алгоритм был реализован в программе на языке Fortran для решения одномерной задачи распада разрыва плотности ионов. Проведенное исследование метода показало его перспективность. Алгоритм вычитания шумовой добавки позволяет существенно уменьшить разброс в средней скорости и плотности в методе частиц в ячейках для поставленной выше задачи. Качество решения зависит от выбора подходящей схемы для уравнений нахождения v, п Алгоритма вычитания шума III.
В предыдущей главе численный шум в методе частиц в ячейках рассматривался как общее понятие. В этой главе рассматривается одна из причин, приводящих к изменению движения частицы, и как следствие, к повышению столкновительности в плазме, а именно наличие пространственной сетки и возникновение так называемой «самосилы». Увеличение количества модельных частиц ведет к уменьшению шумовых эффектов в модельной плазме, но полностью их не устраняет. Одним из эффектов, которые нельзя устранить увеличением количества модельных частиц, является сила, с которой частица действует сама на себя через сетку ("самовоздействие"или "самосила").
Вычисление инкремента неустойчивости
На графиках, представленных ниже, приведено сравнение результатов расчетов инкремента неустойчивости по энергии электрического поля (4.12) и по амплитуде главной моды (4.13) с решением, описанным в [132] (для одномерной задачи). Расчеты проводились при параметрах, соответствующих трем разным режимам развития пучковой неустойчивости. Счетные параметры: сетка по пространству 100 х 4 х 4, шаг по времени г = 0.001. Обозначим 1р число частиц в одной ячейке сетки и рассмотрим вопрос точности получаемого решения в разных режимах в зависимости от 1р, а также вопрос ранней идентификации неустойчивого решения.
1) Гидродинамический реэюим (kAv 7) Основные параметры: отношение плотности пучка к плотности плазмы щ/по = 2-10_3, разброс электронов пучка по скоростям AV/VQ = 0.035, длина области L = 1.2566. При заданных параметрах величина инкремента неустойчивости должна быть 7 = 0.0722 [132].
На рисунке 4.2 показан темп роста энергии электрического поля и значение инкремента 7i, вычисленное по формуле (4.12). Здесь и далее на графиках с индексом а прямой пунктирной линией показан темп роста волны с аналитическим инкрементом, положение этой линии снесено по координате t, максимумы графиков совмещены по времени. На графиках с индексом б прямой пунктирной линией показано точное значение инкремента. На графиках видна одна из особенностей моделирования холодной плазмы методом частиц - саморазогрев модельной плазмы, который выражается в том, что при малых температурах плазмы в самом начале вычислений энергия возрастает до определенного значения [51], величина которого зависит от счетных параметров. Как видно из рисунка 4.2, уровень саморазогрева уменьшается с увеличением числа частиц. Когда амплитуда развивающейся неустойчивой волны становится больше уровня шумов, связанных с саморазогревом, экспоненциальная стадия роста волны становится видна на графиках (линейный рост кривой на рисунке 4.2 а). Далее Рис. 4.2: Временная зависимость логарифма энергии электрического поля (а) и производная логарифма энергии электрического поля (инкремент 1) (б) при разном числе частиц в ячейке. Гидродинамический режим.
Временная зависимость логарифма амплитуды электрического поля (а) и производная логарифма амплитуды (инкремент 2) (б) при разном числе частиц в ячейке. Гидродинамический режим. решение выходит из линейного режима, и темп роста энергии спадает. Наступает стадия насыщения неустойчивости. Мы рассматриваем только линейную стадию развития неустойчивости, на которой и вычисляем значение инкремента. На графике, расположенном справа на рисунке 4.2, хорошо видна сходимость по количеству частиц к требуемому значению в линейной стадии.
Второй способ вычисления инкремента (по амплитуде главной волны электрического поля, формула (4.11)) показан на рисунке 4.3. Здесь, как видно из графиков, не играют роли шумы от саморазогрева модельной плазмы, что позволяет практически сразу наблюдать экспоненциальный рост. Более того, как видно из представленного графика, для той же точности вычисления инкремента можно использовать меньшее количество частиц.
В таблице 4.1 приведено значение инкремента, вычисленного по энергии и по амплитуде, а также характерная величина ошибки в процентах. В гидродинамическом случае ошибка по энергии примерно в 2 раза больше, чем по амплитуде. Из приведенных графиков, а также из таблицы видно, что в гидродинамическом режиме достаточным можно считать число частиц в ячейке, равное 25-50, если использовать вычисление инкремента по амплитуде. Ошибка в этом случае не превышает 10%. 2) Переходный режим (kv 7)
Основные параметры: отношение плотности пучка к плотности плазмы щ/по = 2-10_3, разброс электронов пучка по скоростям v/vo = 0.14, длина области L = 1.1424. При этих параметрах величина инкремента неустойчивости должна быть 7 = 0.0232. Рис. 4.4: Временная зависимость логарифма энергии электрического поля (а) и производная логарифма энергии электрического поля (инкремент 1) (б) при разном числе частиц в ячейке. Переходный режим.
Временная зависимость логарифма энергии электрического поля (а) и производная логарифма энергии электрического поля (инкремент 2) (б) при разном числе частиц в ячейке. Переходный режим.
Величина инкремента неустойчивости по сравнению с гидродинамическим режимом уменьшилась в три раза. Таким образом энергия нарастает медленнее, насыщение неустойчивости происходит позже по времени. Как и в гидродинамическом режиме, на графике логарифма энергии электрического поля (рисунок 4.4) виден саморазогрев модельной плаз 89 мы. На графике амплитуды основной моды (рисунок 4.5) этого эффекта нет. Из приведенных рисунков и таблицы 4.2 видно, что в переходном режиме необходимо использовать гораздо большее число частиц в ячейке по сравнению с гидродинамическим режимом. При увеличении числа частиц значение инкремента приближается к точному. При недостаточном числе частиц (1р = 50) график инкремента 72 (рисунок 4.5 б) не имеет области, близкой к горизонтальной прямой, которая была четко видна в гидродинамическом расчете и которая появляется в переходном режиме при большем числе частиц. При 1р = 250 область экспоненциального роста видна лучше. Далее с ростом числа частиц кривая \п(Е0)/2 в интервале [0,150] приближается к прямой линии, а 72(t) выходит на уровень, примерно соответствующий точному значению инкремента. Таким образом в данном случае для того, чтобы вычислить инкремент, нам необходимо использовать для расчетов как минимум 250 частиц в ячейке.
По ресурсам вычисления становятся более затратными, считать приходится на большие времена с большим числом частиц, чем в первом расчете. Из таблицы 4.2 видно, что ошибка по энергии также в 1.5-2 раза больше, чем по амплитуде. Для того, чтобы ошибка вычисления инкремента по амплитуде была в пределах 10%, число частиц в одной ячейке должно быть 1р « 2500. Чтобы достичь той же точности при вычислении инкремента по энергии электрического поля, потребуется 5000 частиц в одной ячейке.
3) Кинетический реэюими (kAv гу).
Основные параметры: отношение плотности пучка к плотности плазмы щ/по = 2- Ю-4, разброс электронов пучка по скоростям Av/vo = 0.14, длина области L = 1.1424. Аналитическое значение инкремента неустойчивости в этом случае 7 = 0.0027. Рис. 4.6: Временная зависимость логарифма энергии электрического поля (а) и производная логарифма энергии электрического поля (инкремент 72) (б) при разном числе частиц в ячейке. Кинетический режим.
Для вычисления инкремента в кинетическом режиме требуется гораздо больше частиц, чем в переходном. Далее, при рассмотрении картины поведения электронов пучка на фазовой плоскости станет понятно, почему это так. Величина инкремента по сравнению с переходным режимом уменьшилась почти в 9 раз, время насыщения неустойчивости увеличилось примерно в 6 раз, из-за численных шумов вычисление инкремента по энергии электрического поля стало невозможным. Поэтому здесь показан только график амплитуды электрического поля (рисунок 4.6). Несмотря на перечисленные сложности, вычисление инкремента по амплитуде дает удовлетворительный результат и показывает, что развивается именно та неустойчивая волна, которая соответствует выбранным начальным и граничным значениям.
Как и в переходном режиме, при недостаточном числе частиц график инкремента 72 при 1р = 500 не имеет области, близкой к горизонтальной прямой, и, соответственно, выбрать точку для вычисления инкремента невозможно. С ростом числа частиц область экспоненциального роста становится шире. Для вычисления значения инкремента в кинетическом режиме с ошибкой в 25%, требуется около 5000 частиц в ячейке. Отметим, что, в отличие от гидродинамического режима, графики в переходном и кинетическом режиме приведены со сглаживанием.
