Scientific journal
Fundamental research
ISSN 1812-7339
"Перечень" ВАК
ИФ РИНЦ = 1,674

NUMERICAL SIMULATION OF TURBULENT FLOWS IN A PROTOPLANETARY DISK

Kuksa M.M. 1 Marov M.Y. 2
1 Obninsk Institute of Nuclear Power Engineering
2 Vernadsky Institute of Geochemistry and Analytical Chemistry of Russian Academy of Sciences
Under observation is the closed set of equations of compressible magnetic hydrodynamics (MHD) in a mean motion scale. The set of equations underlies a two-dimensional model of axisymmetric protoplanetary disk in cylindrical coordinates. In addition to a conventional probability-theoretical averaging the weigh Favre averaging is used in model construction. For the purpose of accounting the magnetic field impact on turbulent flows the modified coefficient of turbulent kinematic viscosity is applied. The parallel code for numerical calculations of the averaged MHD-equation set is developed. In a series of computing experiments the “free” boundary conditions have been selected, so matter could inflow and outflow from the modeling region. The derived steady-state accretion rate is in agreement with observations of classic T Tauri stars. As a result of numerical simulation the distributions of mean density, hydrodynamic velocity and magnetic induction in the turbulized protoplanetary disk are obtained for a distance of 1 au from the star.
protoplanetary disk
turbulent magnetohydrodynamics
numerical simulation
1. Kolesnichenko A.V., Marov M.Ya. Astronomicheskij vestnik. 2009. Vol. 43, no. 5. pp. 424–448.
2. Kolesnichenko A.V., Marov M.Ya. Astronomicheskij vestnik. 2007. Vol. 41, no. 1. pp. 3–23.
3. Kolesnichenko A.V., Marov M.Ya. Astronomicheskij vestnik. 2008. Vol. 42, no. 1. pp. 1–50.
4. Balbus S.A., Hawley J. F. The Astrophysical Journal. 1991. Vol. 376. pp. 214–222.
5. Brandenburg A., Nordlund A., Stein R. F., Torkelsson U. The Astrophysical Journal. 1995. Vol. 446. pp. 741–754.
6. Campbell C.G. Magnetohydrodynamics in binary stars. Dordrecht: Kluwer Academic Publishers, 1997. 306 p.
7. Kuksa M.M. Odessa Astronomical Publications. 2012. Vol. 25. Issue 2. pp. 104–106.
8. Lin D.N.C., Papaloizou J. Monthly Notices of the Royal Astronomical Society. 1980. Vol. 191. pp. 37–48.
9. Lin D.N.C., Pringle J.E. Monthly Notices of the Royal Astronomical Society. 1987. Vol. 225. pp. 607–613.
10. Lynden-Bell D., Pringle J.E. Monthly Notices of the Royal Astronomical Society. 1974. Vol. 168. pp. 603–637.
11. Shakura N.I., Sunyaev R.A. Astronomy and Astrophysics. 1973. Vol. 24. pp. 337–355.
12. Shakura N.I., Sunyaev R.A., Zilinkevich S.S. Astronomy and Astrophysics. 1978. Vol. 62. pp. 179–187.
13. Velikhov E.P. Sov. Phys. JETP. 1959. Vol. 9. pp. 995–998.

Реконструирование аккреционного протопланетного диска является одной из фундаментальных проблем космогонии Солнечной системы. До сих пор открытым остается вопрос о механизме переноса углового момента в аккреционных дисках, поскольку молекулярная вязкость не может обеспечить темп аккреции, наблюдаемый в дисках вокруг молодых звезд солнечного типа. Н.И. Шакура и Р.А. Сюняев [11], Д. Линден-Белл и Дж. Прингл [10] предположили, что повышенная вязкость может быть вызвана турбулентностью. В качестве источников турбулентности предлагались конвекция [8], нелинейная гидродинамическая неустойчивость [12], гравитационная неустойчивость [9] и возмущения, вызванные внешними воздействиями, однако ни один из них не мог обеспечить перенос углового момента за требуемое время. Прогресс был достигнут с привлечением магниторотационной неустойчивости, открытой Е.П. Велиховым [13] и развитой Ст. Бальбусом и Дж. Хоули [4]. Численные расчеты [5] показали, что магниторотационная неустойчивость порождает турбулентность, которая генерирует и поддерживает магнитное поле в присутствии диссипации.

Существование даже слабого магнитного поля существенно усложняет гидродинамические течения в протопланетном диске. В работах А.В. Колесниченко и М.Я. Марова [1, 3] в приближении одножидкостной магнитной гидродинамики получена замкнутая система магнитогидродинамических уравнений масштаба среднего движения, предназначенная для моделирования сдвиговых и конвективных турбулентных течений слабо ионизованной дисковой среды в присутствии магнитного поля. Данная модель послужила отправной точкой для настоящего исследования.

Цель работы ‒ выполнить численный расчет двумерной модели турбулизованного проводящего протопланетного диска, найти распределения физических параметров в диске и сделать вывод о влиянии магнитного поля на турбулентные течения в нем.

Математическая модель

Одним из подходов к моделированию турбулентности является осреднение физических величин по времени (пространству) или по ансамблю возможных реализаций. При построении модели плазмы в состоянии развитой турбулентности используются два оператора осреднения.

Во-первых, теоретико-вероятностное осреднение по Рейнольдсу

Eqn69.wmf

где мгновенное значение величины A представляется в виде суммы осредненной Eqn70.wmf и пульсационной A′ составляющих

Eqn71.wmf

Во-вторых, средневзвешенное осреднение по Фавру

Eqn72.wmf

где мгновенное значение величины A есть сумма средневзвешенного значения ⟨A⟩ и турбулентной флуктуации A′′

Eqn73.wmf

Осредненная система уравнений сжимаемой магнитной гидродинамики для режима развитой изотропной турбулентности в присутствии гравитационного и магнитного полей имеет вид

Eqn74.wmf (1)

где t – время; Eqn75.wmf и Eqn76.wmf – осредненная плотность и ее начальное значение соответственно; Eqn77.wmf – средневзвешенная скорость; Eqn78.wmf – начальная скорость звука; γ = 5/3 – показатель адиабаты; Φgrav – гравитационный потенциал, создаваемый протосолнцем (самогравитация диска не учитывается); Eqn79.wmf – осредненный симметричный тензор деформаций; νturb – коэффициент турбулентной кинематической вязкости; ηturb – коэффициент турбулентной магнитной диффузии; Eqn80.wmf – осредненное собственное поле диска; Bext – внешнее постоянное поле (например, создаваемое протосолнцем или остаточное поле межзвездной среды) Eqn81.wmf – результирующее осредненное магнитное поле.

Применительно к модели протопланетного диска дифференциальные операторы записываются в цилиндрической системе координат (r, φ, z) с учетом ∂/∂φ = 0, поскольку диск считается симметричным относительно оси вращения Oz.

Для моделирования коэффициентов турбулентного переноса используются два подхода. Первый подход основан на (стандартной) α-модели турбулентной вязкости [11]

Eqn82.wmf (2)

где Eqn83.wmf – средневзвешенная угловая скорость жидкой частицы в дифференциально вращающемся диске.

Второй подход основан на модифицированной модели [2], учитывающей влияние магнитного поля на турбулентное течение через посредство пути смешения

Eqn84.wmf (3)

где путь смешения

Eqn85.wmf

где в свою очередь магнитогидродинамическое число Ричардсона

Eqn86.wmf

Коэффициент турбулентной магнитной диффузии считается постоянным и, согласно [6], определяется по аналогии со стандартной моделью

Eqn87.wmf (4)

где ΩK,1 au – кеплеровская угловая скорость на расстоянии 1 а.е. от звезды.

Эмпирический параметр α принимается равным 0,01, что является общепринятым значением в моделях газовых протопланетных дисков.

Численный метод, начальные и граничные условия

Автором разработаны программные модули для среды Matlab, в которых реализована апробированная ранее явная численная схема третьего порядка аппроксимации [7] и средства визуализации результатов моделирования.

Решение системы дифференциальных уравнений (1) ищется в узлах регулярной ортогональной сетки с числом узлов N = Nr×Nz, где Nr – число узлов по оси Or, Nz – число узлов по оси Oz. С каждой стороны сетки дополняется тремя слоями граничных узлов, в которых задаются граничные условия.

Пространственные производные приближаются конечными разностями шестого порядка аппроксимации

Eqn88.wmf(5)

Eqn89.wmf(6)

где Δr – шаг сетки по радиусу. Производные по z вычисляются аналогичным образом.

Производная по времени вычисляется с помощью явной трехэтапной схемы Рунге‒Кутта третьего порядка аппроксимации. Шаг по времени Δt на каждом временном слое определяется из условия Куранта‒Фридрихса‒Леви, что обеспечивает сходимость явной численной схемы.

В начальный момент времени плазма равномерно распределена в области моделирования с плотностью Eqn90.wmf, что при начальной скорости звука Eqn91.wmf и толщине диска h = 0,05 а.е. дает начальную поверхностную плотность Eqn92.wmf. Начальная радиальная и вертикальная скорости равны нулю Eqn93.wmf. Начальная азимутальная скорость кеплеровская Eqn94.wmf. Собственное магнитное поле диска в начальный момент времени отсутствует Eqn95.wmf. Внешнее магнитное поле имеет небольшую вертикальную составляющую Eqn96.wmf.

Граничные условия, позволяющие веществу свободно проникать и покидать область моделирования, сведены в таблицу.

Результаты исследования и их обсуждение

Результаты численного моделирования с модифицированным коэффициентом турбулентной вязкости (3) представлены на рис. 1, 2.

Темп аккреции газа на внешней границе области моделирования rext = 1,1 а.е. согласуется с наблюдательными данными о классических звездах типа Т Тельца

Eqn97.wmf (7)

Результаты численного моделирования с коэффициентом (2) и с модифицированным коэффициентом (3) имеют следующие общие черты:

Параметр

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

в радиальном направлении

в вертикальном направлении

слева

справа

сверху

снизу

Плотность

односторонняя производная

односторонняя производная

фиксированное значение

фиксированное значение

Радиальная скорость

копирование крайнего значения

копирование крайнего значения

асимметричная экстраполяция

асимметричная экстраполяция

Азимутальная скорость

асимметричная экстраполяция

асимметричная экстраполяция

асимметричная экстраполяция

асимметричная экстраполяция

Вертикальная скорость

асимметричная экстраполяция

асимметричная экстраполяция

копирование крайнего значения

копирование крайнего значения

Радиальная индукция

копирование крайнего значения

копирование крайнего значения

копирование крайнего значения

копирование крайнего значения

Азимутальная индукция

копирование крайнего значения

копирование крайнего значения

копирование крайнего значения

копирование крайнего значения

Вертикальная индукция

копирование крайнего значения

копирование крайнего значения

копирование крайнего значения

копирование крайнего значения

pic_23.tif

Рис. 1. Осредненные плотность и гидродинамическая скорость в диске

– Радиальная и азимутальная компоненты магнитного поля меняют знак при переходе через экваториальную плоскость диска.

– Собственное магнитное поле диска развивается (растет по модулю), однако для выяснения причины этого явления требуются дальнейшие исследования.

– В обоих случаях подтвердилась гипотеза, высказанная в [1] о том, что азимутальная компонента магнитной индукции Eqn98.wmf гораздо сильнее меняется по высоте, чем по радиусу, а значит, дает основание учитывать ее градиент в моделях турбулентного переноса.

Результаты, полученные с модифицированным коэффициентом, отличаются тем, что:

– Вертикальная компонента Eqn99.wmf имеет противоположный знак;

– Магнитное поле развивается медленнее, чем в первом случае.

Выводы

В настоящей работе на основе замкнутой системы МГД-уравнений для случая развитой турбулентности разработана математическая модель осесимметричного околосолнечного протопланетного диска конечной толщины. Реализован подход, позволяющий проследить взаимосогласованную эволюцию турбулентных течений и магнитного поля в диске.

pic_24.tif

Рис. 2. Осредненная собственная магнитная индукция в диске

Автором разработаны программные модули для среды Matlab и выполнено численное моделирование при «свободных» граничных условиях, в результате которого были получены распределения осредненной плотности и скорости, конфигурация собственного магнитного поля в диске, а также оценка темпа аккреции, согласующаяся с наблюдательными данными о классических звездах на стадии Т Тельца.

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

Рецензенты:

Камаев Д.А., д.т.н., заведующий лабораторией, ФГБУ Нпо «Тайфун», г. Обнинск;

Перегуда А.И., д.т.н., профессор, Калужский филиал Московского финансово-юридического университета, г. Малоярославец.

Работа поступила в редакцию 01.08.2013.