Краткий учебный проект: численное моделирование формирования протопланетного диска и взаимодействия с зародышем планеты.
Постановка задачи
- Модель: двумерная тонкодисковая гидродинамика (поверхностная плотность $\Sigma (r,\phi ,t)$, поле скорости $\mathbf{v}=(v_r,v_{\phi })$), локально изотермическое уравнение состояния.
- Основные уравнения (в цилиндрических координатах, вертикально усреднённые):
$$\frac{\partial \Sigma }{\partial t}+\nabla \cdot (\Sigma \mathbf{v})=0,$$ $$\frac{\partial (\Sigma \mathbf{v})}{\partial t}+\nabla \cdot (\Sigma \mathbf{v}\otimes \mathbf{v})=-\nabla P-\Sigma \nabla \Phi +\nabla \cdot \mathbf{T},$$ где $P=c_s^2\Sigma$ (локально изотермично), $\Phi$ — потенциал звезды и планеты, $\mathbf{T}$ — вязкий тензор (альфа-модель ниже).
Набор начальных условий
- Системы единиц: принять $G{M}_{\ast }=1$, радиус планеты в начальный момент $r_p=1$, орбитальный период $T_{\mathrm{orb}}=2\pi$.
- Поверхностная плотность: $\Sigma (r)={\Sigma }_0{\left(\frac{r}{1}\right)}^{-p}$, выбрать $p=1$ и нормировку ${\Sigma }_0$ чтобы диск был тонким (учебный выбор ${\Sigma }_0=1$ в код.ед.).
- Температура/напор: аспект отношение диска $h\equiv H/r=0.05$ (локально изотермично $c_s=h{v}_K$, где $v_K=\sqrt{G{M}_{\ast }/r}$).
- Вязкость (альфа-модель): $\nu =\alpha c_{s}H$ с $\alpha =10^{-3}$ (вариантно $\alpha =10^{-4}$ и $\alpha =10^{-2}$ для сравнения).
- Планета: массовая доля $q\equiv M_p/M_{\ast }$ взять набор: $q\in \{10^{-5},10^{-4},10^{-3}\}$ (соответственно суперальс/нептун/юпитер-подобные случаи). Планету фиксировать на круговой орбите $r_p=1$ или разрешить миграцию.
- Пространственная область: радиально $r\in [r_{\mathrm{in}},r_{\mathrm{out}}]=[0.4,2.5]$ (в единицах $r_p$), по азимуту $\phi \in [0,2\pi ]$.
- Граничные условия: заливка волновыми дампинг-зонами у краёв (движение к начальной плотности в зоне $r<0.5$ и $r>2.2$), неабсорбирующие по углу.
Численная настройка (рекомендации для студентов)
- Код: FARGO, PLUTO или Athena++ (FARGO удобен для долгих эволюций благодаря ускорению по азимуту).
- Сетка: тестовый уровень — $N_r\times N_{\phi }=256\times 512$; продвинутый — $512\times 1024$.
- Тайм-степ: CFL-условие; симуляция на промежуток $t_{\mathrm{final}}=200$ орбит (для начального проекта) или до $t=1000$ орбит для изучения медленной эволюции.
- Потенциал планеты смягчить: мягкая длина $ϵ=0.6H$ для 2D аппроксимации.
Ожидаемые результаты (что анализировать)
- Карты поверхности плотности $\Sigma (r,\phi ,t)$: спиральные волны, формирование зазора для больших $q$.
- Азимутальные усреднения $⟨\Sigma {⟩}_{\phi }(r,t)$: профиль зазора (глубина и ширина) во времени.
- Квазистационарный момент ротора/торк на планету $\Gamma (t)$ и скорость миграции ${\dot{r}}_p(t)$ (если планета подвижна).
- Радиальные потоки массы и массовый расход через внутреннюю границу $\dot{M}(r,t)=-2\pi r⟨\Sigma v_r{⟩}_{\phi }$.
- Фурье-распад по азимутальному числу $m$ для оценки доминирующих мод: амплитуда $A_m(r,t)=|\int \Sigma e^{-im\phi }d\phi |$.
Критерии валидации и тесты корректности
- Консервация:
- Относительная ошибка суммарной массы диска за расчёт: ∣M(t)−M(0)M(0)∣<1%\left|\dfrac{M(t)-M(0)}{M(0)}\right|<1\% M(0)M(t)−M(0) <1% (с учётом потерь через граничные зоны).
- Аналогично для углового момента: погрешность должна быть малой ($<5\%$ в учебном расчёте).
- Сходимость по разрешению: ключевые характеристики (глубина зазора, средний торк) не должны сильно меняться при увеличении разрешения (различие $<20\%$ между $256\times 512$ и $512\times 1024$).
- Сравнение с аналитикой:
- Для малых масс (тип I) средний торк должен быть порядка теоретической скалира: в порядковых величинах $\Gamma \sim {\Gamma }_0$, где ${\Gamma }_0\sim {\left(\frac{q}{h}\right)}^2{\Sigma }_p{r}_p^4{\Omega }_p^2$. Сопоставление знака и порядка величины с линейной теорией (Tanaka et al.) — качественный тест.
- Критерий открытия зазора (Crida et al. 2006): проверить условие
$$\frac{3}{4}\frac{H}{R_H}+\frac{50\nu }{q{r}_p^2{\Omega }_p}\lesssim 1,$$ где $R_H=r_p{\left(\frac{q}{3}\right)}^{1/3}$. Для $q=10^{-3}$, $h=0.05$, $\alpha =10^{-3}$ ожидается открытый зазор; для $q=10^{-5}$ — нет.
- Повторяемость: результаты (глубина зазора, типичный торк) стабилизируются при запуске с теми же начальными условиями и различными схемами времени/ограничений в рамках ожидаемой погрешности.
- Диагностика волн: скорость фазовой скорости спиралей близка к локальной угловой скорости возбуждающего источника; проверка сдвига угловой позиции спиралей во времени.
Ожидаемые численные значения (пример)
- Для параметров $q=10^{-3}$, $h=0.05$, $\alpha =10^{-3}$ ожидаемые явления: плотный чистый зазор с глубиной ${\Sigma }_{\mathrm{gap}}/{\Sigma }_0\ll 1$ через $\sim 100$ орбит; торк порядка $|\Gamma |\sim 10^{-5}÷10^{-4}$ в код.ед. (ориентировочно, измерять в ${\Sigma }_0{r}_p^4{\Omega }_p^2$).
- Для $q\le 10^{-4}$ — линейные спиральные волны, миграция типа I, зазора не будет.
Практические задания для студентов
- Базовый: провести симуляцию без миграции планеты для трёх масс $q$ (см. набор) и сравнить профили $⟨\Sigma {⟩}_{\phi }(r,t)$ через $t=50,100,200$ орбит.
- Продвинутый: разрешить миграцию и измерить ${\dot{r}}_p$; сравнить с теоретическим предсказанием типа I (порядки).
- Вариации: изменить $\alpha$ и $h$ и проследить влияние на зазор (проверить Crida).
Критерии зачёта (кратко)
- Наличие корректной реализации уравнений и граничных условий.
- Соблюдение консервации массы/импульса в пределах заданных допусков.
- Анализ: карты $\Sigma$, профиль зазора, график торка $\Gamma (t)$, сравнение с аналитическими критериями (Crida, порядок величины для торка).
- Краткий отчёт (1–2 стр.) с выводами и обсуждением источников ошибок.
Если нужно, дам готовый шаблон начальных файлов для FARGO/PLUTO, контрольные графики и примеры команд запуска.