CC BY
80
УДК 550.388.2
Л. В. Зинин, С. А. Ишанов, А. А. Шарамет, С. В. Мациевский
МОДЕЛИРОВАНИЕ РАСПРЕДЕЛЕНИЯ ИОНОВ ВБЛИЗИ ЗАРЯЖЕННОГО СПУТНИКА МЕТОДОМ МОЛЕКУЛЯРНОЙ ДИНАМИКИ. 2-D ПРИБЛИЖЕНИЕ
Рассмотрена математическая модель взаимодействия разреженной тепловой плазмы с положительно заряженным космическим аппаратом упрощенной формы. Использовался метод молекулярной динамики. Результаты моделирования показывают наличие скопления протонов в области перед спутником и ионной тени сразу за ним.
Mathematical model of interaction between low density thermal plasma with charged simple shape space probe was considered. Molecular dynamics simulation method was used. Modeling result shows a positive cloud of protons before the satellite and ion shadow behind it.
Ключевые слова: математическое моделирование, заряженный спутник, методы молекулярной динамики, параллельное программирование.
Kew words: mathematical modeling, charged satellite, method of molecular dynamics, parallel programming.
Введение
При постановке и интерпретации экспериментов на космических аппаратах необходимо проводить анализ распределения электрического поля вблизи заряженного спутника ввиду того, что сравнительно небольшой потенциал спутника способен внести существенные искажения в функцию распределения частиц. Потенциал в несколько вольт приводит к тому, что масс-спектрометрические измерения легких ионов, таких как H+, становятся невозможны. Даже в случае небольших потенциалов спутника для анализа экспериментальных измерений необходимо знать распределение электрического поля вблизи спутника, чтобы оценить его влияние на траектории частиц, степень и характер искажения функции распределения ионов.
Например, для бортовых измерений тепловых ионов с высотного спутника, потенциал которого соизмерим или превышает kT измеряемых ионов (где k — постоянная Больцмана; Т — температура), происходит значительное искажение величин измеряемых потоков ионов. Радиус Дебая RD также играет важную роль характерного пространственного размера взаимодействия в плазме. Особенно усложняется задача, когда RD сравним с характерными размерами спутника.
53
Вестник Балтийского федерального университета им. И. Канта. 2012. Вып. 10. С. 53 — 60.
На данный момент существуют несколько подходов, применяемых для расчета пространственного распределения электрического поля вблизи космического аппарата для различных ситуаций в окружающей плазме.
Первая задача, которая в большой степени уже принципиально решена [1-4], заключается в нахождении распределения электрического поля вокруг тела сложной формы в приближении Rd = œ. Для решения этой задачи требуется задать потенциал на поверхности спутника S либо знать распределение зарядов a(S). Обычно для нахождения распределения потенциала в вакууме решают уравнение Лапласа ДФ = 0 с известными граничными условиями: значением потенциала на бесконечности (равного 0) и на поверхности заряженного тела.
Данный подход требует использования граничных условий по потенциалу спутника, а не по плотности зарядов ct(S) на его поверхности. При таком подходе могут быть решены наиболее часто встречающиеся задачи, например случай эквипотенциального спутника. Вместе с тем для решения уравнения Лапласа в окрестностях тела сложной формы достаточно трудоемким оказывается описание границ, то есть создание математической модели космического аппарата. Отметим, что такой подход применим для предельного случая RD = œ и не дает реального распределения частиц вокруг спутника. Такая модель успешно использовалась для анализа измерений прибора Гиперболоид на спутнике Интербол-2 [2; 5; 6].
Вторым предельным случаем является так называемое приближение тонкого слоя, когда RD ^ 0. В случае приближения тонкого слоя главными факторами становятся направление и величина относительной скорости движения плазмы (спутника), значения тепловых скоростей частиц и спутниковый потенциал.
Задача с конечным радиусом Дебая наиболее сложная и актуальная. Прогресс в ее решении связан как с развитием теории численных методов [7; 8], так и со значительно возросшими мощностями компьютеров, позволившими рассчитать параметры плазмы. Затруднения возникают в связи с многомерностью задачи и требованиями высокого временного разрешения.
Предпринимается ряд попыток решить ее в самой упрощенной геометрической интерпретации. Так, для моделирования взаимодействия тела и плазмы применялись, в частности, следующие два подхода.
1. Гидродинамический подход, рассмотренный в работах [9 — 12], может быть использован только для высокой концентрации частиц свыше 105 см-3 и не позволяет рассчитывать траектории реальных частиц. Радиус Дебая в такой постановке является характерным масштабом области интегрирования, а плазма рассматривается как заряженная жидкость.
2. Метод крупных частиц (Particle In Cell или PIC метод), наиболее часто используемый в настоящее время, применяется для более разряженной плазмы и дает вполне удовлетворительные результаты, описанные в работах [10—12].
Класс аналитических решений расчета пространственного распределения электрического поля и параметров плазмы вокруг спутника с известным потенциалом приводится в классической книге [13].
Метод молекулярной динамики
Компьютерная молекулярная динамика (МД) — это один из наиболее мощных вычислительных методов, эффективно применяемых для моделирования физических систем [14]. Метод МД позволяет вычислять траектории отдельных атомов и ионов, исследовать динамику взаимодействия частиц на молекулярном уровне [15 — 18].
МД обладает также высоким пространственным и временным разрешением и позволяет получить информацию о процессах, происходящих в атомно-молекулярных масштабах и на временах порядка нескольких наносекунд [19 — 20]. Этот метод особенно эффективен на таких масштабах, где квантовые эффекты менее существенны, чем электростатические взаимодействия.
Современное развитие мощной вычислительной техники дает возможность моделировать динамику молекулярных систем, состоящих из огромного числа (от десятков тысяч до миллионов), с большим набором параметров и разнообразных условий, имитирующих реальный эксперимент. Существенный плюс метода МД в его применимости для любой концентрации.
При моделировании производится непосредственный учет парного взаимодействия отдельных атомов, или, как в нашем случае, протонов и электронов, что позволяет получить точные результаты. Главный фактор, ограничивающий применение метода МД, — его высокая потребность в вычислительных ресурсах, поэтому для его использования требуются суперкомпьютерные мощности.
Основные положения
Для описания движения атомов и частиц применяется классическая механика.
Силы взаимодействия ион-электрон можно представить в форме классического кулоновского взаимодействия.
Точное знание траектории движения частицы системы на больших промежутках времени не является необходимым для получения результатов макроскопического характера.
Начальное распределение частиц соответствует определённой статистической функции распределения.
Метод молекулярной динамики применим, в случае когда длина волны Де Бройля атома или частицы много меньше, чем межатомное расстояние. Кроме того, при низких температурах квантовые эффекты становятся определяющими, и для рассмотрения таких систем необходимо использовать квантово-химические методы.
Математическая модель
В описании модели движения частиц используются уравнения классической механики:
56
аХі2
(1)
У = У о + Ъу +-ур (2)
где х0 у0 — начальные координатні частицы; ух, ъу — проекции скорости; Ох, ау — проекции ускорения на оси абсцисс и ординат соответственно.
Для описания взаимодействия используется закон Ньютона с силой Лоренца в правой части:
та = цЁ + ц[У х Н], (3)
где т, а, ц, V — масса, ускорение, заряд и скорость частицы соответственно; ё, Н — вектора напряженности электрического и магнитного
полей соответственно.
Геометрия модели
Заменим реальную ситуацию взаимодействия плазмы с движущимся в ней заряженным телом (спутником) следующей модельной ситуацией. Рассмотрим плоский случай, то есть заряженное тело представим кругом, который является проводящим и на границе имеет положительный потенциал Ф0. Спутник движется в невозмущённой плазме с постоянной скоростью и0 (рис. 1). Далее рассмотрим случай, в котором магнитное поле не учитывалось.
X,
Фо
4-' [ V
1 \
V
-6 і э-
Рис. 1. Геометрия модельной задачи
Математическое моделирование осуществим в декартовой системе координат ХУ. Рассмотрим квадрат со стороной I, одна из вершин которого имеет координаты (0, 0), а два ребра лежат на координатных осях. Этот квадрат — окрестность для расположенного в центре заряженного тела-круга и выбирается так, чтобы можно было считать влияние заряженного тела вне этой окрестности пренебрежимо малым. Заряженное тело-круг задается следующим образом: центром в точке (Хс, Ус) и радиусом Я. На границе тела поддерживается постоянный положительный потенциал Ф0. Круг движется со скоростью Ы0 вдоль оси абсцисс. Электрический потенциал невозмущенной плазмы равен 0.
Описание алгоритмов
В ходе моделирования для пересчета положения частицы на новый временной слой применяются (1 — 2), и при отсутствии магнитного поля формула (3) примет вид та = уБ. Исходя из этого получим уравнения для пересчета скорости частицы для нового слоя по времени:
m
= qE, (4)
At
\ /
k — номер временного слоя; E — напряженность электрического поля;
E = tE, (5)
i=1
n — количество частиц, находящихся на расстоянии радиуса Дебая, или меньше (относительно рассчитываемой частицы);
E = *' (6)
ri — расстояние между частицами.
Таким образом, в ходе моделирования для каждой частицы необходимо получить расстояние между ней и остальными, произвести вычисления по формуле (6) и по формулам (1, 2, 4, 5) и получить координаты частиц на новом временном слое.
Для задания начального распределение скоростей частиц используется функция распределения Максвелла по скоростям:
3
w % ( m ^2 ( mv2 |
F( v) = N [ 2kT J exp І,- 2ЇТ J, (7)
где v — скорость; N — количество частиц; m — масса частицы; k — постоянная Больцмана; T — температура.
При достаточно большом количестве частиц в объеме интегрирования (106 — 107) и шагов по времени (104) для применения подобного алгоритма на практике требуется его параллельная реализация, существенно зависящая от архитектуры высокопроизводительной ЭВМ.
В рамках данного исследования использовалась графическая карта на базе архитектуры nVidia Tesla, программирование для которой имеет ряд особенностей [21]. При программировании применялась программно-аппаратная архитектура CUDA — расширение языка C/C++.
Поставим во взаимно-однозначное соответствие частицу и нить (независимый последовательный исполнитель). Все данные о частицах хранятся в общей памяти графической карты. Частицы собираются в блоки — совокупности нитей, запущенных для решения одной задачи, выполняющихся на одном процессоре и логически отделенных от остальных. Далее блоки собираются в сеть — совокупность блоков, запущенных для решения одной задачи. Такая иерархия позволяет любой
57
частице запросить всю необходимую информацию для исполнения последовательного алгоритма, описанного выше. Для ускорения решения использовались возможности разделяемой памяти процессоров.
Рассмотрим подробнее блок. Каждой его нити требуются одинаковые данные из общей памяти, но запросы в память выполняются медленно. Каждая нить в блоке загружает из общей в разделяемую память свою часть порции данных, производится синхронизация (рис. 2). Нить обрабатывает данные, и так как они лежат в разделяемой памяти, процесс проходит намного быстрее (во время обработки происходит множество запросов в память), затем снова следует синхронизация. Алгоритм повторяется, пока все данные не будут обра-Рис. 2. Схема блока Еютаны.
Результаты моделирования
Алгоритм метода молекулярной динамики, несмотря на кажущуюся простоту (используются классические уравнения механики), имеет высокую сложность порядка O(N2). Практические замеры показали, что программа решает задачу за 40 минут, что более чем в 400 раз быстрее классического подхода.
Для оценки погрешности решения программы применялось правило Рунге. Использовалось 105 частиц, рассчитанных на сетках с шагом 2-10-8 и 10-8. Погрешность составила 10-4 м для каждой из координат, что можно считать сравнительно небольшой величиной.
На рисунке 3 показаны результаты моделирования. В связи с тем, что электроны имеют высокую скорость, влияние спутника будет на них минимально (на рисунке 2 они не показаны). Перед спутником образовывается положительно заряженный фронт, за спутником тянется пустой шлейф, который называется ионной тенью.
Рис. 3. Слева — начальное распределение протонов вокруг спутника, справа — по прошествии 5,71-10 — 5 секунд (на 5710-м шаге)
58
->-<Все обработано?
Нет
Загрузка порции Синхронизация Обработка порции Синхронизация
К следующему этапу
Распределения, полученные в ходе моделирования, удовлетворительно согласуются как с физическим представлением движения спутника в плазме, так и с моделями, построенными другими методами.
Заключение
Описанный метод молекулярной динамики позволяет наиболее полно, по сравнению с другими, моделировать процесс движения спутника в плазме. Рассмотренный вычислительный алгоритм и оптимизация программы для расчета на компьютерах с параллельной архитектурой существенно сокращает время расчета.
Метод показал удовлетворительные результаты для модельной задачи, но при переходе к реальным 3О-вычислениям или при увеличении концентрации его применимость наталкивается на технические трудности. Для модельной задачи при количестве частиц порядка 107 необходимо хранить в оперативной памяти не менее 700 Мб данных для каждого узла, поэтому сегодня дальнейшее масштабирование с использованием описанного метода возможно на суперкомпьютерах с общей памятью. Использование кластерной конфигурации ограничивается скоростью обмена между вычислительными узлами.
Работа выполнена при финансовой поддержке РФФИ по проектам № 12-01-00477а, № 11-01-00098а и № 11-01-00558а.
Список литературы
1. Гальперин Ю. И., Григорьев С. А., Зинин Л. В. и др. Расчеты распределения электрического потенциала вокруг спутника «Авроральный зонд» проекта Интербол // Препринт Институт космических исследований РАН. Пр-1879. 1993
2. Зинин Л. В., Гальперин Ю. И., Григорьев С. А. и др. Об измерениях эффектов поляризационного джета во внешней плазмосфере // Космические исследования. 1998. T. 36. C. 42-52.
3. Ридлер В., Торкар К., Веселов М. В. и др. Эксперимент РОН по активному регулированию электростатического потенциала космического аппарата // Космические исследования. 1998. Т. 36. С. 53—62.
4. Torkar K., Veselov M. V., Afonin V. V. et al. An experiment to study and control the Langmuir sheath around Interball-2 // Ann. Geophys. 1998. Vol. 16. P. 1086 — 1096.
5. Bouhram M., Dubouloz N., Hamelin M. et al. Plasma sheath analysis and current calculations from the potential distribution around Interball Auroral probe // 2001: A spacecraft charging odyssey. Proceeding of the 7th Spacecraft Charging Technology Conference, 2001, ESTEC, Noordwijk, The Netherlands. ESA SP-476. P. 557 — 561.
6. Hamelin M., Bouhram M., Dubouloz N. et al. Combined effects of satellite and ion detector geometries and potentials on the measurements of thermal ions. The Hyperboloid instrument on Interball // 2001: A spacecraft charging odyssey. Proceeding of the 7th Spacecraft Charging Technology Conference, 2001, ESTEC, Noordwijk. The Netherlands. ESA SP-476. P. 569 — 574.
7. Елизарова Т. Г., Четверушкин Б. Н. Об одном вычислительном алгоритме для расчета газодинамических течений // Доклады АН СССР. 1989. Т. 279.
8. Граур И. А., Елизарова Т. Г., Четверушкин Б. Н. Использование кинетических алгоритмов для расчета газодинамических задач, моделирующих вязкие течения. Препринт ИПМ им. Келдыша. 1985. № 85.
9. Ma T.-Z., Schunk R. W. A fluid model of high voltage spheres in the ionosphere // Planet. Space Sci. 1989. Vol. 37. P. 21 — 47.
10. Рылина И. В., Зинин Л. В., Григорьев С. А. и др. Гидродинамический подход к моделированию распределения тепловой плазмы вокруг движущегося заряженного спутника // Космические исследования. 2002. Т. 40, № 4. С. 395 — 405.
11. Zinin L., Grigoriev S., Rylina I. The models of electric field distributions near a satellite // Proceedings of the conference in memory of Yuri Galperin / eds: Ze-lenyi L. M., Geller M. A., Allen J. H. CAWSES Handbook-001. P. 76—83. 2004.
12. Котельников В. А., Котельников М. В., Гидаспов В. Ю. Математическое моделирование обтекания тел потоками столкновительной и бесстолкновитель-ной плазмы. М., 2010.
13. Альперт Я. Л., Гуревич А. В., Питаевский Л. П. Искусственные спутники в разреженной плазме. М., 1964.
14. Allen M. P., Tildsley D. J. Computer simulation of liquids. Oxford: Claredon Press, 1989.
15. Балабаев Н. К., Гривцов А. Г., Шноль Э. Э. Численное моделирование движения линейной полимерной цепочки // Докл. АН СССР. 1975. Т. 220, вып. 5. С. 1096—1098.
16. Лагарьков А. Н., Сергеев В. М. Метод молекулярной динамики в статистической физике // УФН. 1978. Т. 125, вып. 3. С. 409 —448.
17. Stillinger F. H., Weber Th. A. Computer simulation of local order in condensed phases of silicon // Phys. Rev. Bd. 1985. Vol. 31, N. 8. P. 5262—5271.
18. Ihara S., Itoh S., Kitakami J. Mechanisms of cluster implantation silicon: A molecular dynamic study // Phys. Rev. B. 1998. Vol. 58, N. 16. P. 10736 — 10744.
19. Cheng H. P. Cluster-surface collisions: Characters of Xe^ and C20-Si surface bombardment // J. Chem. Phys. 1999. Vol. 111, N. 16. P. 7583 — 7592.
20. Qi L., Young W. L., Sinnott S. B. Effect of surface reactivity on the nuclear of hydrocarbon thin film through molecular-cluster beam deposition // Surf. Sci. 1999. Vol. 426. P. 83.
21. Боресков А. В., Харламов А. А. Основы работы с технологией CUDA. М., 2010.
Об авторах
Зинин Леонид Викторович — канд. физ.-мат. наук, доц., Балтийский федеральный университет им. И. Канта
Email: [email protected]
Ишанов Сергей Александрович — д-р физ.-мат. наук, проф., Балтийский федеральный университет им. И. Канта.
Email: [email protected]
Шарамет Александр Александрович — асп., Балтийский федеральный университет им. И. Канта.
Email: [email protected]
Мациевский Сергей Валентинович — канд. физ.-мат. наук, доц., Балтийский федеральный университет им. И. Канта.
Email: [email protected]
Authors
Dr Leonid Zinin — assistant professor, I. Kant Baltic Federal University.
Email: [email protected].
Dr Sergey Ishanov — professor, I. Kant Baltic Federal University.
Email: [email protected]
Dr Alexandr Sharamet — PhD student, I. Kant Baltic Federal University.
Email: [email protected]
Dr Sergey Matsievsky — assistant professor, I. Kant Baltic Federal University.
Email: [email protected]
