Астрономический вестник, 2019, T. 53, № 3, стр. 195-213

Основные допущения задачи

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

i) Первичные пылевые частицы (мономеры) – однородные по составу, твердые и недеформируемые, сферичные по форме и монодисперсны.

ii) Предполагается несжимаемость материала мономеров, ${{\rho }_{0}} = {\text{const}}.$

iii) Истинная плотность материала мономеров много больше истинной плотности газовой составляющей дисковой системы, ${{\rho }_{0}} \gg {{\tilde {\rho }}_{g}}.$

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

v) Несущая фаза – сжимаемый совершенный газ.

vi) Вязкостными эффектами для несущей и дисперсной фазы можно пренебречь.

vii) Предполагается условие термического равновесия газовой и дисперсной фаз, ${{T}_{\operatorname{g} }} = {{T}_{{\text{p}}}} = T.$

viii) Предполагается, что фрактальная среда пылевых кластеров внутри элементарного (совокупного) макрообъема $\delta V$ имеет фрактальную размерность D, а размерность на его границе $\delta W$ равна d (в общем случае размерность d не равна ни 2, ни D – 1).

ix) При исследовании гидродинамической устойчивости газопылевого слоя можно пренебречь турбулентной диффузией пыли, вызванной турбулентностью газовой фазы¹¹.

Основные параметры фрактальных структур

Далее все первичные мелкодисперсные компактные частицы газопылевого протопланетного облака вне зависимости от их реальной формы, размера и материала будем считать твердыми сферами, имеющими один и тот же радиус ${{r}_{0}}$ и массу ${{m}_{0}},$ поскольку форма мономера (сферическая, эллипсоидальная и т.п.) оказывает незначительное влияние (≲2%) на фрактальную размерность образующихся кластеров (см., например, Bertini и др., 2009). На первоначальном этапе роста (на основе механизма “частица–частица”) пылевые образования состоят из небольшого числа первичных мономеров и по этой причине не могут, вообще говоря, считаться фракталами. Но, по мере дальнейшего слипания мономеров в кластеры, механизмом “кластер–кластер” формируются фрактальные пылевые агрегаты относительно крупных размеров. При этом число ${{n}_{0}}$ первичных пылевых частиц (ядер), входящих в состав изотропного фрактального пылевого агрегата (ФПА) и масса кластера ${{m}^{{{\text{cl}}}}}$ определяются следующими асимптотическими формулами (см., например, Смирнов, 2011)

в которых $R = {{\left( {\frac{5}{{6N}}\sum\nolimits_{i,k = {\text{1}}}^N {{{{\left| {{{r}_{i}} - {{r}_{k}}} \right|}}^{2}}} } \right)}^{{1/2}}}$ – радиус гирации (характерный размер изотропного кластера), определяемый как среднеквадратичный радиус агрегата, измеряемый от его центра тяжести (Okuzumi и др., 2009); $N$– число мономеров, принадлежащих кластеру; ${{r}_{i}}$ – радиус вектор $i$-го мономера, входящего в кластер; – фрактальная массовая размерность кластера, являющаяся количественной характеристикой заполнения им евклидова пространства размерности $n\;( = 1,2,3),$ а также характеризующая самоподобие его внутренней структуры. При этом массовая размерность $D$ не зависит от того, является ли упаковка сфер радиуса ${{r}_{0}}$ плотной, случайной или пористой с однородным распределением пустот. Реальная структура кластеров характеризуется чрезвычайно сложной и нерегулярной геометрией. Хотя размерность $D$ и не отражает полностью геометрические свойства фрактала, она, однако, позволяет учитывать основные свойства фрактальных структур при моделировании широкого класса явлений. Дальше мы ограничимся рассмотрением случая гомогенной фрактальной среды, когда степенной закон (1) не зависит от расположения (перемещения) кластеров в пространстве. Кроме этого, будем предполагать, что внутренняя структура всех кластеров, формирующихся в аэродисперсной дисковой среде, одинакова, поскольку они образуются приблизительно в одних и тех же условиях изучаемых областей диска. Именно по этой причине численные значения массовой размерности $D$ кластеров (при однотипном способе их сборки) будем считать одинаковыми для всего ансамбля разномасштабных ФПА в диске. Заметим, что дробная массовая размерность $D$ пылевых кластеров в дисковой среде определяет, в конечном счете, их аэродинамические свойства, стабильность и динамику роста (см. Wiltzius, 1987; Chen и др., 1987), а также пространственно-временную эволюцию этих рыхлых образований в диске (Смирнов, 1991).

При объединении большого числа малых пылевых кластеров (в результате процесса кластер-кластерной ассоциации) образуются обладающие самоподобными свойствами на малых расстояниях^{22однородные ворсистые агрегаты, в которых по мере увеличения их радиуса гирации увеличиваются размеры пустот, а средняя плотность (средняя массовая плотность вещества в объеме, занимаемом кластером) убывает по закону ${{\bar {\rho }}^{{{\text{cl}}}}} = {{\rho }_{0}}\left( {{{{{r}_{0}}} \mathord{\left/ {\vphantom {{{{r}_{0}}} R}} \right. \kern-0em} R}} \right){{\,}^{{3 - D}}},$ где ${{\rho }_{0}} = {{3{{m}_{{\text{0}}}}} \mathord{\left/ {\vphantom {{3{{m}_{{\text{0}}}}} {4\pi r_{0}^{3}}}} \right. \kern-0em} {4\pi r_{0}^{3}}}$ – массовая плотность материала первичных ядер. Отсюда следует вывод, что чем больший объем занимают части фрактального агрегата, тем больше пустот всех размеров он содержит. Таким образом, одним из характерных свойств фрактального агрегата является его способность захватывать большое пространство (за счет создания ажурной, сильно разветвленной структуры) при использовании меньшего по сравнению с плотным (компактным) агрегатом количества вещества, необходимого для его образования. Компактность и физические свойства отдельного кластера зависят как от характера движения первичных мономеров и кластеров (прямолинейное или броуновское) до столкновения, так и от вероятности слипания мономеров и кластеров при их соприкосновении. В зависимости от числовой плотности ${{N}_{1}}(r,t)$ мономеров (не входящих в состав кластеров) в единице объема $\delta V \equiv \delta r( = \delta x\delta y\delta z)$ дисковой среды возможны два механизма роста кластеров с фрактальной структурой: в результате прилипания к кластеру мономеров или благодаря процессу кластер-кластерной агрегации. При этом в первом случае процесс роста кластера может происходить либо в результате присоединения к нему единичных ядер, двигающихся прямолинейно (кинетический режим), либо когда много первичных пылевых мономеров, двигающихся диффузионно, одновременно объединяются с кластером (диффузионный или гидродинамический режим).}

Далее эволюционирующее газопылевое протопланетное ламинарное облако будем рассматривать как гетерогенный термодинамический комплекс, состоящий из двух взаимопроникающих континуумов, которые заполняют одновременно один и тот же объем евклидова пространства непрерывно − газовой фазы солнечного состава и полидисперсной фазы пылевых частиц (фрактальная среда с нецелой массовой размерностью $D,$ меньшей размерности $n$ координатного пространства задачи), которые находятся при общей абсолютной температуре $T(r,t)$ и давлении $P(r,t).$ Газовую фазу, являющуюся несущей средой, будем описывать моделью идеальной жидкости. В свою очередь, полидисперсную пылевую фазу будем считать фрактальной средой, состоящей из нескольких фракций: фракции первичных пылевых конденсированных мономеров и внедренного в них большого числа фракций фрактальных агрегатов, отличающихся друг от друга размерами. Другими словами, элементарный макрообъем $\delta V$ дисковой среды, помимо несущей газовой фазы, описываемой обычными структурными параметрами, такими как числовая ${{n}_{{\text{g}}}}(r,t)$ и массовая ${{\rho }_{{\text{g}}}}(r,t) = {{n}_{{\text{g}}}}{{m}_{{\text{g}}}}$ плотности (“размазанные” по совокупному элементарному объему смеси), давление ${{P}_{{\text{g}}}}(r,t),$ температура ${{T}_{{\text{g}}}}(r,t),$ гидродинамическая скорость ${{\text{v}}_{{\text{g}}}}(r,t)$ и др., содержит еще множество пылевых фрактальных образований (кластеров), которое можно разбить на $j$ ($j = 1,2,...,Q$) фракций с разными размерами кластеров³³. Если пронумеровать эти фракции в порядке возрастания размеров кластеров, то 1-фракция будет содержать первичные мономеры, 2-фракция содержит ассоциации двух мономеров и т.д. В результате мы получим $Q$ фракций со следующими характеристиками:

Здесь $D$ – фрактальная размерность отдельного кластера и фрактальной среды в целом; ${{r}_{0}}$ и ${{m}_{0}}$ – радиус и масса первичных мономеров, из которых составлен фрактальный кластер; ${{\rho }_{0}} = {{3{{m}_{{\text{0}}}}} \mathord{\left/ {\vphantom {{3{{m}_{{\text{0}}}}} {4\pi \,r_{0}^{3}}}} \right. \kern-0em} {4\pi \,r_{0}^{3}}}$ – массовая плотность материала первичных ядер; ${{n}_{{{\text{0}}j}}} = {{({{{{R}_{j}}} \mathord{\left/ {\vphantom {{{{R}_{j}}} {{{r}_{0}}}}} \right. \kern-0em} {{{r}_{0}}}})}^{D}}$ – число первичных мономеров, входящих в состав $j$-го ФПА; ${{R}_{j}},$ $V_{j}^{{{\text{с l}}}} = {4 \mathord{\left/ {\vphantom {4 3}} \right. \kern-0em} 3}\pi r_{0}^{3}n_{{0j}}^{{3/D}},$ $m_{j}^{{{\text{с l}}}} = {{m}_{0}}{{({{{{R}_{j}}} \mathord{\left/ {\vphantom {{{{R}_{j}}} {{{r}_{0}}}}} \right. \kern-0em} {{{r}_{0}}}})}^{D}},$ ${{\text{v}}_{{pj}}}(r,t)$ – соответственно радиус гирации, объем, масса и гидродинамическая скорость ФПА $j$-го сорта; $\rho _{j}^{{{\text{с l}}}}(r,t) = m_{j}^{{{\text{с l}}}}N_{j}^{{{\text{с l}}}},$ $\tilde {\rho }_{j}^{{{\text{с l}}}}(r,t) = {{\rho }_{0}}{{\left( {{{R}_{j}}/{{r}_{0}}} \right)}^{{D - 3}}}$ – соответственно массовая плотность и истинная массовая плотность ФПА $j$-го сорта. $N_{j}^{{{\text{с l}}}}(r,t)$ – числовая плотность ФПА $j$-го сорта (число кластеров в единице объема дисковой среды).

В табл. 1 приведены характерные структурные параметры газовой составляющей диска.

Диффузия пылевых фрактальных агрегатов в газе

Перенос ФПА в газопылевом диске определяется их взаимодействием с молекулами несущего газа и двигающимися вместе с ними мелкодисперсными пылевыми частицами, и это взаимодействие имеет различный характер в зависимости от радиуса гирации кластеров и массовой плотности ${{\rho }_{{\text{g}}}}$ газа. В разреженной аэродисперсной среде, когда ${{l}_{{\text{g}}}} \gg {{R}_{j}}$ (где ${{l}_{{\text{g}}}}$ – длина свободного пробега атомов или молекул в газе), сила торможения движущегося пылевого агрегата создается в результате однократных столкновений с ним газовых частиц, что соответствует кинетическому режиму переноса этого агрегата в газовом потоке. В плотном газе, когда ${{l}_{{\text{g}}}} \ll {{R}_{j}},$ в каждый момент времени большое число частиц газа одновременно взаимодействует с кластером, и движение кластера соответствует диффузионному (гидродинамическому) режиму. В случае кинетического режима движения коэффициент диффузии $D_{{{\text{g}},pj}}^{{{\text{kin}}}}$ фрактальных кластеров в газе определяется формулой

а сила сопротивления (9) движению совокупности пылевых агрегатов в результате их взаимодействия с разреженной газовой фазой описывается обобщенным законом Эпстейна для агрегатов с размерами, меньшими или сравнимыми со средней длиной свободного пробега ${{l}_{{\text{g}}}}$ молекул газа. Здесь ${{c}_{{{\text{g}}s}}} = \sqrt {k{{{\kern 1pt} {{T}_{{\text{g}}}}} \mathord{\left/ {\vphantom {{{\kern 1pt} {{T}_{{\text{g}}}}} {{{m}_{{\text{g}}}}}}} \right. \kern-0em} {{{m}_{{\text{g}}}}}}} $ – изотермическая скорость звука в газе.

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

где ${{\eta }_{{\text{g}}}} = {{\rho }_{{\text{g}}}}{{\nu }_{{\text{g}}}},$ ${{\nu }_{{\text{g}}}} = ({{\sqrt {{8 \mathord{\left/ {\vphantom {8 \pi }} \right. \kern-0em} \pi }} \,} \mathord{\left/ {\vphantom {{\sqrt {{8 \mathord{\left/ {\vphantom {8 \pi }} \right. \kern-0em} \pi }} \,} 3}} \right. \kern-0em} 3})\,{{l}_{{\text{g}}}}{{c}_{{{\text{g}}s}}}$ – соответственно коэффициенты сдвиговой и кинематической вязкости (см. Чепмен, Каулинг, 1960); ${{l}_{{\text{g}}}} = {1 \mathord{\left/ {\vphantom {1 {{{\sigma }_{{\text{g}}}}}}} \right. \kern-0em} {{{\sigma }_{{\text{g}}}}}}{{n}_{{\text{g}}}};$ ${{\sigma }_{{\text{g}}}} = 2 \times {{10}^{{ - 15}}}$ см² – газокинетическое сечение столкновения частиц в газе; $\left| {{{w}_{{{\text{g}}p}}}} \right| = \left| {{{\text{v}}_{g}} - {{\text{v}}_{p}}} \right|;$ ${{C}_{{{\text{g}},j}}}({\text{R}}{{{\text{e}}}_{{{\text{g}},j}}})$ − эффективный коэффициент аэродинамического сопротивления движению пылевого кластера $j$-го сорта (сферы радиуса ${{R}_{j}}$) в газе, который в общем случае достаточно сложным образом зависит от числа Рейнольдса ${\text{R}}{{{\text{e}}}_{{{\text{g}},j}}} = 2{{R}_{j}}{{\left| {{{w}_{{{\text{g}}p}}}} \right|} \mathord{\left/ {\vphantom {{\left| {{{w}_{{{\text{g}}p}}}} \right|} {{{\nu }_{{\text{g}}}}}}} \right. \kern-0em} {{{\nu }_{{\text{g}}}}}}.$ В последнее время в астрофизической литературе получило достаточно широкое распространение следующее выражение для коэффициента ${{C}_{{{\text{g}},j}}}$ (Perets, Murray-Clay, 2011):

В данной работе, как указывалось выше, мы не будем учитывать процессы фрагментации ФПА при столкновениях, поэтому вполне естественно предположить, что относительные скорости ${{w}_{{gj}}}$ столкновения малы (столкновения кластеров с большими относительными скоростями сопровождаются, как известно, их разрушением). Тогда, при малых числах Рейнольдса ${{C}_{{{\text{g}},j}}}({\text{R}}{{{\text{e}}}_{{{\text{g}},j}}}) \cong {{24} \mathord{\left/ {\vphantom {{24} {{\text{R}}{{{\text{e}}}_{{{\text{g}},j}}}}}} \right. \kern-0em} {{\text{R}}{{{\text{e}}}_{{{\text{g}},j}}}}},$ коэффициент диффузии кластеров (при гидродинамическом режиме движения) принимает вид:

а соответствующая сила лобового сопротивления $F_{{{\text{g}},d}}^{{{\text{dif}}}}$ (9) задается законом Стокса

Заметим, что в протопланетном диске на расстоянии $r = 1\,\,{\text{а }}{\text{.}}\,{\text{е }}{\text{.}}$ длина свободного пробега частиц газа ${{l}_{{\text{g}}}} \approx 1\,\,{\text{с м }},$ а на расстоянии $r = 10\,\,{\text{а }}{\text{.}}\,{\text{е }}.$ ${{l}_{{\text{g}}}} \approx 1\,\,{\text{м }}.$ По этой причине пылевые агрегаты с размерами $1 - 100\,\,{\text{с м }}$ находятся в стоксовом режиме торможения на $1\,\,{\text{а }}{\text{.}}\,{\text{е }}.$ и в режиме Эпстейна на $10\,\,{\text{а }}{\text{.}}\,{\text{е }}.$

Формулы (10) и (12) для коэффициентов диффузии кластеров в диффузионном и кинетическом режимах удобно объединить и использовать для коэффициента диффузии пылевых образований в газовом потоке следующее соотношение:

которое переходит в (10) или (12) в пределе малых или больших чисел Кнудсена $K{{n}_{j}} = {{{{l}_{{\text{g}}}}} \mathord{\left/ {\vphantom {{{{l}_{{\text{g}}}}} {{{R}_{j}}}}} \right. \kern-0em} {{{R}_{j}}}} \times {{{{{10}}^{{15}}}} \mathord{\left/ {\vphantom {{{{{10}}^{{15}}}} {2{{n}_{{\text{g}}}}}}} \right. \kern-0em} {2{{n}_{{\text{g}}}}}}{{R}_{j}}$ соответственно. При этом сила взаимодействия (9) между совокупностью ФПА и газовой фазой диска принимает вид:

Важно отметить, что формула (17) при численном моделировании эволюции протопланетного облака позволяет учитывать плавный переход от кинетического режима к диффузионному режиму взаимодействия пылевых агрегатов с газом по мере увеличения как плотности газовых частиц ${{n}_{{\text{g}}}},$ так и радиусов гирации ${{R}_{j}}$ ФПА при их оседании к центральной плоскости диска, причем этот переход управляется параметром ${{n}_{{\text{g}}}}{{R}_{j}}.$

Ограничения на размеры пылевых фрактальных агрегатов, при которых справедлив закон Стокса

При написании формул (12) и (13) мы предположили, что соответствующие числа Рейнольдса малы. Выясним, при каких максимальных размерах пылевых агрегатов это допущение справедливо, на примере их свободного квазистационарного оседания в газе к центральной плоскости диска под влиянием $z$-компоненты силы тяготения Солнца ${{g}_{z}} = \Omega _{k}^{2}\,z$ (здесь ${{\Omega }_{k}}$ – кеплеровская угловая скорость на экваториальной плоскости диска). При указанных условиях уравнение движение $j$-кластера (7) сводится к виду

Отсюда для скорости гравитационного оседания (вдоль оси $z$) одиночного ФПА $j$-го сорта в газопылевом диске будем иметь

Определим число Рейнольдса формулой ${{\operatorname{Re} }_{{{\text{g}},j}}} = \;2{{R}_{j}}{{{{u}_{{zj}}}} \mathord{\left/ {\vphantom {{{{u}_{{zj}}}} {{{\nu }_{{\text{g}}}}}}} \right. \kern-0em} {{{\nu }_{{\text{g}}}}}};$ тогда условие на радиус гирации ФПА, при выполнении которого число Рейнольдса мало (${{\operatorname{Re} }_{{{\text{g}},j}}} \ll 1$), имеет вид: $1 \gg {\text{R}}{{{\text{e}}}_{{{\text{g}},j}}} = {{g}_{z}}\left( {4{{\rho }_{0}}{{r_{0}^{{3 - D}}} \mathord{\left/ {\vphantom {{r_{0}^{{3 - D}}} {9\eta _{{\text{g}}}^{2}}}} \right. \kern-0em} {9\eta _{{\text{g}}}^{2}}}} \right)R_{j}^{D},$ откуда

Это неравенство позволяет оценить максимальный размер стоксовских пылевых агрегатов. Из этой оценки (менее жесткой по сравнению с оценкой ${{R}_{j}} \ll {{\left( {{{9\eta _{{\text{g}}}^{2}} \mathord{\left/ {\vphantom {{9\eta _{{\text{g}}}^{2}} {4{{g}_{z}}{{\rho }_{{\text{g}}}}{{\rho }_{0}}}}} \right. \kern-0em} {4{{g}_{z}}{{\rho }_{{\text{g}}}}{{\rho }_{0}}}}} \right)}^{{1/3}}}$ для компактных тел) следует, в частности, что в гравитационном поле пушистые фрактальные агрегаты оседают значительно медленнее, чем компактные частицы той же массы.

Следует отметить, что в силу структурных особенностей кластеров сила сопротивления газовой среды более точно определяется рассеянием ее частиц на первичных ядрах фрактальных агрегатов (см. Михайлов, Власенко, 1995). Соответственно этому механизму рассеяния длина свободного пробега мономеров ${{\lambda }_{{\text{0}}}}$ внутри кластера (при средней плотности числа мономеров входящих в кластер ${{\langle n\rangle }_{{0,j}}} = {{{{n}_{{0,j}}}} \mathord{\left/ {\vphantom {{{{n}_{{0,j}}}} {V_{j}^{{{\text{с l}}}}}}} \right. \kern-0em} {V_{j}^{{{\text{с l}}}}}} = {{3R_{j}^{{D - 3}}} \mathord{\left/ {\vphantom {{3R_{j}^{{D - 3}}} {4\pi \,r_{0}^{D}}}} \right. \kern-0em} {4\pi \,r_{0}^{D}}}$) определяется соотношением

согласно которому для кластеров с фрактальной размерностью $D < {{2}^{\,}}$ выполняется условие ${{l}_{{j0}}} > {{R}_{j}}.$ Поскольку при таком способе описания условие применимости формулы Стокса для силы сопротивления движению фрактального агрегата имеет вид ${{l}_{{\text{g}}}} \ll {{R}_{j}}$ и ${{l}_{{j0}}} \ll {{R}_{j}},$ то из формулы (21) видно, что для фрактальных агрегатов с $D < 2$ это условие не выполняется, в силу чего формула Стокса для таких кластеров неприменима.

Массовую размерность $D$ определяют, как правило, на основе численного моделирования поведения кластера в гравитационном (или электрическом) поле с помощью “in situ”-методов процесса его сборки. Эти методы различаются специфическими деталями кластер-кластерной агрегации, к которым, в частности, относятся: характер движения кластеров (прямолинейное или броуновское), характер объединения кластеров в зависимости от вероятности слипания ${{\kappa }^{{\text{P}}}}$ при их взаимном касании, наличие или отсутствие полного реструктуринга (при котором кластеры связываются в трех точках), нарушение изотропии объединяемых кластеров, связанное с наведенным электрическим диполем во внешнем электрическом поле или наведенным магнитным моментом во внешнем магнитном поле, несферичность сталкивающихся кластеров, наличие вращательной диффузии ФПА, приводящей к захвату налетающего кластера краями образуемого пылевого агрегата (что способствует уменьшению его фрактальной размерности) и т.п. В табл. 2 представлены значения фрактальной размерности кластеров, образующихся в трехмерном пространстве при различных механизмах роста.

Итак, с учетом формулы (15), уравнения движения (6) и (7) для газа и совокупности пылевых гранул могут быть переписаны следующим образом:

Здесь

– время торможения пылевых гранул $j$-фракции в газовом потоке. Система уравнений (4), (5), (22) и (23) может быть использована для моделирования гидродинамической трехмерной неустойчивости во фрактальной газопылевой дисковой среде.

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

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

Гидродинамическая неустойчивость газопылевой среды в плоскости диска и способность твердых пылевых образований (кластеров, гранул, комков) гравитационно уплотняться управляются некоторым числом безразмерных параметров, среди которых ключевое значение имеет параметр металличности $\Gamma = \sum\nolimits_{j = 1}^Q {{{\Gamma }_{j}}} \,$ = = $\sum\nolimits_{j = 1}^Q {{{{{\Sigma }_{{pj}}}} \mathord{\left/ {\vphantom {{{{\Sigma }_{{pj}}}} {{{\Sigma }_{{\text{g}}}}}}} \right. \kern-0em} {{{\Sigma }_{{\text{g}}}}}}} = {{{{\Sigma }_{p}}} \mathord{\left/ {\vphantom {{{{\Sigma }_{p}}} {{{\Sigma }_{{\text{g}}}}}}} \right. \kern-0em} {{{\Sigma }_{{\text{g}}}}}},$ являющийся суммой отношений поверхностных плотностей пылевых фракций ${{\Sigma }_{{pj{\text{/g}}}}} = \sqrt {2\pi } {{\rho }_{{pj{\text{/g}}}}}{{H}_{{pj{\text{/g}}}}}$ и газовой ${{\Sigma }_{{\text{g}}}}$ компоненты дисковой среды. При значении параметра $\Gamma $ в протопланетном диске порядка ${{10}^{{ - 2}}},$ совокупная пылевая фаза динамически не влияет на движение газовой фазы, в то время как газ сильно воздействует на гидродинамику пыли. Здесь ${{H}_{{pj{\text{/g}}}}} = {{{{c}_{{pj{\text{/g}}}}}} \mathord{\left/ {\vphantom {{{{c}_{{pj{\text{/g}}}}}} {{{\Omega }_{K}}}}} \right. \kern-0em} {{{\Omega }_{K}}}}$ – высоты однородной атмосферы для отдельных пылевых фракций и газовой фазы, причем согласно (Youdin, Lithwick, 2007)

В эту формулу входит второй из ключевых параметров – безразмерное число Стокса для $j$-й пылевой фракции

где ${{\Omega }_{k}} = \sqrt {{{G{{M}_{{{\text{star}}}}}} \mathord{\left/ {\vphantom {{G{{M}_{{{\text{star}}}}}} {{{{\left| r \right|}}^{3}}}}} \right. \kern-0em} {{{{\left| r \right|}}^{3}}}}} $ – кеплеровская угловая скорость (равная характерной частоте колебаний звезды поперек плоскости диска), $r$ – центральный радиус-вектор,${{\tau }_{{{\text{stop}},j}}}$ – время торможения пылевых агрегатов $j$-й фракции в газовом потоке⁵⁵. Точное значение параметра ${\text{S}}{{{\text{t}}}_{j}}$ зависит как от размеров пылевых агрегатов, так и от их радиального местоположения (Chiang, Youdin, 2010). Время торможения ${{\tau }_{{{\text{stop}},j}}}$ пылевых агрегатов определяется формулой (24), а для плоского слоя (см. (29)).

В случае, если число Стокса для $j$-й пылевой фракции $1 \ll {\text{St}} \ll {\alpha }_{{\text{g}}}^{{ - {\text{1}}}},$ то

Здесь

– безразмерный коэффициент турбулентной интенсивности (см. Shakura, Sunyaev, 1973), учитывающий турбулизацию газового потока; $\nu _{{\text{g}}}^{{{\text{turb}}}}$ – коэффициент турбулентной вязкости газа. Аналогично ${{\alpha }_{{pj}}} = {{\nu _{{pj}}^{{{\text{turb}}}}} \mathord{\left/ {\vphantom {{\nu _{{pj}}^{{{\text{turb}}}}} {{{c}_{{s{\text{g}}}}}}}} \right. \kern-0em} {{{c}_{{s{\text{g}}}}}}}{{H}_{{pj}}},$ где $\nu _{{pj}}^{{{\text{turb}}}}$ – коэффициент вязкости для $j$-й фракции пыли в турбулентном потоке газа, причем параметры ${{\alpha }_{{pj}}}$ и ${{\alpha }_{{\text{g}}}}$ связаны соотношением:

В работе (Pinte и др., 2016) показано, что численное значение параметра ${{\alpha }_{p}}$ для пылевого слоя, вычисленное с учетом шкалы высот для пыли в субдиске, отличается от стандартной альфа-параметризации интенсивности α_g при глобальной турбулентности, имеющей место во всем протопланетном диске, равной ${{{\alpha }}_{{\text{g}}}} \approx {{10}^{{ - 2}}},$ тогда как ${{{\alpha }}_{p}} \approx 3 \times {{10}^{{ - 4}}}.$

Третий важный параметр ${{\zeta }_{j}} = {{c_{{pj}}^{2}} \mathord{\left/ {\vphantom {{c_{{pj}}^{2}} {c_{{s{\text{g}}}}^{2}}}} \right. \kern-0em} {c_{{s{\text{g}}}}^{2}}}$ – отношение дисперсий скоростей пылевой и газовой составляющих гетерогенного потока. Заметим, что в первом приближении, пульсации скоростей пылевых частиц связаны со свойствами турбулентности несущей газовой фазы, а не с межагрегатными столкновениями

Четвертый параметр

– параметр Тумре (Toomre, 1964) для газовой компоненты, который определяет классическую гравитационную неустойчивость (GI) газового диска в отсутствии пылевых частиц. При значении этого параметра ниже критического $Q_{{\text{g}}}^{{cr}} = 2$ возникает нестабильный диск, а при значении параметра ${{Q}_{{\text{g}}}}$ ниже единицы происходит фрагментация диска. Аналогичный параметр может быть введен и для пыли, причем ${{Q}_{p}} = {{\Gamma }^{{ - 1}}}{\kern 1pt} {{\zeta }^{{1/2}}}{{Q}_{{\text{g}}}}.$

Базовые уравнения гетерогенной механики в плоскости тонкого диска

Будем далее считать, что исходные уравнения двухфазной газопылевой гидродинамики осреднены по вертикальной толщине газопылевого диска. Будем также пренебрегать возможными осложнениями (например, сдвиговой неустойчивостью и турбулентностью), возникающими на масштабах, меньших шкалы высот для пыли ${{H}_{{pj}}} = {{{{c}_{{pj}}}} \mathord{\left/ {\vphantom {{{{c}_{{pj}}}} {{{\Omega }_{K}}.}}} \right. \kern-0em} {{{\Omega }_{K}}.}}$ Кроме этого будем предполагать, что дисперсия скорости газа ${{c}_{{\text{g}}}}$ и соответствующая шкала высот ${{H}_{{\text{g}}}} = {{c}_{{s{\text{g}}}}}/{{\Omega }_{K}}$ конечны, а слои пылевых фракций в субдиске находятся внутри газовой оболочки, так что $c_{{pj}}^{{}} < {{c}_{{sg}}}$ и ${{H}_{{pj}}} < {{H}_{{\text{g}}}}.$

Принимая локальную модель “обрезанной” коробки для диска (см. Goldreich, Lynden-Bell, 1965), мы сосредоточимся на рассмотрении окрестности точки $(r,\theta ) = ({{r}_{0}},{{\Omega }_{K}}t)$ в цилиндрических координатах. При этом радиальная и азимутальная координаты имеют вид: $(x,y) = (r - {{r}_{0}},\,\,{{r}_{0}}(\theta - {{\Omega }_{K}}t)),$ а кеплеровская частота ${{\Omega }_{K}} = \sqrt {{{G{{M}_{ \odot }}} \mathord{\left/ {\vphantom {{G{{M}_{ \odot }}} {r_{0}^{3}}}} \right. \kern-0em} {r_{0}^{3}}}} .$ Ограничимся стационарным фоном газопылевой среды с равномерным распределением плотностей пыли и газа (${{\Sigma }_{{pj0}}} = {\text{const}},$ ${{\Sigma }_{{{\text{g}}0}}} = {\text{const}}$). Тогда осредненные по высоте $z$ уравнения гетерогенной механики (5), (6), (22) и (23), используемые далее для изучения дисперсионных свойств возмущений в плоскости тонкого диска, принимают вид:

где ${{\nabla }_{ \bot }} \equiv {{i}_{x}}\partial /\partial x + {{i}_{y}}\partial /\partial y;$ $r = {{i}_{x}}x + {{i}_{y}}y$ – центральный радиус-вектор; ${{\Sigma }_{{pj{\text{/g}}}}} = \sqrt {2\pi } \,{{H}_{{pj{\text{/g}}}}}{{\rho }_{{pj{\text{/g}}}}},$ ${{p}_{{pj{\text{/}}s{\text{g}}}}} = c_{{pj{\text{/}}s{\text{g}}}}^{2}{{\Sigma }_{{pj{\text{/g}}}}}$ – соответственно поверхностные плотность и плоское давление отдельных пылевых фракций и газа в тонком диске; – время торможения пылевых гранул $j$-го сорта в газовой компоненте плоского диска.

Система уравнений (25)–(28) должна быть дополнена осредненным по высоте уравнением Пуассона

Здесь $\delta z$ – дельта-функция.

Линейные возмущения движения дисковой газопылевой среды

Для изучения динамики малых возмущений линеаризуем систему (25)–(28), для чего представим входящие в нее переменные в виде сумм равновесных и возмущенных величин:

слабо нарушающих невозмущенный однородный фон гидродинамических параметров ($\left| {\delta \chi } \right| \ll \left| {{{\chi }_{0}}} \right|$ для всех переменных $\chi = {{\chi }_{0}} + \delta \chi $). Для простоты ограничимся осевой симметрией задачи. При этом, по предположению, пылевые агрегаты считаются захваченными газовым потоком, и обе фазы, согласно кеплеровскому закону, двигаются с одинаковой скоростью ${{\text{v}}_{{(pj/g)}}}_{0} = - {3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0em} 2}{{\Omega }_{K}}\,x{{e}_{y}}.$

В случае пренебрежения нелинейными членами (квадратичными по амплитуде) и всеми пространственными производными от стационарных величин, малыми по сравнению с производными от возмущенных величин в системе (25)–(28), можно получить систему линеаризованных уравнений гетерогенной механики для описания малых колебаний вращающегося газопылевого диска. Будем далее считать, что радиальные возмущения осесимметричного газопылевого слоя имеют вид монохроматической волны, т.е. зависимость возмущенных величин $\delta \chi $ от радиальной координаты $x$ и времени $t$ пропорциональна $\exp ( - i\omega t + ikx),$ где $\omega $ – круговая частота монохроматической волны с длиной волны $\lambda = {{2\pi } \mathord{\left/ {\vphantom {{2\pi } {\left| k \right|}}} \right. \kern-0em} {\left| k \right|}};$ $k$ – волновое число; $\left| k \right| = \sqrt {k_{x}^{2} + k_{y}^{2}} $ (частота $\omega $ и волновое число $k$ могут быть комплексными величинами). В этом случае линеаризованная указанным выше способом система (25)–(29) перейдет в следующую систему алгебраических уравнений:

Заметим, что если учитывать влияние на гравитационный потенциал $\Phi $ толщин пылевого и газового слоев, то уравнение (37) следует переписать в виде (Shu, 1984):

где поправочный множитель $T(\left| k \right|H)$ = = ${1 \mathord{\left/ {\vphantom {1 {(1 + \left| k \right|H) < 1}}} \right. \kern-0em} {(1 + \left| k \right|H) < 1}}$ учитывает уменьшение самогравитации газопылевого слоя с ростом толщины диска для вертикально осредненной модели, причем в пределе тонкого диска для длинных волн, $\left| k \right|H \ll 1.$

Система (31)–(37) является исходной для дальнейшего анализа динамики малых линейных возмущений в рассматриваемой модели газопылевого диска с произвольным распределением величин ${{\Sigma }_{{pj,0}}},$ ${{\Sigma }_{{{\text{g0}}}}},$ ${{c}_{{pj{\text{/g}}}}},$ ${{\Phi }_{0}}.$ На основе этой системы уравнений возможно, вообще говоря, получить в явном виде дисперсионное соотношение для изучения неустойчивости тонкого слоя. Однако вывод этого соотношения в общем случае – процедура чрезвычайно трудоемкая. По этой причине ограничимся здесь (для простоты) рассмотрением двухжидкостной дисковой среды, состоящей из газа и совокупной пылевой фазы, полагая, что все пылевые фракции двигаются с одинаковой скоростью, ${{\text{v}}_{{pj}}} = {{\text{v}}_{p}}.$

Дисперсионное соотношение

Далее будем использовать обозначение $n = - i\omega .$ Очевидно, что неустойчивость возникает (возмущения растут), если величина $\operatorname{Re} [n]$ становится положительной. Приравнивая к нулю определитель системы алгебраических уравнений (31)–(37), после несложных преобразований получим следующее дисперсионное уравнение 5-го порядка по $n{\text{:}}$

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

С помощью уравнений, полученных на основании системы (25)–(30), возможно изучение как классических неустойчивостей (таких как гравитационная неустойчивость (GI), вековая гравитационная неустойчивость (SGI), потоковая неустойчивость (SI)), так и ряда новых неустойчивостей двухфазного газопылевого потока, например, резонансных неустойчивостей сопротивления (Resonant Drag Instabilities), некоторые из них могут иметь важное значение для образования планетезималей (Squire, Hopkins, 2018a; 2018b). При этом нахождение критериев возникновения тех или других гидродинамических неустойчивостей в газопылевом слое основано на существенном упрощении дисперсионного уравнения (38), возможного при использовании численных значений безразмерных параметров (таких, как ${\text{St}},$ ${{\Gamma }_{0}},$ $\zeta ,$ ${{Q}_{{p{\text{/g}}}}},$ ${{\alpha }_{{p{\text{/g}}}}},$ ${{H}_{{p{\text{/g}}}}}$ и др.), характерных для той или иной области протопланетного диска (см., например, Youdin, 2005а; 2005b; 2011; Youdin, Lithwick, 2007; Takahashi, Inutsuka, 2014; 2016; Carrera и др., 2015; 2017; Зиглина, Макалкин, 2016).

Оставляя строгий количественный анализ системы (38)–(43) на будущее численное моделирование, проиллюстрируем возможности этого подхода на простом примере аналитического определения условий возникновения гидродинамических неустойчивостей в тонком субдиске в случае, когда нет заметного обратного влияния пыли на орбитальное движение газа.

Линейные возмущения движения пыли в субдиске

Рассмотрим плоский пылевой субдиск, в котором возникает гидродинамическая неустойчивость. Характерная шкала расстояний в субдиске составляет приблизительно 5% от газовой шкалы высот ${{H}_{{\text{g}}}}.$ Будем далее предполагать, что ${{{{H}_{p}}} \mathord{\left/ {\vphantom {{{{H}_{p}}} {{{H}_{{\text{g}}}} < 1}}} \right. \kern-0em} {{{H}_{{\text{g}}}} < 1}}$ и газ движется устойчивым образом по круговой орбите с кеплеровой скоростью ${{\text{v}}_{{\text{g}}}} = - {3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0em} 2}{{\Omega }_{K}}x{{i}_{y}}.$ Заметим, что в реальных дисках из-за турбулентности газовой компоненты имеется пульсирующая компонента скорости газа, которая приводит к дисперсии скорости пыли ${{c}_{p}}$. Поскольку целью этого пункта является качественное объяснение гидродинамической неустойчивости, скорее, чем подробное моделирование динамики возмущений пылевой фазы в плоскости субдиска, будем для простоты предполагать, что жидкостной пылевой континуум является изотермическим, и таким образом, для пылевого давления справедливо соотношение ${{p}_{p}} = c_{p}^{2}{{\Sigma }_{p}},$ где дисперсия скорости пыли ${{c}_{p}}$ постоянна. В принятых предположениях система уравнений, описывающая возмущение гидродинамических параметров пылевого континуума, включает в себя линеаризованные указанным выше способом уравнения (26), (28) и (29):

где $k$ – радиальное волновое число. Условие существования нетривиальных решений для системы (44)–(47) приводит к следующему дисперсионному уравнению 3-го порядка относительно комплексной величины $n\,(k){\text{:}}$где – функция Варда (Ward, 2000).

Методом Кардана можно получить точное решение этого кубического уравнения. Однако это решение не приводит к простым формулам для различных показателей роста. Вместе с тем качественный анализ соотношения (48) возможен на основе рациональной аппроксимации отдельных его членов с учетом их предварительной оценки.

Вводя с этой целью безразмерные переменные (см. Youdin, 2005a)

перепишем соотношение (48) в следующем безразмерном виде: где

Ясно, что в случае ${{s}^{2}} < 0$ имеет место устойчивость решения (если частота $\omega $ – действительное число, то это означает, что волна распространяется по диску – система устойчива); в случае ${{s}^{2}} = 0$ имеет место нейтральность; для случая ${{s}^{2}} > 0$ имеем неустойчивость, если $Re\,\left[ s \right] > 0.$

Поскольку анализ корней дисперсионного уравнения (50) (особенно определение его неустойчивых мод) для “классического” случая компактности и одноразмерности пылевых агрегатов выполнен для всей области образования планетезималей в целом ряде работ (см., например, Youdin, 2005; 2011; Youdin, Lithwick, 2007; Chiang, Youdin, 2010; Shariff, Cuzzi, 2011; Takahashi, Inutsuka, 2014; 2016; Carrera и др., 2015; 2017), то мы не будем воспроизводить его здесь. Заметим только, что даже при массовом разнообразии параметров ${{c}_{p}},{{\alpha }_{p}},\Gamma $ и ${\text{St}}$ всегда можно найти волновое число κ такое, что система (44)–(47) будет неустойчива (Chiang, Youdin, 2010; Youdin, 2011).

Предварительный качественный анализ точного решения (найденного методом Кардана) кубического дисперсионного уравнения (50), выполненного для ряда комбинаций параметров кольцевых областей околосолнечного протопланетного диска (с радиальными расстояниями $r \in \left\{ {0.1\,\,{\text{a}}{\text{.}}\,{\text{e}}{\text{.}},\,\;1\,\,{\text{a}}{\text{.}}\,{\text{e}}{\text{.}},\;10\,\,{\text{a}}{\text{.}}\,{\text{e}}{\text{.}}} \right\},$ значениями параметров металличности $\Gamma \in $ {0.01, 0.05, 10} и турбулентной интенсивности $\alpha \in $ {${{10}^{{ - 2}}},{{10}^{{ - 4}}},{{10}^{{ - 5}}}$}), показал, что учет фрактальной размерности ($D \in $ {2, 2.11, 2.9}) и неодинаковости геометрических размеров (радиусов гирации $R \in $ {0.1 см, 10 см, 100 см}) пылевых агрегатов оказывает существенное влияние на условия, при которых возникают различные гидродинамические неустойчивости в пылевом слое. В последующих публикациях предполагается определение критических условий возникновения гидродинамических неустойчивостей на основе численной реализации уравнений (38)–(43), описывающих линейные возмущения гидродинамических параметров в плоскости вращающегося газопылевого диска.

Плотность состояний и обобщенные дифференциальные операторы для фрактальных сред

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

Для описания фрактальных сред с помощью дробно-интегральных моделей обычно используется так называемая функция плотности разрешенных состояний^{77${{c}_{n}}(D,R),$ описывающая, как плотно упакованы разрешенные состояния в $n$‑мерном ($n = 1,2,3$) евклидовом пространстве (см., например, Tarasov, 2005; 2010). При этом выражение ${{c}_{n}}(D,R)d{{V}_{n}}$ представляет собой число разрешенных мест в элементарном объеме $d{{V}_{n}}.$ Здесь $R = {r \mathord{\left/ {\vphantom {r {{{L}_{{{\text{hydr}}}}}}}} \right. \kern-0em} {{{L}_{{{\text{hydr}}}}}}}$ − безразмерная координата. Вид и свойства функции ${{c}_{n}}(D,R)$ определяются свойствами и симметриями фрактального распределения. Далее мы будем предполагать, что дисковая фрактальная среда обладает сферически симметричным распределением вещества в ней; в этом случае можно использовать следующее выражение для плотности состояний:}

являющиеся ядром интеграла дробного порядка $D$ (см. Tarasov, 2005), задающего интегрирование Рисса с точностью до числового множителя (Самко и др., 1987). Отметим, что если $D = 3,$ а $d = 2,$ то ${{c}_{3}}(3,2,R) = 1,$ если $d = 2,$ то ${{c}_{{\text{2}}}}(2,R) = 1.$

Отметим также, что в случае иных симметрий среды с дробной размерностью, для функции плотности разрешенных состояний ${{c}_{n}}(D,R)$ нужно использовать выражения для ядер интегралов дробного порядка, отвечающих этой симметрии, т.е. фигурирующих в дробно-интегральной модели данной фрактальной среды (в частности, в случае распределения вещества с фрактальной размерностью, соответствующей симметрии параллелепипеда, нужно использовать плотность состояний по Риману–Лиувиллю) (см., например, Самко и др., 1987).

Введем обобщенный дифференциальный оператор ${{\nabla }^{D}},$ с помощью которого удобно записать обобщенные гидродинамические уравнения для фрактальных сред. Используемый далее обобщенный оператор ${{\nabla }^{D}}$ в форме Рисса (см. Самко и др., 1987) имеет вид

Заметим, что правило почленного дифференцирования ${{\nabla }^{D}}(ab) = a{{\nabla }^{D}}b + \,b{{\nabla }^{D}}a$ для оператора набла ${{\nabla }^{D}}$ не справедливо; этот оператор удовлетворяет следующему правилу:

где

С помощью оператора ${{\nabla }^{D}}$ обобщенная дивергенция записывается в виде

а обобщенная (на пылевую фрактальную среду) полная производная по времени принимает вид

Если обозначить $\Xi (D,d) \equiv \frac{{{{2}^{{D - d}}}}}{{\sqrt \pi }}\frac{{{\Gamma (}{D \mathord{\left/ {\vphantom {D 2}} \right. \kern-0em} 2}{\text{)}}}}{{{\Gamma (}{d \mathord{\left/ {\vphantom {d 2}} \right. \kern-0em} 2}{\text{)}}}}$ и учесть, что $\nabla {{\left| R \right|}^{{d - 2}}} = (d - 2){{\left| R \right|}^{{d - 4}}}R,$ то

Базовые уравнения для моделирования эволюции фрактального газопылевого диска

Пусть газ и твердые частицы образуют два континуума, которые взаимодействуют друг с другом через аэродинамическое сопротивление. В общем случае необходимо вводить в рассмотрение объемные содержания отдельных фаз s, поскольку в отличие от гомогенной смеси, где каждая компонента может рассматриваться как занимающая весь элементарный объем пространства непрерывно и равномерно с другими компонентами, в гетерогенной системе каждая фракция занимает лишь часть элементарного объема. В частности, ${{s}_{{\text{g}}}} \equiv {{{{\rho }_{{\text{g}}}}} \mathord{\left/ {\vphantom {{{{\rho }_{{\text{g}}}}} {{{{\tilde {\rho }}}_{{\text{g}}}}}}} \right. \kern-0em} {{{{\tilde {\rho }}}_{{\text{g}}}}}},$ где ${{\tilde {\rho }}_{{\text{g}}}}$ – истинная плотность газа, равная отношению массы газа в элементарном объеме $\delta V$ к части этого объема $\delta {{V}_{{\text{g}}}},$ которую занимает масса газа. Объемное содержание ФПА $j$-го сорта определяется аналогично: $s_{j}^{{{\text{с l}}}}(r,t) = {{\rho _{j}^{{{\text{с l}}}}} \mathord{\left/ {\vphantom {{\rho _{j}^{{{\text{с l}}}}} {\tilde {\rho }_{j}^{{{\text{с l}}}}}}} \right. \kern-0em} {\tilde {\rho }_{j}^{{{\text{с l}}}}}},$ где ${{\rho _{j}^{{{\text{с l}}}}} \mathord{\left/ {\vphantom {{\rho _{j}^{{{\text{с l}}}}} {\tilde {\rho }_{j}^{{{\text{с l}}}}}}} \right. \kern-0em} {\tilde {\rho }_{j}^{{{\text{с l}}}}}}$ – соответственно массовая плотность и истинная массовая плотность. Тогда объемное содержание ${{s}_{{\text{g}}}}$ несущей фазы выражается в виде ${{s}_{{\text{g}}}} = 1 - s_{p}^{{{\text{с l}}}},$ где

В случае если суммарный гетерогенный континуум рассматривать в однодавленческом приближении, ${{P}_{\operatorname{g} }} = {{P}_{p}} = P({{\rho }_{g}},T),$ то система уравнений Эйлера, описывающая движение газа и фрактальной пылевой среды, имеет вид (см. Нигматулин, 1987; Tarasov, 2010):

Важно при этом отметить, что для исчерпывающего решения заявленной в этой статье проблемы необходимо в общем случае моделировать процессы коагуляции пылевых частиц и потоковую неустойчивость одновременно. Поэтому строгое решение задачи образования протопланетозималей в рамках рассматриваемого механизма включает, наряду с оценкой скоростей газопылевых кластеров, также и нахождение функции распределения (спектра) кластеров по размерам (массам) в произвольный момент времени $t$ в точке $r.$ Другими словами, необходимо совместное решение уравнений движения для фрактальной пылевой среды и обобщенного кинетического уравнения Смолуховского для пространственно неоднородной коагуляции. Заметим, что указанный подход использовался ранее в работах (Kolesnichenko, Marov, 2013; Колесниченко, 2017), в которых моделирование эволюции допланетного облака проводилось в рамках фрактальной газопылевой среды, использующей для учета фрактальности пылевого континуума дробные интегралы, порядок которых определялся массовой размерностью $D$ пылевых кластеров.

Работа выполнена в рамках госзаданий ИПМ им. М.В. Келдыша и ГЕОХИ им. В.И. Вернадского РАН при частичной поддержке грантами РФФИ № 17-02-00507, 18-01-00064 и Программы Президиума РАН № 12.