Реконструирование аккреционного протопланетного диска является одной из фундаментальных проблем космогонии Солнечной системы. До сих пор открытым остается вопрос о механизме переноса углового момента в аккреционных дисках, поскольку молекулярная вязкость не может обеспечить темп аккреции, наблюдаемый в дисках вокруг молодых звезд солнечного типа. Н.И. Шакура и Р.А. Сюняев [11], Д. Линден-Белл и Дж. Прингл [10] предположили, что повышенная вязкость может быть вызвана турбулентностью. В качестве источников турбулентности предлагались конвекция [8], нелинейная гидродинамическая неустойчивость [12], гравитационная неустойчивость [9] и возмущения, вызванные внешними воздействиями, однако ни один из них не мог обеспечить перенос углового момента за требуемое время. Прогресс был достигнут с привлечением магниторотационной неустойчивости, открытой Е.П. Велиховым [13] и развитой Ст. Бальбусом и Дж. Хоули [4]. Численные расчеты [5] показали, что магниторотационная неустойчивость порождает турбулентность, которая генерирует и поддерживает магнитное поле в присутствии диссипации.
Существование даже слабого магнитного поля существенно усложняет гидродинамические течения в протопланетном диске. В работах А.В. Колесниченко и М.Я. Марова [1, 3] в приближении одножидкостной магнитной гидродинамики получена замкнутая система магнитогидродинамических уравнений масштаба среднего движения, предназначенная для моделирования сдвиговых и конвективных турбулентных течений слабо ионизованной дисковой среды в присутствии магнитного поля. Данная модель послужила отправной точкой для настоящего исследования.
Цель работы ‒ выполнить численный расчет двумерной модели турбулизованного проводящего протопланетного диска, найти распределения физических параметров в диске и сделать вывод о влиянии магнитного поля на турбулентные течения в нем.
Математическая модель
Одним из подходов к моделированию турбулентности является осреднение физических величин по времени (пространству) или по ансамблю возможных реализаций. При построении модели плазмы в состоянии развитой турбулентности используются два оператора осреднения.
Во-первых, теоретико-вероятностное осреднение по Рейнольдсу
где мгновенное значение величины A представляется в виде суммы осредненной и пульсационной A′ составляющих
Во-вторых, средневзвешенное осреднение по Фавру
где мгновенное значение величины A есть сумма средневзвешенного значения 〈A〉 и турбулентной флуктуации A′′
Осредненная система уравнений сжимаемой магнитной гидродинамики для режима развитой изотропной турбулентности в присутствии гравитационного и магнитного полей имеет вид
(1)
где t – время; и
– осредненная плотность и ее начальное значение соответственно;
– средневзвешенная скорость;
– начальная скорость звука; γ = 5/3 – показатель адиабаты; Φgrav – гравитационный потенциал, создаваемый протосолнцем (самогравитация диска не учитывается);
– осредненный симметричный тензор деформаций; νturb – коэффициент турбулентной кинематической вязкости; ηturb – коэффициент турбулентной магнитной диффузии;
– осредненное собственное поле диска; Bext – внешнее постоянное поле (например, создаваемое протосолнцем или остаточное поле межзвездной среды)
– результирующее осредненное магнитное поле.
Применительно к модели протопланетного диска дифференциальные операторы записываются в цилиндрической системе координат (r, φ, z) с учетом ∂/∂φ = 0, поскольку диск считается симметричным относительно оси вращения Oz.
Для моделирования коэффициентов турбулентного переноса используются два подхода. Первый подход основан на (стандартной) α-модели турбулентной вязкости [11]
(2)
где – средневзвешенная угловая скорость жидкой частицы в дифференциально вращающемся диске.
Второй подход основан на модифицированной модели [2], учитывающей влияние магнитного поля на турбулентное течение через посредство пути смешения
(3)
где путь смешения
где в свою очередь магнитогидродинамическое число Ричардсона
Коэффициент турбулентной магнитной диффузии считается постоянным и, согласно [6], определяется по аналогии со стандартной моделью
(4)
где ΩK,1 au – кеплеровская угловая скорость на расстоянии 1 а.е. от звезды.
Эмпирический параметр α принимается равным 0,01, что является общепринятым значением в моделях газовых протопланетных дисков.
Численный метод, начальные и граничные условия
Автором разработаны программные модули для среды Matlab, в которых реализована апробированная ранее явная численная схема третьего порядка аппроксимации [7] и средства визуализации результатов моделирования.
Решение системы дифференциальных уравнений (1) ищется в узлах регулярной ортогональной сетки с числом узлов N = Nr×Nz, где Nr – число узлов по оси Or, Nz – число узлов по оси Oz. С каждой стороны сетки дополняется тремя слоями граничных узлов, в которых задаются граничные условия.
Пространственные производные приближаются конечными разностями шестого порядка аппроксимации
(5)
(6)
где Δr – шаг сетки по радиусу. Производные по z вычисляются аналогичным образом.
Производная по времени вычисляется с помощью явной трехэтапной схемы Рунге‒Кутта третьего порядка аппроксимации. Шаг по времени Δt на каждом временном слое определяется из условия Куранта‒Фридрихса‒Леви, что обеспечивает сходимость явной численной схемы.
В начальный момент времени плазма равномерно распределена в области моделирования с плотностью , что при начальной скорости звука
и толщине диска h = 0,05 а.е. дает начальную поверхностную плотность
. Начальная радиальная и вертикальная скорости равны нулю
. Начальная азимутальная скорость кеплеровская
. Собственное магнитное поле диска в начальный момент времени отсутствует
. Внешнее магнитное поле имеет небольшую вертикальную составляющую
.
Граничные условия, позволяющие веществу свободно проникать и покидать область моделирования, сведены в таблицу.
Результаты исследования и их обсуждение
Результаты численного моделирования с модифицированным коэффициентом турбулентной вязкости (3) представлены на рис. 1, 2.
Темп аккреции газа на внешней границе области моделирования rext = 1,1 а.е. согласуется с наблюдательными данными о классических звездах типа Т Тельца
(7)
Результаты численного моделирования с коэффициентом (2) и с модифицированным коэффициентом (3) имеют следующие общие черты:
Параметр |
Граничные условия |
|||
в радиальном направлении |
в вертикальном направлении |
|||
слева |
справа |
сверху |
снизу |
|
Плотность |
односторонняя производная |
односторонняя производная |
фиксированное значение |
фиксированное значение |
Радиальная скорость |
копирование крайнего значения |
копирование крайнего значения |
асимметричная экстраполяция |
асимметричная экстраполяция |
Азимутальная скорость |
асимметричная экстраполяция |
асимметричная экстраполяция |
асимметричная экстраполяция |
асимметричная экстраполяция |
Вертикальная скорость |
асимметричная экстраполяция |
асимметричная экстраполяция |
копирование крайнего значения |
копирование крайнего значения |
Радиальная индукция |
копирование крайнего значения |
копирование крайнего значения |
копирование крайнего значения |
копирование крайнего значения |
Азимутальная индукция |
копирование крайнего значения |
копирование крайнего значения |
копирование крайнего значения |
копирование крайнего значения |
Вертикальная индукция |
копирование крайнего значения |
копирование крайнего значения |
копирование крайнего значения |
копирование крайнего значения |
Рис. 1. Осредненные плотность и гидродинамическая скорость в диске
– Радиальная и азимутальная компоненты магнитного поля меняют знак при переходе через экваториальную плоскость диска.
– Собственное магнитное поле диска развивается (растет по модулю), однако для выяснения причины этого явления требуются дальнейшие исследования.
– В обоих случаях подтвердилась гипотеза, высказанная в [1] о том, что азимутальная компонента магнитной индукции гораздо сильнее меняется по высоте, чем по радиусу, а значит, дает основание учитывать ее градиент в моделях турбулентного переноса.
Результаты, полученные с модифицированным коэффициентом, отличаются тем, что:
– Вертикальная компонента имеет противоположный знак;
– Магнитное поле развивается медленнее, чем в первом случае.
Выводы
В настоящей работе на основе замкнутой системы МГД-уравнений для случая развитой турбулентности разработана математическая модель осесимметричного околосолнечного протопланетного диска конечной толщины. Реализован подход, позволяющий проследить взаимосогласованную эволюцию турбулентных течений и магнитного поля в диске.
Рис. 2. Осредненная собственная магнитная индукция в диске
Автором разработаны программные модули для среды Matlab и выполнено численное моделирование при «свободных» граничных условиях, в результате которого были получены распределения осредненной плотности и скорости, конфигурация собственного магнитного поля в диске, а также оценка темпа аккреции, согласующаяся с наблюдательными данными о классических звездах на стадии Т Тельца.
В ходе исследования подтвердилось предположение о существенном изменении азимутальной компоненты магнитной индукции с высотой, что дает основание учитывать ее градиент в модели коэффициента турбулентной кинематической вязкости через посредство пути смешения.
Рецензенты:
Камаев Д.А., д.т.н., заведующий лабораторией, ФГБУ Нпо «Тайфун», г. Обнинск;
Перегуда А.И., д.т.н., профессор, Калужский филиал Московского финансово-юридического университета, г. Малоярославец.
Работа поступила в редакцию 01.08.2013.
Библиографическая ссылка
Кукса М.М., Маров М.Я. ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ТУРБУЛЕНТНЫХ ТЕЧЕНИЙ В ПРОТОПЛАНЕТНОМ ДИСКЕ // Фундаментальные исследования. – 2013. – № 10-1. – С. 35-39;URL: https://fundamental-research.ru/ru/article/view?id=32210 (дата обращения: 17.08.2024).