Астрономический журнал, 2023, T. 100, № 4, стр. 305-318

1. ВВЕДЕНИЕ

Горячими экзопланетами называются планеты, орбиты которых расположены близко к родительской звезде (большая полуось орбиты не превосходит 0.1 а.е.) [1]. Одним из ключевых факторов, определяющих состояние атмосферы таких планет, является нагрев излучением звезды, который происходит за счет поглощения излучения в диапазонах мягкого рентгена (1–10 нм) и жесткого ультрафиолета (10–100 нм). Кроме того, важный вклад в нагрев верхней атмосферы могут давать частицы плазмы звездного ветра и магнитосферы. В частности, необходимо учитывать процесс высыпания электронов магнитосферного происхождения [2]. В формировании экзосфер горячих экзопланет определяющую роль играют процессы взаимодействия со звездой. Высокая температура верхней атмосферы сопровождается ее ионизацией и расширением и может приводить к значительному притоку плазмы в магнитосферу. В результате возникают интенсивные газодинамические потоки. Вещество плазмосферы в этом случае уже не является бесстолкновительным, а форма плазмосферы может значительно отличаться от сферической.

Примерами горячих экзопланет являются горячие юпитеры и теплые нептуны. Первый горячий юпитер был открыт еще в 1995 г. [3]. Вследствие близкого расположения к родительской звезде и относительно больших размеров газовые оболочки горячих экзопланет-гигантов могут переполнять свои полости Роша. Это приводит к формированию интенсивных истечений газа не только на ночной, но и на дневной стороне (окрестность внутренней точки Лагранжа ${{L}_{1}}$) планеты [4, 5]. Такую расширяющуюся верхнюю атмосферу можно назвать протяженной оболочкой горячей экзопланеты-гиганта. Она расположена за пределами полости Роша, имеет достаточно большие размеры, относительно высокую плотность и характеризуется значительными отклонениями от сферической формы. На наличие таких протяженных оболочек косвенно указывает избыточное поглощение излучения в ближнем ультрафиолетовом диапазоне, наблюдаемое у некоторых горячих юпитеров и нептунов во время их прохождения на фоне диска родительской звезды [6–12]. Эти выводы подтверждаются непосредственными численными расчетами в рамках аэрономических моделей [13–17].

Трехмерная газодинамическая структура протяженных оболочек горячих юпитеров исследовалась в работах [18–23], где было показано, что в зависимости от параметров в процессе взаимодействия звездного ветра с расширяющейся оболочкой планеты могут формироваться структуры трех основных типов: замкнутые, квазизамкнутые и открытые. При этом тип протяженной оболочки определяет значение темпа потери массы горячего юпитера [18]. В двумерных моделях [24, 25], а также в трехмерной модели [26] эффектами звездного ветра пренебрегалось. В последующих работах этих авторов [27, 28] влияние звездного ветра было учтено. Однако при этом в аэрономической модели задавалось достаточно высокое значение концентрации на фотометрическом радиусе, что приводило, как правило, к формированию оболочки открытого типа с примерно одинаковым значением темпа потери массы.

Оценки величины магнитного поля горячих юпитеров (см., напр., [29]) показывают, что оно является относительно слабым, но вполне достаточным, чтобы оказывать влияние на динамику оболочек. Неоднократно предпринимались попытки учета магнитного поля в одномерных [30–33], двумерных [34] и трехмерных [32, 35–37] численных аэрономических моделях атмосфер горячих юпитеров. В работах с участием автора также исследовалось влияние собственного магнитного поля планеты [38, 38 ] и поля ветра [40–45] на динамику оболочек горячих юпитеров (см. обзор [1] и монографию [2]).

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

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

2. ТЕРНАРНАЯ СФЕРИЧЕСКАЯ СЕТКА

Координаты $x$, $y$, $z$ глобальной декартовой системы координат связаны с координатами $r$, $\theta $, $\varphi $ соотношениями:

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

Наиболее строгим с математической точки зрения подходом такого типа является использование кубизированной сферы [47]. Сетка строится с помощью проекции шести граней куба на поверхность описанной вокруг него сферы. В результате можно получить точное разбиение сферы на шесть непересекающихся одинаковых областей. Локальные угловые координаты в каждой области получаются с помощью проекции декартовых координат на соответствующей грани куба. При этом, однако, получающаяся криволинейная система координат не является сферической и к тому же оказывается неортогональной. Поэтому использование такого подхода не всегда является удобным. Тем не менее его часто используют в астрофизических вычислениях (см., напр., [48, 49]).

Альтернативным подходом является формирование композитной сетки [50] на основе локальных сферических координат. В работе [51] для моделирования трехмерной структуры солнечного ветра использованы три области: одна экваториальная и две полярные. В этом случае достигается минимальное отклонение от глобальной сферической системы координат, поскольку размеры полярных областей можно варьировать. С другой стороны, конфигурации областей отличаются между собой. В композитной сетке Инь-Янь [52] используется разбиение сферы на две одинаковые области, повернутые друг относительно друга на 90°. В этом подходе получается минимальное количество областей с одинаковой конфигурацией. В качестве одного из недостатков иногда указывают на слишком длинную границу между областями, что может замедлять параллельную программу. В работах [53, 54] предложено использовать композитную сферическую сетку, состоящую из шести одинаковых областей, соответствующих шести граням куба. В этом случае сетка строится на основе регулярной конфигурации, а области имеют короткие границы. Однако количество областей относительно велико.

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

Будем нумеровать секторы индексом $l$, пробегающим значения 0, 1, 2. В каждом секторе используются локальные сферические координаты $r$, $\theta $, $\varphi $. При этом угловые переменные изменяются в пределах:

Внутри каждого сектора удобно ввести локальную декартову систему координат ${{x}^{{(l)}}}$, ${{y}^{{(l)}}}$, ${{z}^{{(l)}}}$: Тогда связь между локальными и глобальными декартовыми координатами определяется следующими простыми формулами: Иными словами, сектору 0 соответствует тождественная перестановка глобальных координат $[xyz]$, сектору 1 соответствует перестановка $[yzx]$, а сектору 2 – перестановка $[zxy]$. Преобразования для декартовых компонентов векторов получаются аналогичным образом. Заметим, что в секторе 0 локальные координаты совпадают с глобальными.

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

а для второй точки получаем Назовем эти две точки точками перехлеста.

Наличие этих точек приводит к тому, что граница данного сектора по угловым переменным имеет структуру, изображенную на рис. 2. Граница по координате $\theta $ (соответствует минимальному или максимальному значению $\theta $ при произвольном значении $\varphi $) везде примыкает к одному и тому же сектору (вверху и внизу на рисунке). Граница по координате $\varphi $ (соответствует минимальному или максимальному значению $\varphi $ при произвольном значении $\theta $) примыкает сразу к двум секторам (слева и справа на рисунке). Эти участки границы разделяются точками перехлеста. Нумерация секторов на рисунке осуществляется с помощью циклической перестановки индексов 0, 1, 2. Например, если рассматриваются границы для сектора 1, то сверху окажется сектор 0, а снизу – сектор 2. Точки перехлеста необходимо учитывать в процедуре сшивки секторов для получения полного решения.

Для ускорения трехмерных расчетов вычислительную программу адаптируют для выполнения на нескольких параллельных процессах. В данном подходе параллельная программа организуется следующим образом. Вначале разбиваем сетку в радиальном направлении на несколько слоев. Обозначим через ${{N}_{r}}$ исходное число ячеек в радиальном направлении, а через $M$ – соответствующее число слоев. Тогда ${{N}_{r}} = M{{n}_{r}}$, где ${{n}_{r}}$ – число ячеек в одном радиальном слое. Один процесс обрабатывает один радиальный слой в одном секторе. Иначе говоря, непараллельная программа решает задачу только в одном секторе. Если обозначить через ${{n}_{\theta }}$ и ${{n}_{\varphi }}$ количество ячеек в координатных направлениях $\theta $ и $\varphi $, то полное число ячеек окажется равным

Полное число процессов при этом равно $P = 3{\kern 1pt} M$.

Таким образом, данный процесс ранга $p$ ($0 \leqslant p \leqslant P - 1$) будет решать задачу в некотором радиальном слое с номером $m$ ($0 \leqslant m \leqslant M - 1$) внутри текущего сектора $l$. Нетрудно видеть, что эти числа связаны простым соотношением

Отсюда, зная ранг процесса $p$, можно найти индекс сектора $l$ и номер радиального слоя $m$: где функция ${\text{int}}(x)$ возвращает целую часть вещественного числа $x$.

Поскольку в композитных сетках области могут перекрываться, то в таких случаях для сшивки значений на границах областей необходимо использовать интерполяцию величин [50]. В данном подходе задача сводится к двумерной интерполяции по угловым переменным $\theta $ и $\varphi $ в ячейках, попадающих в области перекрытия двух секторов. При этом скалярные величины (плотность, давление и т.п.) просто интерполируются, а векторные величины (скорость, магнитное поле) нужно еще преобразовать в новую локальную систему координат. В численном коде, описанном ниже, реализовано несколько видов интерполяционных процедур (ближайший сосед, билинейная, биквадратичная и бикубическая). Однако, как показали численные эксперименты, результаты расчетов с использованием билинейной, биквадратичной и бикубической интерполяций различаются слабо. В тестовых расчетах (см. Приложение 8 ) использовалась бикубическая интерполяция, а в расчетах структуры протяженной оболочки горячего юпитера (см. ниже) использовалась билинейная интерполяция.

3. ОПИСАНИЕ МОДЕЛИ

Физическая и математическая постановка задачи полностью соответствует численной модели, описанной в работе [46]. Не будем здесь приводить подробное описание этой численной модели, ограничившись лишь упоминанием основных деталей. Для описания структуры течения в окрестности горячей экзопланеты используется приближение многокомпонентной магнитной гидродинамики. Это, в частности, позволяет учесть сложный химический состав верхней атмосферы и истекающей оболочки. Уравнения записываются в неинерциальной системе отсчета, вращающейся вместе с двойной системой, состоящей из планеты и звезды, вокруг их общего центра масс. Удобно ввести глобальную декартову систему координат ($x$, $y$, $z$), начало которой совпадает с центром планеты. Ось $x$ направим вдоль линии, соединяющей центры звезды и планеты, а ось $y$ направим вдоль орбитального движения планеты. Тогда с учетом выбранных ориентаций осей $x$ и $y$ третья ось $z$ будет совпадать по направлению с вектором орбитальной угловой скорости.

Полное магнитное поле B удобно представлять в виде суперпозиции фонового поля H и возмущения b, которое определяет магнитное поле, индуцируемое токами в самой плазме [1, 38, 40, 55, 56]. Такой подход обеспечивает бóльшую устойчивость численного алгоритма и повышает качество вычислений в разреженных областях с сильным магнитным полем, например, в магнитосфере планеты. В рамках исследуемой задачи фоновое поле H создается токами, распределенными за пределами расчетной области (корона звезды и недра планеты). Поэтому фоновое поле удовлетворяет условию потенциальности, $\nabla \times {\mathbf{H}} = 0$. Это позволяет его частично исключить из уравнений многокомпонентной магнитной гидродинамики [57, 58]. Кроме того, при расчетах структуры течения в окрестности горячего юпитера в нашей модели предполагается, что фоновое магнитное поле является стационарным, $\partial {\mathbf{H}}{\text{/}}\partial t = 0$. Это, в частности, обусловлено тем, что собственное вращение горячего юпитера, для которого мы проводим расчеты, из-за сильных приливных взаимодействий со стороны близко расположенной родительской звезды оказывается синхронизованным с орбитальным вращением. Следовательно, во вращающейся системе отсчета, связанной с орбитальным движением планеты, ориентация ее собственного магнитного поля не будет изменяться со временем.

Вид уравнений многокомпонентной магнитной гидродинамики в сферических координатах приведен в Приложении A . Разностная схема строится на основе метода конечного объема, некоторые детали которого приведены в Приложении B . Алгоритм численного решения уравнений многокомпонентной магнитной гидродинамики состоит из нескольких последовательных этапов, возникающих в результате применения метода расщепления по физическим процессам. На первом этапе численно решаются уравнения многокомпонентной магнитной гидродинамики без учета фонового поля. При этом используется разностная схема Роу-Эйнфельдта-Ошера повышенного порядка аппроксимации, описанная в работе [46]. На втором этапе производится учет фонового магнитного поля [40]. На третьем этапе осуществляется очистка дивергенции магнитного поля с помощью метода обобщенного множителя Лагранжа [59]. На четвертом этапе учитываются внешние силы и эффекты нагрева-охлаждения.

Учет звездного ветра осуществлялся на основе осесимметричной модели [60], опирающейся на хорошо изученные свойства солнечного ветра [61]. В рамках этой модели ветра можно рассчитать профили плотности, скорости и других величин. Это позволяет получать параметры ветра в любой точке ${\mathbf{r}} = (x,y,z)$ расчетной области. При переходе в неинерциальную систему отсчета магнитное поле ветра не меняется, а скорость ветра пересчитывается с учетом орбитальной скорости планеты. Полученные распределения используются для задания начальных условий в области, занимаемой звездным ветром, а также для реализации граничных условий. В модели ветра предполагается политропное уравнение состояния. Это означает, что с точки зрения адиабатической модели в звездном ветре присутствуют эффективные процессы нагрева. Для согласования этих двух моделей в уравнение энергии добавлена соответствующая функция нагрева, которая обеспечивает механизм ускорения ветра.

В начальный момент времени вокруг планеты задавалась сферически симметричная изотермическая атмосфера, распределение плотности в которой определялось из условия гидростатического равновесия. Начальная толщина атмосферы определялась из условия равновесия по полному давлению с веществом звездного ветра в точке лобового столкновения. При моделировании структуры протяженных оболочек горячих юпитеров предполагается водородно-гелиевый состав верхней атмосферы. В начальный момент времени задавался однородный химический состав, когда параметр $\chi = [{\text{He/H}}]$, равный отношению числа ядер гелия к числу ядер водорода, остается постоянным в любом выделенном объеме атмосферы. Учитывались следующие значимые компоненты атмосферы [1, 2]: H, H⁺, H₂, He, He⁺, а также протоны звездного ветра. Концентрация электронов определяется из условия квазинейтральности плазмы, а уравнение движения электронной жидкости в приближении МГД играет роль закона Ома и учитывается в уравнении индукции при определении самосогласованного электрического поля.

Граничные условия на фотометрическом радиусе задавались на основе результатов расчетов, проведенных в рамках одномерных аэрономических моделей с учетом надтепловых частиц [1, 2]. Из аэрономических расчетов следует, что в верхней атмосфере горячего юпитера под воздействием жесткого излучения звезды формируется планетный ветер, который определяет темп потери массы ${{\dot {M}}_{{\text{a}}}} \approx {{10}^{9}}$ г/с. Отсюда можно найти скорость истечения атмосферы на фотометрическом радиусе. Отметим, что эта скорость оказывается достаточно малой (порядка 50 см/с) и поэтому в рамках рассматриваемой модели не приводит к заметному отклонению от условия гидростатического равновесия. Внутри атмосферы в начальный момент времени задавалась радиальная скорость (планетный ветер) с помощью соотношения ${{r}^{2}}\rho {{{v}}_{r}} = {\text{const}}$, выражающего закон сохранения массы, где константа определяется граничными условиями на фотометрическом радиусе.

Для повышения пространственного разрешения в области атмосферы планеты использовалась сетка, экспоненциально сгущающаяся по радиальной координате к центру планеты. Для угловых переменных $\theta $ и $\varphi $ использовалась однородная сетка. Характерная толщина ячейки на фотометрическом радиусе ${{R}_{{{\text{pl}}}}}$ планеты составляла величину $0.005{\kern 1pt} {{R}_{{{\text{pl}}}}}$. Этого оказывалось вполне достаточно для численной реализации условия гидростатического равновесия. Поэтому в приводимых ниже расчетах не использовались никакие процедуры релаксации в атмосфере, чтобы загасить возникающие в ней возмущения. Такие процедуры (см., напр., [18]) применялись нами во всех расчетах, проводимых ранее с использованием декартовых координат. В рамках новой численной модели, основанной на использовании сферических координат, структура атмосферы моделировалась самосогласованно.

4. РЕЗУЛЬТАТЫ РАСЧЕТОВ

В качестве примера приведем результаты расчета структуры течения в окрестности горячего юпитера HD 209458b. Планета имеет массу ${{M}_{{{\text{p}}l}}} = 0.71{\kern 1pt} {{M}_{{\text{J}}}}$ и фотометрический радиус ${{R}_{{{\text{pl}}}}} = 1.38{\kern 1pt} {{R}_{{\text{J}}}}$, где ${{M}_{{\text{J}}}}$ и ${{R}_{{\text{J}}}}$ – масса и радиус Юпитера. Родительская звезда относится к спектральному классу G0 и характеризуется следующими параметрами: масса ${{M}_{{{\text{st}}}}} = 1.1{\kern 1pt} {{M}_{ \odot }}$, радиус ${{R}_{{{\text{st}}}}} = 1.2{\kern 1pt} {{R}_{ \odot }}$, период собственного вращения ${{P}_{{{\text{rot}}}}}{{ = 14.4}^{d}}$. Большая полуось орбиты планеты $A = 10.2{\kern 1pt} {{R}_{ \odot }}$, что соответствует периоду обращения вокруг звезды ${{P}_{{{\text{orb}}}}} = 84.6$ ч.

Температура атмосферы задавалась равной 7500 K, а полная концентрация частиц на фотометрическом радиусе ${{10}^{{11}}}$ см^–3. Химический состав атмосферы определялся значением $\chi = 0.05$, что соответствует массовому содержанию водорода 0.83 и гелия 0.17. Степени диссоциации молекулярного водорода, ионизации атомарного водорода и ионизации атомарного гелия задавались равными 0.25, 0.75 и 0.25 соответственно. Химические реакции, а также реакции ионизации, рекомбинации и диссоциации не учитывались. Следовательно, как и в работе [46], все компоненты плазмы рассматривались в качестве пассивных примесей, переносимых вместе с веществом.

Предполагалось, что величина магнитного момента ${{\mu }_{{{\text{pl}}}}}$ планеты составляет 0.1 от магнитного момента Юпитера. При этом ось магнитного диполя была наклонена на угол 30° к оси вращения, а угол с направлением на центр звезды составлял 130°. Среднее магнитное поле на поверхности звезды задавалось равным 0.01 Гс. Такое значение соответствует сверхальфвеновскому режиму обтекания планеты звездным ветром, при котором формируется протяженная оболочка квазиоткрытого типа (модель “3” из статьи [46]).

Сферическая система координат выбиралась таким образом, чтобы ее центр совпадал с центром планеты. Расчеты проводились для следующих параметров тернарной сферической сетки: ${{N}_{r}} = 256$, ${{n}_{\theta }} = 32$, ${{n}_{\varphi }} = 64$. По радиальной координате шаг сетки задавался неравномерно. Минимальный шаг $\Delta {{r}_{{\min }}} = 0.005{\kern 1pt} {{R}_{{{\text{pl}}}}}$ достигался на внутренней границе, совпадающей с фотометрическим радиусом горячего юпитера. В параллельной программе использовалось 32 радиальных слоя. Расчет проводился с использованием 96 процессоров. Отметим, что использование сферической сетки в данном расчете позволило достичь пространственного разрешения вблизи планеты в 10 раз лучше, чем, например, в расчетах из работы [38], где было задействовано порядка 1000 процессоров.

Результаты моделирования демонстрируют рис. 3 и 4. На этих рисунках представлены распределения плотности (цвет, изолинии), скорости (стрелки) и магнитного поля (сплошные линии) в орбитальной $xy$ (левые панели) и вертикальной $xz$ (правые панели) плоскостях. Плотность выражена в единицах плотности ветра в окрестности планеты, ${{\rho }_{{\text{w}}}} = 2.3 \times {{10}^{{ - 21}}}$ г/см³. Граница полости Роша показана пунктирной линией. Планета расположена в центре расчетной области. Полученное решение соответствует моменту времени порядка трети орбитального периода от начала счета. В течение этого промежутка времени происходит формирование устойчивой квазистационарной структуры. На рис. 3 представлена общая картина течения, а на рис. 4 более подробно показана область вблизи планеты.

В результате обтекания звездным ветром на ночной стороне планеты формируется широкий (порядка нескольких радиусов планеты) водородно-гелиевый турбулентный шлейф. Взаимодействие звездного ветра с оболочкой горячего юпитера приводит к возникновению отошедшей ударной волны. На дневной стороне формируется мощный поток вещества из окрестности внутренней точки Лагранжа ${{L}_{1}}$. Этот поток направлен в сторону звезды и движется против ветра под действием ее гравитации. Вдоль поверхности струи развивается неустойчивость Кельвина-Гельмгольца. Величина потока имеет не слишком большую интенсивность и поэтому струя в значительной степени отклоняется и запирается звездным ветром. В результате происходит формирование протяженной оболочки квазиоткрытого типа. Полученное решение хорошо согласуется с аналогичным решением [46] для случая декартовой сетки. Следует обратить внимание на то, что решение на тернарной сферической сетке (см. рис. 4) является гладким и не содержит численных артефактов ни вдоль полярной оси, ни на границах секторов.

На рис. 5 показаны радиальные профили плотности (левая панель) и радиального компонента скорости ${{{v}}_{r}}$ (правая панель). Красный цвет соответствует дневной стороне (направление на звезду), а синий цвет – ночной стороне (направление от звезды). Пунктирными линиями показаны начальные профили, а жирными линиями – профили, полученные на момент времени, равный половине орбитального периода от начала счета. Следует обратить внимание, что абсцисса определяет высоту над поверхностью планеты ($r = {{R}_{{{\text{pl}}}}}$), выраженную в фотометрических радиусах. Нулевое значение абсциссы соответствует расстоянию в один фотометрический радиус от центра планеты, а высоте, равной ${{R}_{{{\text{pl}}}}}$, соответствует расстояние от центра в два фотометрических радиуса. Взаимодействие со звездным ветром приводит к возникновению возмущений во внутренних частях атмосферы. Однако они незначительны и не разрушают гидростатическое равновесие. Скорость вещества внутри атмосферы всюду остается близкой к нулю за исключением самых поверхностных слоев. Отсюда можно сделать вывод, что пространственного разрешения сетки оказывается вполне достаточно для самосогласованного расчета структуры атмосферы горячего юпитера, чего невозможно было реализовать ранее в расчетах на декартовых сетках. Это позволяет в дальнейшем использовать новую модель, основанную на тернарной сферической сетке, в том числе и для трехмерных аэрономических расчетов.

5. ЗАКЛЮЧЕНИЕ

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

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

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