Космические исследования, 2023, T. 61, № 6, стр. 447-453

Прототип службы прогноза спокойного солнечного ветра на основе МГД-моделирования и граничных условий модели WSA

С. Арутюнян¹, А. Кодуков¹, М. Субботин¹, Д. Павлов^1, 2, *

¹ Санкт-Петербургский государственный электротехнический университет “ЛЭТИ”
Санкт-Петербург, Россия

² Санкт-Петербургский государственный университет
Санкт-Петербург, Россия

^* E-mail: dapavlov@etu.ru

Поступила в редакцию 01.03.2023
После доработки 12.03.2023
Принята к публикации 13.03.2023

EDN: CAXOKB
DOI: 10.31857/S0023420623600113

Полный текст (PDF)

Аннотация

Создан прототип службы МГД-моделирования спокойного солнечного ветра и прогнозирования скорости и плотности частиц солнечного ветра в межпланетном пространстве, аналогичной службам NOAA и ESA. Служба состоит из МГД-симулятора, модуля обработки результатов симуляции и веб-интерфейса. Симулятор основан на реализации метода TVDLF в пакете PLUTO. Граничные условия модели (плотность, радиальная скорость, магнитное поле, температура) на расстоянии 0.1 а. е. от начала координат получаются регулярно из соответствующей службы NOAA, в которой они рассчитаны по модели WSA на основании магнитограмм сети GONG. Доступны два режима граничных условий: постоянные и суточные. Симуляции проводились на равномерной сетке в диапазоне 0.1–1.7 а. е. по расстоянию (512 элементов), –60°…+60° по широте (60 элементов), 0°–360° по долготе (180 элементов). Проведено сравнение рассчитанных карт скорости и плотности частиц с расчетами NOAA SWPC и NASA CCMC при одинаковых граничных условиях. Проведено ретроспективное сравнение получаемых прогнозов с данными прямых измерений (OMNI).

ВВЕДЕНИЕ

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

Опасность, которую несет в себе солнечный ветер, привела к необходимости разработки служб прогнозирования солнечной погоды, а развитие методов численного моделирования и технологий мониторинга солнечной активности сделало создание таких служб возможным. Так, в США в составе Национального управления океанических и атмосферных исследований работает Центр прогноза космической погоды (англ. National Oceanic and Atmospheric Administration – Space Weather Prediction Center, NOAA – SWPC) (https://www.swpc. noaa.gov/), выпускающий регулярные прогнозы на основе трехмерной МГД-модели (магнитной гидродинамики) солнечного ветра ENLIL [1]. В 2012 г. в Европейском космическом агентстве (англ. European Space Agency, ESA) было создано Управление космической погоды (англ. Space Weather Office, SWO) (https://swe.ssa.esa.int/current-space-weather/), которое к настоящему моменту получает прогнозы аналогичные SWPC с помощью программы EUHFORIA [2]. В России служба прогнозирования солнечного ветра разработана в кисловодской Горной астрономической станции Главной (Пулковской) астрономической обсерватории (ГАС ГАО) (http://solarstation.ru/sun-service/forecast), прогнозы в этой службе основаны на кинетической модели [3]. В настоящей работе описан прототип новой службы прогноза солнечного ветра, аналогичной SWPC и SWO и основанной на трехмерной МГД-модели солнечной плазмы.

МЕТОД МОДЕЛИРОВАНИЯ

МГД-уравнения

Моделирование солнечной плазмы как жидкости, подчиняющейся МГД-уравнениям, осуществляется в секторе полого шара, внутренняя граница которого имеет радиус 0.1 а. е. (22.5 R_S). Уравнения МГД выглядят следующим образом:

(1)

$\frac{{\partial p}}{{\partial t}} + \nabla \cdot (p{\mathbf{v}}) = 0,$

(2)

$\frac{{~\partial {\mathbf{m}}}}{{\partial t}} + \nabla \cdot {{\left[ {{\mathbf{mv}} - {\mathbf{BB}} + I\left( {p + \frac{{{{{\mathbf{B}}}^{2}}}}{2}} \right)} \right]}^{T}} = ~ - p\nabla \Phi ,$

(3)

$\frac{{\partial {\mathbf{B}}}}{{\partial t}} + \nabla \times ( - {\mathbf{v}}\, \times \,{\mathbf{B}}) = 0,$

(4)

$\begin{gathered} \frac{{\partial \left( {{{E}_{{\text{t}}}}} \right)}}{{\partial t}} + \\ + \,\,\nabla \left[ {\left( {\frac{{\rho {{{\mathbf{v}}}^{2}}}}{2} + \rho e + p + \rho \Phi } \right){\mathbf{v}} + ( - {\mathbf{v}} \times {\mathbf{B}}) \times {\mathbf{B}}} \right] = 0, \\ \end{gathered} $

(5)

${{E}_{t}} = \frac{\rho }{{\Gamma - 1}} + \frac{{{{m}^{2}}}}{{2\rho }} + \frac{{{{B}^{2}}}}{2},$

где $\rho $ – плотность массы протонов; v – скорость; ${\mathbf{m}}~ = \rho {\mathbf{v}}~$ – плотность импульса; p – давление; B – вектор магнитной индукции; $\Phi $ – гравитационный потенциал Солнца; $\Gamma $ – показатель адиабаты, $\Gamma = 1.5$ [4]. Температура T не входит в уравнения (1)–(5), но связана с плотностью и давлением уравнением состояния: $p = k\rho {T \mathord{\left/ {\vphantom {T {({{m}_{u}}\mu )}}} \right. \kern-0em} {({{m}_{u}}\mu )}},$ где k – постоянная Больцмана; ${{m}_{u}}~$ – атомная единица массы; $\mu $ – средняя молярная масса частиц (принято значение $\mu = 0.616$). Здесь (1) – уравнение непрерывности, (2) – уравнение движения, (3) – уравнение магнитного поля, (4) – уравнение полной энергии, (5) – формула плотности полной энергии, являющейся суммой тепловой, кинетической и магнитной энергии.

В ENLIL для решения системы уравнений, аналогичной (1)–(4), используется метод уменьшения полной вариации с разностной схемой Лакса–Фридрихса (англ. Total Variation Diminishing Lax – Friedrich, TVDLF) [5]. В EUHFORIA используется метод Харта–Лакса–ван Леера (англ. Harten – Lax – van Leer, HLL), основанный на решении задачи Римана. Коды обеих программ не являются публично доступными. В настоящей работе представлена программа, основанная на пакете PLUTO (http://plutocode.ph.unito.it/) [6], который содержит реализации численных методов решения систем уравнений в частных производных, в том числе разновидностей МГД для широкого спектра астрофизических моделей. В настоящей работе используется наиболее простая система (1)–(4), не содержащая сопротивления, вязкости и диссипации.

Граничные условия

Внутренние граничные условия уравнений (1)–(4) рассчитываются на основе эмпирической модели Wang–Sheeley–Arge (WSA) [7]. Входными данными модели WSA служит синоптическая карта магнитограммы Солнца, а выходными – корональные карты (граничные условия) скорости частиц и напряженности магнитного поля на границе внешней короны (21.5 R_S). Выходные данные WSA преобразуются в граничные условия ($\rho ,$ ${{v}_{r}},$ ${{B}_{r}},$ ${{B}_{\varphi }},$ $T$) с использованием условия постоянства момента количества движения и некоторых предварительно заданных базовых значений $\rho ,$ ${{v}_{r}}$ и $T,$ с помощью программы wsa2bc. Сама программа отсутствует в свободном доступе, но итоговые граничные условия, рассчитанные с помощью WSA 2.2 на основе магнитограмм, предоставляемых сетью солнечных обсерваторий GONG (англ. Global Oscillation Network Group) (https://gong.nso.edu/); эти граничные условия используются в настоящей работе.

Несмотря на то, что действительная скорость вращения Солнца зависит от широты, в используемых граничных условиях, основанных на модели WSA, период вращения Солнца считается равным кэррингтоновкому периоду 25.38 сут в инерциальной системе координат (СК) или 27.2753 сут (в СК HEEQ). В связи с тем, что используемые граничные условия представлены в СК HEEQ, моделирование также происходит в СК, совпадающей с HEEQ в начальный момент времени и вращающейся с постоянной угловой скоростью ${{\Omega }_{z}}_{~} = {{2\pi } \mathord{\left/ {\vphantom {{2\pi } {365.25}}} \right. \kern-0em} {365.25}},$ что является достаточно хорошим приближением HEEQ. При этом в уравнения (1)–(4) вносятся члены, ответственные за центробежное и кориолисово ускорение. Граничные условия в ходе моделирования вращаются с указанным синодическим кэррингтоновским периодом 27.2753 сут.

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

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

(6)

$\begin{gathered} \rho (r) = \frac{{0.01~\bar {\rho }}}{{{{r}^{2}}}},\,\,\,\,B(r) = \frac{{0.01~\bar {B}}}{{{{r}^{2}}}}, \\ T(r) = \frac{{0.01~\bar {T}}}{{{{r}^{2}}}},\,\,\,\,{{v}_{r}}(r) = \bar {v}, \\ \end{gathered} $

где $\bar {\rho },$ $\bar {B},$ $\bar {T},$ $\bar {v}$ – средние значения граничных условий на дату начала прогноза.

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

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

ТЕХНИЧЕСКИЕ ПОДРОБНОСТИ

Используется равномерная трехмерная сетка в сферической СК в диапазоне 0.1–1.7 а. е. по расстоянию (512 элементов), –60°…+60° по широте (60 элементов), 0°–360° по долготе (180 элементов). Шаг по времени постоянный, равен 0.0005 сут (при увеличении шага численный метод расходится).

В PLUTO реализована поддержка параллельных вычислений посредством интерфейса передачи сообщений (англ. Message Passing Interface, MPI). При расчете на шести ядрах процессора AMD EPYC 7742 выталкивание начальных условий занимает 187 мин, а дальнейшая симуляция на 5 сут занимает 91 мин.

Результаты моделирования сохраняются в файлы формата dbl (внутренний формат PLUTO; одному файлу соответствует трехмерная карта значений основных переменных модели). Сохранение происходит каждый час модельного времени; таким образом, прогноз на 5 сут состоит из 121 карты. Также есть возможность сохранения карт в формате VTK.

ОБЗОР СОЗДАННОГО ПРОТОТИПА СЛУЖБЫ

Связь компонент службы

На рис. 1 изображена схема связи модулей разработанной службы. Модуль загрузки граничных условий запрашивает данные с веб-сайта NOAA. Модуль запуска PLUTO ожидает новых входных данных; как только данные будут получены, запускается численное моделирование, основанное на PLUTO. Затем на основе полученных трехмерных карт строятся двумерные изображения плотности и скорости частиц; далее двумерные изображения объединяются в видеофайл и передаются на веб-сайт. Схема устроена таким образом, чтобы можно было в случае ошибок вручную запускать отдельные модули.

Рис. 1.

Схема взаимодействия компонент службы.

Построение двумерных изображений

Модуль построения изображений использует значения плотности и скорости частиц на дискретной сетке размера 512 × 60 × 180, полученные в ходе МГД-моделирования на различные моменты времени. Изображения представляют собой срезы в двух плоскостях: плоскости, параллельной плоскости солнечного экватора и пересекающей Землю и плоскости YZ СК HEEQ. Для получения изображений, разрешающая способность которых превосходит разрешающую способность сетки, применяется трилинейная интерполяция.

Рис. 2.

Разработанный веб-сайт прогноза спокойного солнечного ветра.

Широта Земли в СК HEEQ рассчитывается на основе эфемериды EPM2021 [8], разработанной в Институте прикладной астрономии РАН. Там же разработана онлайн-служба (https://iaaras.ru/dept/ephemeris/online/), которая позволяет рассчитать эфемериду геоцентрических положений Солнца ES^eq в экваториальной СК. Гелиоцентрические координаты Земли в HEEQ рассчитываются по следующей формуле:

$S{{E}^{{{\text{HEEQ}}}}} = {{R}_{z}}(\alpha )~{{R}_{x}}( - 26.13)~{{R}_{z}}( - 16.13)~[ - E{{S}^{{eq}}}],$

где $\alpha $ выбирается таким образом, чтобы выполнялось условие $SE_{y}^{{{\text{HEEQ}}}} = 0.$

Веб-сайт

Построенные 121 растровые изображения с помощью пакета FFmpeg конвертируются в видеофайл формата mp4, при этом задается количество кадров в секунду, который передается на разработанный веб-сайт прогноза солнечного ветра (http://solarwind.entroforce.ru/). На рис. 3 представлена основная часть сайта с прогнозом.

Рис. 3.

Карты скорости и плотности частиц в PLUTO (сверху) и SWPC (снизу) за 26.XII.2022.

Большинство браузеров поддерживает функцию перехода в видео на определенный кадр; информация о количестве кадров в секунду позволяет реализовать на сайте полосу прокрутки для выбора отдельных кадров в видеофайле. Данный подход достаточно эффективен, так как предоставляет возможность почасового просмотра прогноза при объеме данных, передаваемых на клиентскую систему, около 1.5 Мбайт за 5 сут.

АНАЛИЗ РЕЗУЛЬТАТОВ МОДЕЛИРОВАНИЯ

Сравнение с SWPC

Служба SWPC позволяет запрашивать архивные данные о симуляции на конкретную дату, однако полное сравнение трехмерной сетки для SWPC провести невозможно. В предоставляемых архивах содержатся графические изображения и исходные данные этих изображений (срезы в горизонтальной плоскости широты Земли и плоскости YZ СК HEEQ, а также графики скорости и плотности частиц в точках, соответствующих положению Земли и космических обсерваторий STEREO A/B).

Сравнение проводилось в плоскости широты Земли. На рис. 3 представлено качественное сравнение карт нормированной плотности и скорости солнечного ветра. В верхней части изображены карты из PLUTO, а в нижней части находятся карты из SWPC. За исключением цветовой палитры изображения сопоставимы. В табл. 1 приведено численное сравнение значений в рассматриваемой плоскости, можно заметить, что средняя разность значений составляет около 3.7% от среднего для скорости и 14.3% для нормированной плотности.

Таблица 1.

Сравнение результатов моделирования солнечного ветра PLUTO и ENLIL (моделирование выполнено в SWPC)

Величина	МГД-Модель	min	max	mean	Средняя абсолютная разность
$v$, км	PLUTO	282.1	639.0	403.1	14.9
$v$, км	ENLIL (SWPC)	272.9	647.4	400.9	14.9
ρ, см^–3	PLUTO	1.8	13.9	4.8	0.7
ρ, см^–3	ENLIL (SWPC)	1.9	14.2	4.9	0.7

Сравнение с CCMC

Центр моделирования, координируемый сообществом – CCMC (англ. Community Coordinated Modeling Center) (https://ccmc.gsfc.nasa.gov/) – служба, предоставляющая доступ к различным реализациям алгоритмов моделирования космической погоды. Данный ресурс позволяет запрашивать результаты расчетов программ, входные данные для которых определяются пользователем. В том числе CCMC предоставляет доступ к запуску ENLIL, благодаря чему существует возможность сверить результат ее работы с вычислениями, проведенными в настоящей работе. В отличие от SWPC, CCMC предоставляет численные результаты симуляции в виде трехмерной сетки. Граничные условия, используемые в симуляциях CCMC, также доступны в составе расширенного архива данных (raw). Последний доступный цикл вращения Солнца для симуляции в CCMC на момент написания статьи 7.VI.2021–4.VII.2021, поэтому для сравнения был выбран именно он.

Сравнение с CCMC проводилось аналогично сравнению с SWPC. Карты плотности и скорости представлены на рис. 4, в верхней части изображения находятся карты PLUTO, а в нижней – карты из CCMC. Численные результаты представлены в табл. 2, в этом случае средняя разность значений составляет 32.4% от среднего значения плотности и 11.2% от среднего значения скорости.

Рис. 4.

Карты скорости и плотности частиц в PLUTO (сверху) и CCMC (снизу) за 7.VI.2021.

Таблица 2.

Сравнение результатов моделирования солнечного ветра PLUTO и ENLIL (моделирование выполнено в CCMC)

Величина	МГД-Модель	min	max	mean	Средняя абсолютная разность
$v$, км	PLUTO	199.8	650.4	373.1	41.8
$v$, км	ENLIL (CCMC)	200.6	650.4	374.8	41.8
ρ, см^–3	PLUTO	1.9	116.1	9.9	3.2
ρ, см^–3	ENLIL (CCMC)	2.3	95.5	10	3.2

Сравнение с OMNI

OMNI (https://omniweb.gsfc.nasa.gov/) – регулярно обновляемый набор данных, содержащий информацию о параметрах солнечного ветра, собираемый различными космическими аппаратами (КА). В частности измерения КА ACE и Wind позволяют оценить плотность и скорость частиц солнечной плазмы вблизи Земли.

Был проведен расчет прогноза скорости и плотности частиц солнечного ветра на основе моделирования в течение 5 сут с 26.XII.2022 по 31.XII.2022 как с постоянными граничными условиями, так и с суточными. Результаты моделирования представлены на рис. 5.

Рис. 5.

Скорость и плотность частиц вблизи Земли с 26.XII.2022 по 31.XII.2022.

ЗАКЛЮЧЕНИЕ

Разработан прототип службы регулярного прогнозирования солнечного ветра по трехмерной МГД-модели на основе математического пакета PLUTO. Ежедневно обновляемые прогнозы на 5 сут представлены на сайте (http://solarwind.entroforce.ru/). Полученные карты солнечного ветра в целом соответствуют зарубежным аналогам при одинаковых граничных условиях моделирования. Есть различия, которые, по всей видимости, обусловлены следующими причинами:

• модель ENLIL в настоящее время учитывает эффект объемного нагрева плазмы (англ. volumetric heating). В реализации PLUTO данный эффект отсутствует; значение показателя адиабаты Γ = 1.5 позволяет частично компенсировать этот недостаток [2], но не в полной мере;

• значения средней молярной массы μ в реализациях ENLIL, использующихся в CCMC и SWPC, не приводятся; они могут отличаться от значения μ = 0.616, использованного в данной работе;

• при математической эквивалентности МГД-уравнений и численных схем для их решения, используемых в ENLIL и PLUTO, различия в реализациях могут оказывать влияние на получаемые результаты.

Данный вопрос требует дальнейшего исследования.

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

• моделирование корональных выбросов массы;

• учет в МГД-модели дополнительных эффектов, в частности сопротивления;

• использование непосредственных данных модели WSA и восстановление граничных условий по этим данным (сейчас служба использует вторичные данные, рассчитанные программой wsa2bc);

• разработка совместной с Кисловодской ГАС службы прогноза, использующей граничные условия, рассчитанные по отечественным наблюдательным данным [9].

Еще одним интересным направлением представляется использование рассчитанных карт электронной плотности для учета задержки в солнечной плазме сигнала между Землей и КА на орбите Марса или Меркурия, что важно для построения высокоточных эфемерид Солнечной системы [10].

Авторы благодарны Душану Одстрчилу (англ. Dusan Odstrcil, George Mason University) за ценные советы. Авторы также благодарны Алексею Печёркину и Наталье Странниковой (ЛЭТИ) за деятельное участие на раннем этапе работы. Авторы благодарны коллективу CCMC за возможность проведения симуляций (симуляция, проведенная в ходе настоящей работы, доступна на сайте CCMC по идентификатору Aleksandr_Kodukov_101722_SH_1).

Исследование Д. А. Павлова выполнено в Санкт-Петербургском международном математическом институте имени Леонарда Эйлера при финансовой поддержке Министерства науки и высшего образования Российской Федерации (соглашение № 075–15–2022–287 от 06.04.2022).

Список литературы

Odstrcil D. Modeling 3-D solar wind structure // Advances in Space Research. 2003. V. 32. Iss. 4. P. 497–506. https://doi.org/10.1016/S0273-1177(03)00332-6
Pomoell J., Poedts S. EUHFORIA: European heliospheric forecasting information Asse 2.0t // J. Space Weather Space Climate. 2018. V. 8. Art. ID. A35. https://doi.org/10.1051/swsc/2018020
Tlatov A.G., Berezin I.A., Strelkov M.A. Simulation of Coronal Mass Ejection Propagation Based on Data from Ground-Based Patrol Observations // Geomagnetism and Aeronomy. 2019. V. 59. Iss. 7. P. 843–845. https://doi.org/10.1134/S0016793219070247
Odstrcil D., Mays M.L., Hess Ph. et al. Operational Modeling of Heliospheric Space Weather for the Parker Solar Probe // The Astrophysical J. Supplement Ser. 2020. V. 246. Iss. 2. Art. ID. 73. 19 p. https://doi.org/10.3847/1538-4365/ab77cb
Tóth G., Odstrčil D. Comparison of Some Flux Corrected Transport and Total Variation Diminishing Numerical Schemes for Hydrodynamic and Magnetohydrodynamic Problems // J. Computational Physics. 1996. V. 128. Iss. 1. P. 82–100. https://doi.org/10.1006/jcph.1996.0197
Mignone A., Bodo G., Massaglia S. et al. PLUTO: a numerical code for computational astrophysics // The Astrophysical J. Supplement Series. 2007. V. 170. Iss. 1. Art. ID. 228. https://doi.org/10.1086/513316
Arge C.N., Pizzo V.J. Improvement in the prediction of solar wind conditions using near-real time solar magnetic field updates // J. Geophysical Research: Space Physics. 2000. V. 105. Iss. A5. P. 10465–10479. https://doi.org/10.1029/1999JA000262
Pitjeva E., Pavlov D., Aksim D., Kan M. Planetary and lunar ephemeris EPM2021 and its significance for Solar system research // Proc. Intern. Astronomical Union. 2019. V. 15. Iss. S364. P. 220–225. https://doi.org/10.1017/S1743921321001447
Berezin I., Tlatov A. Coronal Field Geometry and Solar Wind Speed // Universe 2022. V. 8. Iss. 12. Art. ID. 646. https://doi.org/10.3390/universe8120646
Aksim D., Pavlov D. Improving the solar wind density model used in processing of spacecraft ranging observations // Monthly Notices of the Royal Astronomical Society. 2022. V. 514. Iss. 3. P. 3191–3201. https://doi.org/10.1093/mnras/stac1229

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска содержание выпуска

Космические исследования

Архивы выпусков Информация о журнале Отправить рукопись в журнал