Астрономический вестник, 2020, T. 54, № 3, стр. 208-224

Параметрический метод моментов решения уравнения коагуляции Смолуховского в теории аккумуляции пылевых тел в допланетном диске

А. В. Колесниченко *

Институт прикладной математики им. М.В. Келдыша РАН
Москва, Россия

* E-mail: kolesn@keldysh.ru

Поступила в редакцию 26.11.2019
После доработки 26.12.2019
Принята к публикации 13.01.2020

Полный текст (PDF)

Аннотация

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

Ключевые слова: уравнение коагуляции Смолуховского, метод моментов, протопланетный диск, процессы коагуляции в пылевых средах

ВВЕДЕНИЕ

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

Создание адекватных космогонических моделей связано с изучением динамической и тепловой эволюции гетерогенного газопылевого вещества дифференциально вращающегося допланетного диска при учете турбулентных, магнитогидродинамических и радиационных процессов, протекающих с участием фазовых переходов, химических реакций и коагуляции (см., например, Сафронов, 1960; 1969; Goldreich, Ward, 1973; Nakamoto, Nakagawa, 1994; Weidenschilling, 1980; nakagawa и др., 1983; 1986; Dubrulle и др., 1995; Youdin, Shu, 2002; Колесниченко, 2004; 2005; Marov, Kolesnichenko, 2001; 2013). От пространственно-временного распределения гидро-термодинамических параметров дисковой среды зависят агрегатное состояние основных компонентов допланетного вещества, расположение их фронтов конденсации-сублимации, а также химический состав планет, их спутников, астероидов и комет.

К сожалению, все еще большое число проблем, связанных с данным направлением исследований, пока остается нерешенным. К ним, в первую очередь, относятся вопросы, связанные с ранними этапами эволюции Солнечной системы и причинами ее уникальности по сравнению с известными планетными системами у других звезд. В частности, условия, при которых начинается образование планет в протопланетном диске, все еще остаются плохо понятыми, несмотря на несколько десятилетий астрофизических исследований. В связи с этим, приоритетным по-прежнему остается ключевой вопрос космогонии – каким образом в протопланетном облаке происходит процесс роста твердых тел от пылевых частиц до крупных агломератов (с большой начальной массой порядка массы астероидов ~1015–1019 г и размерами в пределах 0.1–10 км) и от планетезималей до планет?

Важно отметить, что до последнего времени в большинстве теоретических моделей объединения пылевых частиц в допланетном диске изначально принималась компактная структура возникающих пылевых кластеров. Однако, как теперь стало ясно, растущие благодаря взаимным столкновениям частиц пылевые образования могут иметь весьма ажурную структуру и чрезвычайно низкую объемную плотность (см., например, Blum, 2004; Ormel и др., 2007; Suyama и др., 2008; Wada и др., 2008; Okuzumi и др., 2011; Suyama и др., 2012). Для подобных ворсистых агрегатов, имеющих по сравнению с компактными пылевыми частицами относительно большие геометрические поперечные сечения, меняется весь режим движения в газопылевой космической среде, в частности, из-за значительного изменения силы трения. Следовательно, для адекватного моделирования эволюции пылевых агрегатов в диске и, в конечном счете, механизма образования рыхлых протопланетезималей, нужно привлекать к рассмотрению их фрактальные свойства и внутреннюю структуру.

В отличие от ряда классических исследований (см., например, Сафронов, 1969; Weidenschilling, 1980; Nakagawa и др., 1981; 1986), в которых моделирование велось в рамках “обычной” сплошной среды, в работах (Колесниченко, 2000; Колесниченко, Маров, 2014а) предлагалось рассматривать совокупность пылевых агрегатов как особый тип сплошной среды – фрактальной среды, для которой существуют точки и области, не заполненные ее частицами. При этом гидродинамическое моделирование такой среды, обладающей нецелой массовой размерностью $D$, необходимо проводить, в общем случае, в рамках дробно-интегральной модели, использующей различные интегралы дробных порядков, для которых порядок дробного интегрирования определяется массовой размерностью пылевых кластеров, а сам интеграл интерпретируется как особый тип интеграла на фрактале с точностью до числового множителя (см. Tarasov, 2005; 2010).

Поскольку вид коэффициентов коагуляции пылевых кластеров (в частности, коэффициенты мономер-кластерной коагуляции и кластер-кластерной коагуляция) существенно изменяется во фрактальной среде (см. Kolesnichenko, 2001; Kolesnichenko, Marov, 2006; Колесниченко, Маров, 2014б), то это обстоятельствоусложняет такжеи решение кинетического уравнения Смолуховского, моделирующего процессы аккумуляции допланетных тел.

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

В качестве таких аппроксимирующих распределений в представленной работе предложено использовать кривые Пирсона (Pearson, Hartley, 1954). Решение поставленной задачи базируется на приближенных методах получения первых моментов функции распределения, которые позволяют найти центральные моменты, а затем и параметры кривых Пирсона. В работе проводится сравнение подобных методов и обсуждается практическая применимость каждого из них. На основе полученных моментовнайденысоотношения, позволяющие, в частности, оценить среднее значение, рассеяние, симметрию и островершинность вероятностного распределения искомых характеристик стохастической системы. Эти соотношения используются затем при подборе полуэмпирических моделей распределения по диаграмме Пирсона (Hahn, Shapiro, 1967). В качестве справочного материала приведены точные и полные формулы для основных типов кривых Пирсона. Полученные результаты применимы для любых распределений частиц по размерам, для которых найдены первые четыре начальных момента.

РЕШЕНИЕ УРАВНЕНИЯ СМОЛУХОВСКОГО МЕТОДОМ МОМЕНТОВ

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

Интегральная форма записи кинетического уравнения коагуляции Смолуховского (Smoluchowski, 1917) в этом случае записывается в виде:

(1)
$\begin{gathered} \frac{{\partial f({v},t)}}{{\partial t}} = \frac{1}{2}\int\limits_0^{v} {K({v} - w,w)f({v} - w,t)f(w,t){\text{d}}w} - \\ - \,\,f({v},t)\int\limits_0^\infty {K({v},w)f(w,t){\text{d}}w} , \\ \end{gathered} $
где $f({v},t)$ – непрерывная функция распределения (спектр) по объемам $v$ коагулирующего роя пылевых частиц; $K({v},w)$ – ядро коагуляции, описывающее вероятность столкновения и объединения частиц с объемами ${v}$ и $w.$ Коэффициент $K({v},w)$ является неотрицательной функцией, симметричной относительно своих аргументов, $K({v},w)$ = $ = K(w,{v}) \geqslant 0.$ Для решения интегро-дифференциального кинетического уравнения (1) необходимо задание начального условия
(2)
$f({v},0) = {{f}_{0}}({v}),$
а также выполнение условий $f({v},t) \to 0$ при ${v} \to 0$ и ${v} \to \infty .$

Точные методы решения уравнения коагуляции. Обсудим кратко некоторые методы решения кинетического уравнения коагуляции (1), которое по своей структуре является интегро-дифференциальным уравнением. Решение подобных уравнений представляет значительные математические трудности. Имеющиеся методы решения уравнения (1) для упрощения их анализа можно разбить на группы (см. рис. 1). Точные аналитические методы решений образуют небольшую группу, обзор их можно найти в работах (Волощук, Седунов, 1975; Попель и др., 1975; Волощук, 1984). При этом класс точных решений, основанных на применении интегрального преобразования Лапласа, может быть получен только для ядер коагуляции, являющихся линейными функциями по каждому из аргументов в отдельности, т.е. для ядер вида

(3)
$\begin{gathered} {{K}_{0}}({v},w) = {{G}_{0}},\,\,\,\,{{K}_{1}}({v},w) = {{G}_{1}}({v} + w), \\ {{K}_{2}}({v},w) = {{G}_{2}}{v}w. \\ \end{gathered} $
Рис. 1.

Методы решения уравнения коагуляции Смолуховского.

Наиболее теоретически продвинутыми к настоящему времени являются исследования процессов коагуляции для ядер типа ${{K}_{0}}({v},w) = {{G}_{0}},$ не зависящих от объемов коагулирующих частиц. Решение уравнения (1) с ядром ${{K}_{2}}({v},w) = {{G}_{2}}{v}w$ нельзя считать физически реализуемым, поскольку, начиная с некоторого момента времени, число частиц в системе становится отрицательным (см. Сафронов, 1969). В связи с исследованиями эволюции протопланетного газопылевого облака Сафроновым (1969) было получено аналитическое решение уравнения (1) с ядром, пропорциональным сумме объемов сталкивающихся частиц ${{K}_{1}}({v},w)$ = $ = {{G}_{1}}({v} + w),$ которое дает качественно правильный общий ход зависимости $K({v},w)$ от объемов тел. Однако, к сожалению, до сих пор не найдено ни одной дисперсной системы, для которой микрофизика коагуляционного процесса в точности приводила бы к ядрам подобного типа.

В монографии (Волощук, Седунов, 1975) рассмотрены решения кинетических уравнений с ядрами (3) и с начальным распределением в виде гамма-распределения

(4)
$\begin{gathered} f({v},0) = \\ = \left\{ \begin{gathered} {{N}_{0}}\frac{{{{\lambda }^{{\beta }}}}}{{\Gamma (\beta )}}{{{v}}^{{{\beta } - 1}}}{\text{exp}}( - \lambda {v}),\,\,\,\,{v} \geqslant 0,\,\,\,\,\lambda > 0,\,\,\,\,\beta > 0; \hfill \\ 0\,\,\,\,{\text{в}}\,\,\,\,{\text{остальных}}\,\,\,\,{\text{случаях,}} \hfill \\ \end{gathered} \right. \\ \end{gathered} $
где ${{N}_{0}}$ – суммарное количество частиц в начальный момент времени; $\Gamma (\beta ) = \int_0^\infty {{{x}^{{{\beta } - 1}}}{\text{exp}}( - x){\text{d}}x} $ – гамма-функция; в случае, когда $\beta $ – положительное целое число, $\Gamma (\beta ) = (\beta - 1)!$

Для начального условия (4) общее решение кинетического уравнения (1) с ядром ${{K}_{0}}({v},w) = {{G}_{0}}$ можно записать в виде

(5)
$f(X,{{\xi }_{0}}) = \frac{{{{{(1 - {{\xi }_{0}})}}^{2}}}}{{{v}{{\xi }_{0}}}}{\text{exp}}( - \beta X)\sum\limits_{k = 0}^\infty {\frac{{{{\beta }^{{k{\beta }}}}}}{{\Gamma (k\beta )}}} \xi _{0}^{k}{{X}^{{{\beta }k}}},$
где $X \equiv {{{v}{{N}_{0}}} \mathord{\left/ {\vphantom {{{v}{{N}_{0}}} S}} \right. \kern-0em} S};$ ${{\xi }_{0}} \equiv {{1 - N({{\tau }_{0}})} \mathord{\left/ {\vphantom {{1 - N({{\tau }_{0}})} {{{N}_{0}}}}} \right. \kern-0em} {{{N}_{0}}}};$ суммарное количество частиц $N({{\tau }_{0}})$ определяется равенством $N({{\tau }_{0}}) = {{N}_{0}}{{({{1 + {{N}_{0}}{{\tau }_{0}}} \mathord{\left/ {\vphantom {{1 + {{N}_{0}}{{\tau }_{0}}} 2}} \right. \kern-0em} 2})}^{{ - 1}}};$ ${{\tau }_{0}} \equiv {{G}_{0}}t;$ $S$ – количество пылевой фазы в единице объема.

Для случая $\beta = 1,$ соответствующего экспоненциальному начальному распределению, выражение (5) преобразуется к виду

(5.1)
$f(X,{{\xi }_{0}}) = \frac{X}{{v}}{{(1 - {{\xi }_{0}})}^{2}}{\text{exp}}\left[ { - X(1 - {{\xi }_{0}})} \right].$

Для $\beta = 2$ выражение (5) можно записать следующим образом

(5.2)
$f(X,{{\xi }_{0}}) = \frac{X}{{v}}\frac{{(1 - {{\xi }_{0}})}}{{\xi _{0}^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}}}{\text{exp}}( - 2X){\text{sh}}\left( {2X\xi _{0}^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}} \right).$

Общее решение уравнения (1) для случая ядра ${{K}_{1}}({v},w)$ = ${{G}_{1}}({v} + w)$ с начальным условием (4) записывается в виде

(6)
$\begin{gathered} f(X,{{\xi }_{0}}) = \frac{X}{{v}}{{\xi }_{1}}(1 - {{\xi }_{1}}){\text{exp}}\left[ { - (\beta + {{\xi }_{1}})X)} \right] \times \\ \times \,\,\sum\limits_{k = 0}^\infty {\frac{{{{\beta }^{{{\beta }(k + 1)}}}{{X}^{{{\beta } - 1 + k({\beta } + 1)}}}}}{{(k + 1)!\Gamma {\kern 1pt} {\kern 1pt} [\beta (k + 1)]}}} , \\ \end{gathered} $
где ${{\xi }_{1}} \equiv {{1 - N({{\tau }_{1}})} \mathord{\left/ {\vphantom {{1 - N({{\tau }_{1}})} {{{N}_{0}}}}} \right. \kern-0em} {{{N}_{0}}}};$ ${{\tau }_{1}} \equiv {{G}_{1}}t;$ суммарное количество частиц $N({{\tau }_{1}})$ определяется равенством $N({{\tau }_{1}}) = {{N}_{0}}\exp ( - S{{\tau }_{1}}).$

При $\beta = 1,$ т.е. для случая экспоненциального начального распределения, из (6) следует

(6.1)
$f(X,{{\xi }_{0}}) = \frac{{(1 - {{\xi }_{1}})}}{{{v}\xi _{1}^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}}}{\text{exp}}\left[ { - (\beta + {{\xi }_{1}})X)} \right]{{I}_{1}}\left( {2X\xi _{1}^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}} \right).$

Метод моментов. В тех практически важных случаях, когда полного знания закона распределения частиц не требуется, а достаточно лишь информации о поведении во времени нескольких целых моментов ${{m}_{j}} \equiv \int_0^\infty {{{{v}}^{j}}f({v},t){\text{d}}{v}} $ $(j = 0,1,2,...)$ функции распределения $f({v},t),$ задача существенно упрощается (см., например, Marov, Kolesnichenko, 2001; Wright и др., 2001; Estrada, Cuzzi, 2008; Mingzhou и др., 2008; Chenи др., 2014). Для ядер вида (3) она решается точно.

Рассмотрим, в качестве примера, ядро ${{K}_{1}}({v},w)$ = = ${{G}_{1}}({v} + w).$ Для нахождения решения воспользуемся методом моментов, смысл которого состоит в сведении уравнения (1) к бесконечной системе обыкновенных дифференциальных уравнений относительно моментов ${{m}_{j}}$$(j = 0,1,2,...).$ Для получения этой системы умножим обе части уравнения (1) на ${{{v}}^{j}}$ и проинтегрируем результат по ${v}$ в пределах от $0$ до $\infty ;$ в итоге получим следующую систему уравнений относительно моментов ${{m}_{j}}{\text{:}}$

(7)
$\begin{gathered} \frac{{\partial {{m}_{0}}}}{{\partial t}} \equiv \frac{{\partial N}}{{\partial t}} = - \frac{1}{2}\int\limits_0^\infty {\int\limits_0^\infty {K({v},w)f({v},t)f(w,t){\text{d}}wd{v}} } , \\ \frac{{\partial {{m}_{1}}}}{{\partial t}} \equiv \frac{{\partial S}}{{\partial t}} = \frac{1}{2}\int\limits_0^\infty {\int\limits_0^\infty {(w - {v})K({v},w)} } \times \\ \times \,\,f({v},t)f(w,t){\text{d}}w{\text{d}}{v} = 0, \\ \frac{{\partial {{m}_{j}}}}{{\partial t}} = \int\limits_0^\infty {\int\limits_0^\infty {\left[ {{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}{{{({v} + w)}}^{j}} - {{{v}}^{j}}} \right]K({v},w)} } \times \\ \times \,\,f({v},t)f(w,t){\text{d}}w{\text{d}}{v},\,\,\,\,(j = 2,3,...). \\ \end{gathered} $

Нулевой момент равен суммарному количеству частиц, а первый момент – количеству дисперсной фазы в единице объема:

(8)
$\begin{gathered} {{m}_{0}} \equiv N(t) = \int\limits_0^\infty {f({v},t){\text{d}}{v}} , \\ {{m}_{1}} = S(t) = \int\limits_0^\infty {{v}f({v},t){\text{d}}{v}} . \\ \end{gathered} $

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

Например, для ядра ${{K}_{1}}({v},w) \equiv {{G}_{1}}({v} + w)$ первые пять моментных уравнения (1) запишутся в виде:

(9)
$\begin{gathered} \frac{1}{{{{G}_{1}}}}\frac{{\partial {{m}_{0}}}}{{\partial t}} = - {{m}_{1}}{{m}_{0}},\,\,\,\,\frac{1}{{{{G}_{1}}}}\frac{{\partial {{m}_{1}}}}{{\partial t}} = 0, \\ \frac{1}{{{{G}_{1}}}}\frac{{\partial {{m}_{2}}}}{{\partial t}} = 2{{m}_{1}}{{m}_{2}},\,\,\,\,\frac{1}{{{{G}_{1}}}}\frac{{\partial {{m}_{3}}}}{{\partial t}} = 3\left( {{{m}_{1}}{{m}_{3}} + m_{2}^{2}} \right), \\ \frac{1}{{{{G}_{1}}}}\frac{{\partial {{m}_{4}}}}{{\partial t}} = 4{{m}_{1}}{{m}_{4}} + 10{{m}_{2}}{{m}_{3}}. \\ \end{gathered} $

Здесь предполагается, что дисперсная фаза не удаляется из рассматриваемого объема. Если ввести безразмерное время $\tau = {{G}_{1}}St,$ то точное решение системы (9) может быть записано следующим образом:

(10)
$\begin{gathered} {{m}_{1}} = S,\,\,\,\,{{m}_{2}} = {{m}_{2}}(0)exp2\tau , \\ {{m}_{3}} = \exp 3\tau \left[ {{{m}_{3}}(0) - \frac{{m_{2}^{2}(0)}}{S}(1 - \exp \tau )} \right], \\ {{m}_{4}} = \exp 4\tau \left\{ {10\frac{{{{m}_{2}}(0)}}{S}(exp\tau - 1)} \right. \times \\ \left. { \times \,\,\left[ {\frac{{m_{2}^{2}(0)}}{S}(exp\tau - 1) + {{m}_{3}}(0)} \right] + {{m}_{4}}(0)} \right\}. \\ \end{gathered} $

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

Поскольку в настоящее время не найдены аналитические методы решения кинетических уравнений (1) с произвольными ядрами, то многие исследования последних лет были посвящены разработке численных методов решения этих уравнений. Основное внимание в этих работах уделялось построению рациональных схем расчета уравнений (1), оценке точности получаемых результатов и развитию общих методик анализа результатов. Обзор по имеющимся реализациям численных методов можно найти, например, в работах (Волощук, 1984; Chan и др., 2006). Основная трудность при использовании численных методов решения уравнения (1) возникает из-за большого диапазона изменения частиц по размерам. Объемы, по которым производится интегрирование правой части уравнения (1), изменяются в общем случае в широких пределах, что затрудняет выбор рационального шага интегрирования. Для сужения области интегрирования применяют различные преобразования переменных и неравномерные сетки интегрирования (см. Волощук, 1984). Поэтой причине, несмотря на кажущуюся универсальность численных методов, они применяются сравнительно редко.

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

(11)
${{m}_{j}} \equiv \int\limits_0^\infty {{{{v}}^{j}}f({v},t){\text{d}}{v}} ,\,\,\,\,\left( {j = 0,1,2,...} \right).$
являющихся числовыми характеристиками случайной функции $f({v},t).$ Наиболее важной характеристикой центра распределения является среднее значение (математическое ожидание) $E({v}) \equiv {{m}_{1}}$ = $ = \int_0^\infty {{v}f({v},t){\text{d}}{v}} .$ Тогда $j$-й момент относительно математического ожидания (или центральный момент) определяется как ${{\mu }_{j}} \equiv E{{[{v} - {{m}_{1}}]}^{j}}$ = $ = \int_0^\infty {{{{({v} - {{m}_{1}})}}^{j}}f({v},t){\text{d}}{v}} .$ Первый центральный момент всегда равен нулю, ${{\mu }_{1}} = 0.$ Второй центральный момент (дисперсия) является показателем рассеяния и определяется соотношением

(12)
${{\mu }_{2}} \equiv {{\sigma }^{2}}({v}) \equiv E{{[{v} - {{m}_{1}}]}^{2}} = \int\limits_0^\infty {{{{({v} - {{m}_{1}})}}^{2}}f({v},t){\text{d}}{v}} .$

Легко показать, что ${{\mu }_{2}} = {{m}_{2}} - m_{1}^{2}.$ Третий момент относительно среднего, связанный с асимметрией распределения, определяется как

(13)
${{\mu }_{3}} \equiv E{{[{v} - {{m}_{1}}]}^{3}} = {{m}_{3}} - 3{{m}_{2}}{{m}_{1}} + 2m_{1}^{3}.$

Одновершинное распределение с ${{\mu }_{3}} < 0$ имеет левостороннюю (отрицательную) асимметрию. Если ${{\mu }_{3}} > 0,$ то распределение имеет правостороннюю (положительную) асимметрию. Для симметричного распределения ${{\mu }_{3}} = 0.$ Величина

(14)
$\sqrt {{{\beta }_{1}}} \equiv {{{{\mu }_{3}}} \mathord{\left/ {\vphantom {{{{\mu }_{3}}} {{{{({{\mu }_{2}})}}^{{{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0em} 2}}}}}}} \right. \kern-0em} {{{{({{\mu }_{2}})}}^{{{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0em} 2}}}}}}$
измеряет отношение асимметрии распределения к мере рассеяния. Этот нормированный показатель позволяет сравнивать асимметрию двух распределений, имеющих различный масштаб (см., например, Hahn, Shapiro, 1967). Наконец четвертый момент относительно среднего связан с островершинностью распределения и называется эксцессом. Он определяется как

(15)
${{\mu }_{4}} \equiv E{{({v} - {{m}_{1}})}^{4}} = {{m}_{4}} - 4{{m}_{3}}{{m}_{1}} + 6{{m}_{2}}m_{1}^{2} - 3m_{1}^{4}.$

Величина

(16)
${{\beta }_{2}} \equiv {{{{\mu }_{4}}} \mathord{\left/ {\vphantom {{{{\mu }_{4}}} {\mu _{2}^{2}}}} \right. \kern-0em} {\mu _{2}^{2}}}$
является относительным показателем эксцесса. Для нормального распределения, имеющего колоколообразную форму, значение ${{\beta }_{2}} = 3.0.$ Нормальное распределение обычно используется как стандарт, по которому сравнивается островершинность распределений.

Известно, что по первым четырем моментам распределения можно сделать параметрическую оценку самого распределения на основе диаграммы Пирсона, представленной на рис. 2 (см. Hahn, Shapiro, 1967). Далее будет продемонстрирована подобная методика.

Рис. 2.

Области в плоскости (${{\beta }_{1}},$ ${{\beta }_{2}}$) для различных распределений.

ПАРАМЕТРИЧЕСКИЕ МЕТОДЫ РЕШЕНИЯ УРАВНЕНИЯ КОАГУЛЯЦИИ

Многие исследования последних лет были направлены на разработку различных параметрических методов решения кинетических уравнений Смолуховского с произвольными ядрами (см., например, Волощук, Седунов, 1975; Lee, 1983; Волощук, 1984; Логунов, 1979; Frenklach, Harris, 1987; McGraw, 1997; Barrett, Jheeta, 1996; Синайский, 1997; Barrett, Webb, 1998; Колесниченко, 2000; 2001; Marov, Kolesnichenko, 2001; Frenklach, 2002; Marchisio, Fox, 2005; Chanи др., 2006; Mingzhou и др., 2008; Yu и др., 2008; 2011).

Параметрические методы основаны на предположении, что решение (точное или приближенное) принадлежит некоторому наперед заданному параметрическому классу функций, в котором ищется приближение. Теоретически распределение$f({v},t)$ всегда может быть представлено в виде бесконечного ряда11 по какой-либо полной ортогональной системе функций ${{\varphi }_{k}}({v}){\text{:}}$

(17)
$f({v},t) = \sum\limits_{k = 0}^\infty {{{\alpha }_{k}}(t)} {{\varphi }_{k}}({v}),$
что позволяет использовать это разложение как параметрическое представление распределения $f({v},t)$ с параметрами ${{\alpha }_{1}},{{\alpha }_{2}},....$ При этом для определения параметров ${{\alpha }_{1}},{{\alpha }_{2}},...,{{\alpha }_{k}}$ исходное кинетическое уравнение (1) следует заменить системой моментных уравнений (7), число которых совпадает с числом неизвестных коэффициентов. В результате при подстановке (17) в (7) можно получить следующую систему нелинейных уравнений относительно моментов:

(18)
$\begin{gathered} {{m}_{j}}\frac{{\partial {{m}_{j}}}}{{\partial t}} = \sum\limits_{k,l = 0} {{{\alpha }_{k}}{{\alpha }_{l}}{{\beta }_{{klj}}}} , \\ {{m}_{j}} = \sum\limits_{k = 0} {{{\alpha }_{k}}} \int\limits_0^\infty {\int\limits_0^\infty {{{{v}}^{j}}{{\varphi }_{k}}({v}){\text{d}}{v}} } , \\ {{\beta }_{{klj}}} = \int\limits_0^\infty {\int\limits_0^\infty {K({v},w)} } \times \\ \times \,\,\left[ {\frac{1}{2}{{{({v} + w)}}^{j}} - {{{v}}^{j}}} \right]{{\varphi }_{k}}({v}){{\varphi }_{l}}(w){\text{d}}{v}{\text{d}}w. \\ \end{gathered} $

В том случае, когда система функций ${{\varphi }_{j}}({v})$ ортогональна и нормирована на единицу, то при умножении исходного кинетического уравнения (1) на ${{\varphi }_{j}}({v})$ и последующем интегрировании по ${v},$ можно получить следующую систему уравнений

(19)
$\begin{gathered} \frac{{\partial {{\alpha }_{j}}}}{{\partial t}} = \sum\limits_{k,l = 0} {{{\alpha }_{k}}} {{\alpha }_{l}}\beta _{{klj}}^{ * },\,\,\,\,\beta _{{klj}}^{ * } = \int\limits_0^\infty {\int\limits_0^\infty {K({v},w)} } \times \hfill \\ \times \,\,\left[ {\frac{1}{2}{{\varphi }_{j}}({v} + w) - {{\varphi }_{j}}({v})} \right]{{\varphi }_{k}}({v}){{\varphi }_{l}}(w){\text{d}}{v}{\text{d}}w, \hfill \\ \end{gathered} $
которая, в отличие от системы (18), разрешена относительно производных по ${{\alpha }_{j}}.$ Коэффициенты ${{\beta }_{{klj}}}$ и $\beta _{{klj}}^{ * }$ в (18) и (19) постоянны и должны определятьсядо решениясоответствующей системы. К сожалению,системы уравнений (18) и (19), в силу большого объема и громоздкости необходимых вычислений, могут быть решены только численно с помощью ЭВМ.

Важно иметь в виду, что определение коэффициентов ${{\beta }_{{klj}}}$ существенно зависит от выбора базисных функций ${{\varphi }_{j}}({v}),$ в качестве которых могут быть использованы ряды Фурье, полиномы Эрмита, полиномы Чебышева, полиномы Лагерра и т.д. (см., например, Кендал, Стьюарт, 1966; Левин, 1974). Вместе с тем, исследования точности аппроксимации, которой можно достигнуть при использовании конечного числа членов разложения по полиномам Эрмита (так называемого ряда Грамма–Шарлье) показали следующее (см. Кендал, Стьюарт, 1966):

– Эти ряды могут вести себя нерегулярно, – сумма $k$членов может дать худшее приближение, чем сумма ($k - 1$) членов.

– Сумма конечного числа членов этих рядов может привести к отрицательным значениям на “хвосте” распределения.

– Использование конечных рядов дает удовлетворительные результаты только в случае распределений с умеренной асимметрией (${{\beta }_{1}} < 0.7$).

Этот негативный результат является довольно серьезным препятствием при применении нескольких начальных членов ряда (17) для аппроксимации решения кинетического уравнения (1), когда требуется определить не только моменты искомого распределения, но и само распределение.

По счастью, для некоторых ядер специальных видов параметрический метод решения уравнения Смолуховского может быть несколько видоизменен и эффективно использован. К таким ядрам относятся симметричные ядра со степенью однородности $\eta $, которые имеют вид:

(20)
$K({v},w) = G\sum\limits_{j = 0}^n {{{\alpha }_{j}}} \left( {{{{v}}^{{{\eta } - {{{\delta }}_{j}}}}}{{w}^{{{{{\delta }}_{j}}}}} + {{{v}}^{{{{{\delta }}_{j}}}}}{{w}^{{{\eta } - {{{\delta }}_{j}}}}}} \right),$
где $G$ – множитель, определяемый специфическими условиями, при которых происходит коагуляция;${{\alpha }_{j}}$и ${{\delta }_{j}}$ – заданные числа. К этому классу ядер относятся и ядра (3), при которых кинетическое уравнение допускает точные решения.

При подстановке (20) в (7) получается бесконечная система дифференциальных уравнений для моментов:

(21)
$\begin{gathered} \frac{{\partial {{m}_{0}}}}{{\partial t}} = \frac{{\partial N}}{{\partial t}} = - G\sum\limits_{j = 0}^n {{{\alpha }_{j}}} {{m}_{{{\eta } - {{{\delta }}_{j}}}}}{{m}_{{{{\delta }_{j}}}}}, \\ \frac{{\partial {{m}_{1}}}}{{\partial t}} = \frac{{\partial }}{{\partial t}} = 0, \\ .................................................. \\ \frac{{\partial {{m}_{p}}}}{{\partial t}} = \frac{G}{2}\sum\limits_{j = 0}^n {{{\alpha }_{j}}} \times \\ \times \,\,\sum\limits_{k = 1}^{p - 1} {C_{p}^{k}(} {{m}_{{{\eta } - {{{\delta }}_{j}} + p - k}}}{{m}_{{{{{\delta }}_{j}} + k}}} + {{m}_{{{\eta } - {{{\delta }}_{j}} + k}}}{{m}_{{{{{\delta }}_{j}} + p - k}}}), \\ (p = 2,3,...), \\ \end{gathered} $
где $C_{n}^{j} = {{n!} \mathord{\left/ {\vphantom {{n!} {j!(n - j)!}}} \right. \kern-0em} {j!(n - j)!}}$ – биномиальные коэффициенты.

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

Для уменьшения числа отыскиваемых параметров в качестве такого класса может быть выбрана совокупность известных распределений $f({v},t),$ включающая начальное распределение. Наиболее распространенным подходом является использование нормального распределения, которое дает приемлемое описание многих (хотя далеко не всех) реальных явлений. Аналогично, гамма-распределение и логарифмическое нормальное распределение могут быть использованы для описания случайных величин, ограниченных только с одной стороны (сверху, или снизу), в то время как бета-распределение – для описания случайных величин, ограниченных как сверху, так и снизу.

Хотя эти модели приводят к распределениям самой различной формы, все же и они не дают той степени обобщения, которая часто бывает необходима. Это иллюстрируется на рис. 2, где показаны области в плоскости $({{\beta }_{1}},{{\beta }_{2}})$ для различных распределений – нормального, бета-распределения (частный случай – равномерное распределение), гамма-распределения (частный случай – экспоненциальное распределение) и логарифмически нормального, где ${{\beta }_{1}}$ – квадрат нормированного показателя асимметрии, а ${{\beta }_{2}}$ – нормированный показатель островершинности. Сюда входит также $t$-распределение Стьюдента – симметричное распределение, сходящееся к нормальному, когда его параметр (число степеней свободы) произвольно увеличивается.

Пример 1. Предположим для простоты, что в результате процессов коагуляции распределение $f({v},t)$ остается в классе распределений, к которым принадлежит начальное распределение, а со временем меняются только его статистические параметры: среднее значение, дисперсия и т.п. В качестве начального распределения пылевых частиц по размеру (диаметру) ${\text{d}}({v})$ в газопылевом облаке, по аналогии с атмосферным аэрозолем, выберем двухпараметрическое логнормальное распределение.

Плотность вероятности логарифмически нормального закона зависит от среднего значения $\langle \ln {\text{d}}\rangle $ и показателя рассеяния (дисперсии) $\sigma _{L}^{2} \equiv \left\langle {{{{(\ln {\text{d}} - \langle \ln {\text{d}}\rangle )}}^{2}}} \right\rangle $логарифма диаметра ${\text{d:}}$

(22)
$\begin{gathered} f(d;\mu {\text{*}},{{\sigma }_{L}}) = \frac{N}{{{{\sigma }_{L}}{\text{d}}\sqrt {2\pi } }}\exp \left\{ { - \frac{{{{{(\ln {\text{d}} - \langle \ln {\text{d}}\rangle )}}^{2}}}}{{2\sigma _{L}^{2}}}} \right\} = \\ = \frac{N}{{{{\sigma }_{L}}{\text{d}}\sqrt {2\pi } }}\exp \left\{ { - \frac{{{{{\ln }}^{2}}({{\text{d}} \mathord{\left/ {\vphantom {{\text{d}} {\mu {\text{*}}}}} \right. \kern-0em} {\mu {\text{*}}}})}}{{2\sigma _{L}^{2}}}} \right\}. \\ \end{gathered} $

Медиана распределения определяется, как нетрудно убедиться из (22), соотношением $\mu * = \exp \langle \ln {\text{d}}\rangle ,$ а средние значения самого диаметра и его дисперсии соответственно равны

(23)
$\begin{gathered} \langle {\text{d}}\rangle = \exp \left( {{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}\sigma _{L}^{2} + \ln \mu {\text{*}}} \right), \\ {{\sigma }^{2}} \equiv \langle {{({\text{d}} - \langle {\text{d}}\rangle )}^{2}}\rangle = {{\langle {\text{d}}\rangle }^{2}}\left[ {\exp \sigma _{L}^{2} - 1} \right]. \\ \end{gathered} $

Используя эти соотношения, можно получить формулы для статистических параметров ($\sigma _{L}^{2}$ и $\mu {\text{*}}$) логнормального распределения (22) только через средний диаметр частиц $\langle {\text{d}}\rangle $ и относительную дисперсию ${{\beta }^{2}} \equiv {{\langle {{{({\text{d}} - \langle {\text{d}}\rangle )}}^{2}}\rangle } \mathord{\left/ {\vphantom {{\langle {{{({\text{d}} - \langle {\text{d}}\rangle )}}^{2}}\rangle } {{{{\langle {\text{d}}\rangle }}^{2}}}}} \right. \kern-0em} {{{{\langle {\text{d}}\rangle }}^{2}}}}$ их размера:

(24)
$\sigma _{L}^{2} = \ln (1 + {{\beta }^{2}}),\,\,\,\,\mu {\text{*}} = {{\langle d\rangle } \mathord{\left/ {\vphantom {{\langle d\rangle } {\sqrt {1 + {{\beta }^{2}}} }}} \right. \kern-0em} {\sqrt {1 + {{\beta }^{2}}} }}.$

Используем для этого формулу перехода $f({v}) = f\left[ {{\text{d}}({v})} \right]\left| {\frac{{\partial {\text{d}}({v})}}{{\partial {v}}}} \right|$ (справедливую для строго возрастающей функции ${v} = {v}({\text{d}})$ случайной величины ${\text{d}}$ (см. Hahn, Shapiro, 1967)) и распределение (22) для определения плотности начального распределения объема ${v} = (\pi {\text{/}}6){{{\text{d}}}^{3}}$ пылевых частиц; тогда

(25)
$\begin{gathered} f({v};{{\sigma }_{L}},\mu ) = \frac{N}{{3\sqrt {2\pi } {{\sigma }_{L}}{v}}}\exp \left[ { - \frac{{{{{\ln }}^{2}}({v}{\text{/}}\mu )}}{{18\sigma _{L}^{2}}}} \right], \\ \left( {\mu = (\pi {\text{/}}6){{\mu }^{{ * 3}}}} \right). \\ \end{gathered} $

Пусть теперь процесс коагуляции пылевых частиц в диске не меняет этого распределения, а со временем меняются только параметры $\mu (t)$ и $\sigma _{L}^{2}(t).$ Введем моменты логнормального распределения

(26)
${{m}_{p}}(t) = \frac{N}{{3\sqrt {2\pi } {{\sigma }_{L}}(t)}}\int\limits_0^\infty {{{{v}}^{{p - 1}}}} \exp \left\{ { - \frac{{{{{\ln }}^{2}}[{v}{\text{/}}\mu (t)]}}{{18\sigma _{L}^{2}(t)}}} \right\}{\text{d}}{v}.$

Согласно (Lee, 1983) для момента $p$-го порядка справедливо представление

(27)
$\begin{gathered} {{m}_{p}} = {{m}_{1}}{{\mu }^{{p - 1}}}\exp \left[ {{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0em} 2}({{p}^{2}} - 1)\sigma _{L}^{2}} \right], \\ {{m}_{1}} = S = {\text{const}}, \\ \end{gathered} $
позволяющее выразить дробные моменты, входящие в (21) через ${{m}_{1}},$ $\mu ,$ $\sigma _{L}^{2}.$ В итоге получим следующую параметрическую систему двух обыкновенных дифференциальных уравнений (число уравнений должно совпадать с числом неизвестных коэффициентов) для определения параметров $\mu (t),$ $\sigma _{L}^{2}(t)$ по заданным начальным значениям $\mu (0),$ $\sigma _{L}^{2}(0){\text{:}}$

(28)
$\begin{gathered} N \equiv {{m}_{0}} = S{{\mu }^{{ - 1}}}\exp \left( {{{ - 3} \mathord{\left/ {\vphantom {{ - 3} 2}} \right. \kern-0em} 2}\sigma _{L}^{2}} \right), \\ S \equiv {{m}_{1}} = {\text{const}},\,\,\,\,{{m}_{2}} = S\mu \exp \left( {{9 \mathord{\left/ {\vphantom {9 2}} \right. \kern-0em} 2}\sigma _{L}^{2}} \right); \\ \end{gathered} $
(29)
$\begin{gathered} \frac{{\partial {{m}_{0}}}}{{\partial t}} = - G\sum\limits_{j = 0}^K {{{\alpha }_{j}}{{m}_{{{\eta } - {{{\delta }}_{j}}}}}{{m}_{{{{{\delta }}_{j}}}}}} = - Gm_{1}^{2}{{\mu }^{{{\eta } - 2}}}\sum\limits_{j = 0}^n {{{\alpha }_{j}}} \times \\ \times \,\,\exp \left\{ {{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0em} 2}\left[ {\delta _{j}^{2} + {{{(\eta - {{\delta }_{j}})}}^{2}} - 2} \right]\sigma _{L}^{2}} \right\} = \\ = - G{{\mu }^{{\alpha }}}m_{0}^{2}\sum\limits_{j = 0}^K {{{\alpha }_{j}}} \exp \left\{ {{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0em} 2}\left[ {\delta _{j}^{2} + {{{(\eta - {{\delta }_{j}})}}^{2}}} \right]\sigma _{L}^{2}} \right\}, \\ \end{gathered} $
(30)
$\begin{gathered} \frac{{\partial {{m}_{2}}}}{{\partial t}} = 2G\sum\limits_{j = 0}^n {{{\alpha }_{j}}{{m}_{{{\eta } - {{{\delta }}_{j}} + 1}}}{{m}_{{{{{\delta }}_{j}} + 1}}}} = 2Gm_{1}^{2}{{\mu }^{{\alpha }}}\sum\limits_{j = 0}^n {{{\alpha }_{j}}} \times \\ \times \,\,\exp \left\{ {{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0em} 2}[{{{({{\delta }_{j}} + 1)}}^{2}} + {{{(\eta - {{\delta }_{j}} + 1)}}^{2}} - 2]\sigma _{L}^{2}} \right\} = \\ = 2G{{\mu }^{{{\eta } + 2}}}\sum\limits_{j = 0}^n {{{\alpha }_{j}}} \exp \left\{ {{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0em} 2}\left[ {\delta _{j}^{2} + {{{(\eta - {{\delta }_{j}})}}^{2}}} \right]\sigma _{L}^{2}} \right\}. \\ \end{gathered} $

Эта параметрическая система уравнений нелинейна и может быть решена только численно. Результаты подобного моделирования в связи с задачей осаждения пылевых частиц к центральной плоскости диска представлены в монографиях (Marov, Kolesnichenko, 2001; Колесниченко, Маров, 2014б). Здесь же отметим, что изменение во времени среднего числа частиц $N(t)$ можно оценить, предположив, что дисперсия $\sigma _{L}^{2}$ остается постоянной. В этом случае, ограничившись двумя первыми моментами, из (21) будем иметь

(31)
$\frac{{\partial N}}{{\partial t}} = - G{{\mu }^{{\alpha }}}{{N}^{2}}\sum\limits_{j = 0}^K {{{\alpha }_{j}}} \exp \left\{ {{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0em} 2}[\delta _{j}^{2} + {{{(\eta - {{\delta }_{j}})}}^{2}}]\sigma _{L}^{2}} \right\}.$

Решение этого уравнения, полученное при использовании начального условия $N(0) \equiv {{N}_{0}} = {S \mathord{\left/ {\vphantom {S {{\tilde {v}}(0)}}} \right. \kern-0em} {{\tilde {v}}(0)}},$ имеет вид

(32)
$\begin{gathered} N(t) = \frac{S}{{{\tilde {v}}(0)}}\frac{1}{{1 + qt}},\,\,\,\,q = G{{\mu }^{{\alpha }}}\frac{S}{{{\tilde {v}}(0)}}\sum\limits_{j = 0}^n {{{\alpha }_{j}}} \times \\ \times \,\,\exp \left\{ {{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0em} 2}\left[ {\delta _{j}^{2} + {{{(\eta - {{\delta }_{j}})}}^{2}}} \right]\sigma _{L}^{2}} \right\}, \\ \end{gathered} $
где ${\tilde {v}}(0) = (\pi {\text{/}}6){{\langle {\text{d}}\rangle }^{3}}$ = $\mu (0)\exp \left( {{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0em} 2}\sigma _{L}^{2}} \right)$ – начальное значение среднего объема. Отсюда, при использовании соотношения ${\tilde {v}}(t) = S{\text{/}}N,$ можно найти изменение во времени среднего объема частиц. Для относительно больших значений времени, когда $qt \gg 1,$ из (32) следует

(33)
$N(t) = {1 \mathord{\left/ {\vphantom {1 {G{{\mu }^{{\alpha }}}t}}} \right. \kern-0em} {G{{\mu }^{{\alpha }}}t}}\sum\limits_{j = 0}^n {{{\alpha }_{j}}} \exp \left\{ {{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0em} 2}\left[ {\delta _{j}^{2} + {{{(\eta - {{\delta }_{j}})}}^{2}}} \right]\sigma _{L}^{2}} \right\}.$

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

Пример 2. Пусть, например, исходная плотность распределения частиц по объемам является гамма-распределением (4). И пусть на начальном участке процесса коагуляции $f({v},t)$ остается в двухпараметрическом классе этого распределения, которое задается равенством

(34)
$p({v},0) = {{N}_{0}}\frac{{{{\lambda }^{{\beta }}}}}{{\Gamma (\beta )}}{{{v}}^{{{\beta } - 1}}}{\text{exp}}( - \lambda {v}),\,\,\,\,\lambda > 0,\,\,\,\,\beta > 0.$

Тогда $f({v},t)$ запишется в виде

(35)
$f({v},t) = \frac{S}{{\langle {v}\rangle }}p({v},t),\,\,\,\,\langle {v}\rangle = \int\limits_0^\infty {{v}p({v},t)d{v}} = \frac{\beta }{\lambda }.$

Определяя моменты $f({v},t),$ получим

(36)
$\begin{gathered} {{m}_{0}} = S\frac{\lambda }{\beta },\,\,\,\,{{m}_{2}} = {{m}_{0}}\frac{{(\beta + 1)\beta }}{{{{\lambda }^{2}}}}, \\ {{m}_{\alpha }} = \frac{{{{m}_{0}}}}{{{{\lambda }^{{\alpha }}}}}\frac{{\Gamma (\beta + \alpha )}}{{\Gamma (\beta )}},\,\,\,\,{{m}_{{1 + \alpha }}} = \frac{{{{m}_{0}}}}{{{{\lambda }^{{{\alpha } + 1}}}}}\frac{{(\beta + \alpha )\Gamma (\beta + \alpha )}}{{\Gamma (\beta )}}. \\ \end{gathered} $

При больших значениях $\beta $ и $\beta \gg \alpha ,$ раскрывая гамма-функцию по асимптотическому разложению Стирлинга, получим

(37)
${{\Gamma (\beta + \alpha )} \mathord{\left/ {\vphantom {{\Gamma (\beta + \alpha )} {\Gamma (\beta )}}} \right. \kern-0em} {\Gamma (\beta )}} \simeq A{{\beta }^{{\alpha }}},\,\,\,\,A = (1 + \alpha ){\text{exp}}( - \alpha ),$
что позволяет переписать третье и четвертое равенства в (36) в виде

(38)
${{m}_{{\alpha }}} = Am_{0}^{{1 - {\alpha }}}{{S}^{{\alpha }}},\,\,\,\,{{m}_{{1 + {\alpha }}}} = Am_{0}^{{ - {\alpha }}}{{S}^{{1 + {\alpha }}}}.$

Рассмотрим случай, когда ядро коагуляции определяется равенством

(39)
$K({v},w) = G{{{v}}^{{{\eta /2}}}}{{w}^{{{\eta /}2}}}.$

Соответствующие ему первые три уравнения для моментов запишем в виде

(40)
$\frac{{\partial {{m}_{0}}}}{{\partial t}} = - \frac{1}{2}Gm_{{{\eta /2}}}^{2},\,\,\,\,\frac{{\partial {{m}_{1}}}}{{\partial t}} = 0,\,\,\,\,\frac{{\partial {{m}_{2}}}}{{\partial t}} = 2Gm_{{1 + {\eta /}2}}^{2}.$

Учитывая (38), получим

(41)
$\begin{gathered} \frac{{\partial {{m}_{0}}}}{{\partial t}} = - \frac{1}{2}G{{A}^{2}}{{S}^{{\eta }}}m_{0}^{{2(1 + {\eta /}2)}},\,\,\,\,{{m}_{1}} = S, \\ \frac{{\partial {{m}_{2}}}}{{\partial t}} = - 2G{{A}^{2}}{{S}^{{2(1 + {\eta /}2)}}}m_{0}^{{\eta }}. \\ \end{gathered} $

Дополняя эти уравнения начальными условиями и решая их, получим

(42)
$\begin{gathered} \frac{{{{m}_{0}}(\tau )}}{{{{m}_{0}}(0)}} = {{\left\{ {1 + \frac{{1 - \eta }}{2}{{{[{{m}_{0}}(0)]}}^{{1 - {\eta }}}}\tau } \right\}}^{{1{\text{/}}({\eta } - 1)}}}, \\ \tau = G{{A}^{2}}{{S}^{{\eta }}}t, \\ \end{gathered} $
(43)
$\frac{{{{m}_{2}}(\tau )}}{{{{m}_{2}}(0)}} = 1 + \frac{{2{{S}^{2}}}}{{{{m}_{2}}(0){{m}_{0}}(0)}}\left[ {\frac{{{{m}_{0}}(0)}}{{{{m}_{0}}(\tau )}} - 1} \right].$

НЕПАРАМЕТРИЧЕСКИЕ МЕТОДЫ РЕШЕНИЯ УРАВНЕНИЯ КОАГУЛЯЦИИ

Параметрические методы регуляризации системы моментных уравнений (21), несмотря на математическую простоту получения решения, основаны на очень сильном исходном допущении об a priori известном виде искомого распределения, которое обычно задают “волевым” способом. Этот изъян в выборе доопределяющих соотношений может быть ликвидирован, если воспользоваться непараметрическим методом интерполяции для определения связей между целыми и дробными моментами на временном интервале $\left[ {t,t + \Delta t} \right].$ С этой целью перейдем к безразмерным переменным при помощи нормирования всех моментов на их значения в начале интервала и введем интерполяционный полином Лагранжа $L_{j}^{n}(x)$ для оценки дробного момента ${{m}_{{i + \alpha }}}$ в следующем виде (см., например, Логунов, 1979; Press и др., 1992; Marov, Kolesnichenko, 2001)

(44)
$\begin{gathered} m_{{i + {\alpha }}}^{'}(t) = {{\sum\limits_{j = k}^{k + n} {\left[ {m_{j}^{'}(t)} \right]} }^{{L_{j}^{{(n)}}}}},\,\,\,\,k \leqslant t \leqslant k + n, \\ m_{{\gamma }}^{'}(t) \equiv \frac{{{{m}_{{\gamma }}}(t)}}{{{{m}_{{\gamma }}}(0)}}, \\ \end{gathered} $
где

(45)
$\begin{gathered} L_{j}^{n}(x) = \frac{1}{{n!}}\prod\nolimits_{n + 1} {\frac{{{{{( - 1)}}^{{n - j}}}C_{n}^{j}}}{{x - j}}} ,\,\,\,\,C_{n}^{j} = \frac{{n!}}{{j!(n - j)!}}, \\ \prod\nolimits_{n + 1} {(x) = x(x - 1)(x - 2)...(x - n)} , \\ x = i + \alpha - j. \\ \end{gathered} $

Интерполяционная формула (44) справедлива для любого времени $t$ и обеспечивает совпадение оценки интерполируемого момента с его точным значением в любой момент времени, если $\alpha $ – целое, а при $t = 0,$ если $\alpha $ – дробное число (см. Демидович, Марон, 1963). Определяятеперь моменты распределения из соответствующей системы дифференциальных уравнений (при выбранном начальном распределении) можно по первым пяти моментам определить вид искомого распределения частиц по размерам с помощью диаграммы Пирсона (см. Hahn, Shapiro, 1967).

Отметим, что для оценки точности интерполяции по формуле (44) в литературе были проведены многочисленные численные расчеты дробных моментов гамма- и логнормальных распределений. Результаты показали, что относительная ошибка интерполяции монотонно убывает с ростом параметра $n$, при увеличении дисперсий исследуемых распределений, а также при увеличении порядка интерполируемого момента. Так, относительная ошибка интерполяции для начальных дробных моментов гамма-распределения зависит только от его параметра формы $\beta $ и имеет порядок ${{10}^{{ - 2}}}$ при $\beta = 2$ и порядок ${{10}^{{ - 3}}}$ при $\beta = 5.$ С ростом $\beta $ ошибка монотонно убывает. Знак ошибки в общем случае зависит от $\beta $ и порядка интерполируемого момента. Для логнормального распределения соотношение (44) становится точным при $n \geqslant 2.$

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

Пусть $K({v},w)$ – симметричная функция двух переменных со степенью однородности $\eta ,$ которую надо аппроксимировать рядом вида (20). Учитывая условие однородности, запишем $K({v},w)$ в виде (Логинов, 1979)

(46)
$\begin{gathered} K({v},w) = {{{v}}^{{\eta }}}K(1,X) = {{{v}}^{{\eta }}}\sum\limits_{j = 0}^n {{{\beta }_{j}}\left( {{{X}^{{{{{\delta }}_{j}}}}} + {{X}^{{{\eta } - {{{\delta }}_{j}}}}}} \right)} , \\ X \equiv \frac{w}{{v}}. \\ \end{gathered} $

Неизвестные коэффициенты ${{\beta }_{j}}$ определим из условия выполнения равенства между правой и левой частью (46) при фиксированных значениях $X,$ которые являются точками привязки. Чтобы при последующем интегрировании разложения (46) не появлялись моменты порядка выше степени однородности $\eta ,$ коэффициенты ${{\delta }_{j}}$ можно определить равенством

(47)
${{\delta }_{j}} = \eta j{{(n + 1)}^{{ - 1}}}.$

С учетом введенных обозначений систему уравнений для нахождения коэффициентов ${{\beta }_{j}}$ можно переписать в виде

(48)
$K(1,{{X}_{k}}) = \sum\limits_{j = 0}^n {{{\beta }_{j}}\left( {X_{k}^{{{{{\delta }}_{j}}}} + X_{k}^{{{\eta } - {{{\delta }}_{j}}}}} \right)} ,\,\,\,\,k = 1,2,3,.....$

При умножении правой и левой частей (48) на ${{{v}}^{{\eta }}},$ получим следующее: при $n = 0$ аппроксимация ядра $K({v},w)$ отыскивается в классе функций

(49)
$K({v},w) = \beta _{0}^{{(0)}}({{{v}}^{{\eta }}} + {{w}^{{\eta }}});$
при $n = 1$ – в классе функций
(50)
$K({v},w) = \beta _{0}^{{(1)}}({{{v}}^{{\eta }}} + {{w}^{{\eta }}}) + \beta _{1}^{{(1)}}{{{v}}^{{{\eta /2}}}}{{w}^{{{\eta /2}}}};$
при $n = 2$ – в классе функций
(51)
и т.д.

Таким образом, решая уравнения (48) относительно коэффициентов ${{\beta }_{j}}$ можно получить: при $n = 0$

(52)
$\beta _{0}^{{(0)}} = {{K(1,{{X}_{1}})} \mathord{\left/ {\vphantom {{K(1,{{X}_{1}})} {\left( {1 + X_{1}^{{\eta }}} \right)}}} \right. \kern-0em} {\left( {1 + X_{1}^{{\eta }}} \right)}};$
при $n = 1$
(53)
$\begin{gathered} \beta _{0}^{{(1)}} = {{\left( {K(1,{{X}_{1}})X_{2}^{{{\eta /2}}} - K(1,{{X}_{2}})X_{1}^{{{\eta /2}}}} \right)} \mathord{\left/ {\vphantom {{\left( {K(1,{{X}_{1}})X_{2}^{{{\eta /2}}} - K(1,{{X}_{2}})X_{1}^{{{\eta /2}}}} \right)} {\left( {X_{2}^{{{\eta /2}}}(1 + X_{1}^{{\eta }}) - X_{1}^{{{\eta /2}}}(1 + X_{2}^{{\eta }})} \right)}}} \right. \kern-0em} {\left( {X_{2}^{{{\eta /2}}}(1 + X_{1}^{{\eta }}) - X_{1}^{{{\eta /2}}}(1 + X_{2}^{{\eta }})} \right)}}, \\ \beta _{1}^{{(1)}} = {{\left( {K(1,{{X}_{2}})(1 + X_{1}^{{\eta }}) - K(1,{{X}_{1}})(1 + X_{2}^{{\eta }})} \right)} \mathord{\left/ {\vphantom {{\left( {K(1,{{X}_{2}})(1 + X_{1}^{{\eta }}) - K(1,{{X}_{1}})(1 + X_{2}^{{\eta }})} \right)} {\left( {X_{2}^{{{\eta /2}}}(1 + X_{1}^{{\eta }}) - X_{1}^{{{\eta /2}}}(1 + X_{2}^{{\eta }})} \right)}}} \right. \kern-0em} {\left( {X_{2}^{{{\eta /2}}}(1 + X_{1}^{{\eta }}) - X_{1}^{{{\eta /2}}}(1 + X_{2}^{{\eta }})} \right)}}; \\ \end{gathered} $
при $n = 2$

(54)

Аппроксимацию ядра $K({v},w)$ называют $m$-точечной, если правая и левая части соотношения (46) совпадают в $m$ точках при изменении параметров ${{X}_{i}}$ от $0$ до $\infty .$ Когда выбрана $m$-точечная схема привязки, то последовательность точек ${{X}_{i}}$ для системы (48) определяется следующим рекуррентным соотношением:

(55)
$\begin{gathered} (m - 1)\left[ {{{{({{X}_{i}} + 1)}}^{{ - 1}}} - {{{({{X}_{{i + 1}}} + 1)}}^{{ - 1}}}} \right] = 1, \\ i = 2,3,....,\,\,\,\,{{X}_{1}} = 0. \\ \end{gathered} $

Рассчитанные по формуле (55) значения величин ${{X}_{i}}$ при различном числе точек привязки приведены в табл. 1.

Таблица 1.  

Значения величин ${{X}_{i}}$ при различном числе точек привязки (Логинов, 1979)

X m = 2 m = 3 m = 4 m = 5 m = 6 m = 7 m = 8 m = 9
X0 0 0 0 0 0 0 0 0
X1   1 1/2 1/3 1/4 1/5 1/6 1/7
X2       1 2/3 1/2 2/5 1/3
X3           1 3/4 3/5

АППРОКСИМАЦИЯ РЕШЕНИЯ УРАВНЕНИЯ СМОЛУХОВСКОГО ЭМПИРИЧЕСКИМИ РАСПРЕДЕЛЕНИЯМИ

Подход Пирсона. Рассмотрим теперь подход Пирсона к описанию распределений, основанный на отыскании семейства кривых, при помощи которого можно удовлетворительно представить встречающиеся на практике распределения. На рис. 3 представлены графики для определения типа кривой Пирсона в зависимости от параметра ${{\beta }_{1}},$ измеряющего отношение асимметрии распределения к мере рассеяния, и параметра ${{\beta }_{2}},$ являющегося относительным показателем эксцесса. Эти графики заимствованы из работ (Крамер, 1948; Elderton, 1953; Кендал, Стьюарт, 1966).

Рис. 3.

Графики для определения типа кривой Пирсона в зависимости от ${{\beta }_{1}}$ и ${{\beta }_{2}}.$

Пирсон предложил для описания статистического распределения $f(x)$ случайной величины $x$ использовать решения дифференциального уравнения

(56)
$\frac{{{\text{d}}f}}{{{\text{d}}x}} = - \frac{{(x - a)f}}{{{{b}_{0}} + {{b}_{1}}x + {{b}_{2}}{{x}^{2}}}}.$

Семейство плотностей, определяемых этой формулой, известно под названием “семейства распределений Пирсона”. Если иметь в виду унимодальные распределения, то представляется интересным изучить тот класс плотностей, которые: (а) имеют единственную моду, т.е. ${{{\text{d}}f} \mathord{\left/ {\vphantom {{{\text{d}}f} {{\text{d}}x}}} \right. \kern-0em} {{\text{d}}x}} = 0$ в некоторой точке $x = a,$ где $a$ – мода; (б) имеют гладкое соприкосновение с осью $x$ на концах интервала, где сосредоточено распределение, т.е. ${{{\text{d}}f} \mathord{\left/ {\vphantom {{{\text{d}}f} {{\text{d}}x}}} \right. \kern-0em} {{\text{d}}x}} = 0,$ когда $f = 0.$ Нетрудно видеть, что решения уравнения (56) удовлетворяют этим условиям. Следует заметить, что среди распределений семейства (56) существуют и такие распределения, которые имеют J- и U-образную форму.

Прежде, чем переходить к нахождению семейства точных решений уравнения (56), рассмотрим некоторые свойства, присущие этому семейству в целом. Имеем: ${{x}^{n}}({{b}_{0}} + {{b}_{1}}x + {{b}_{2}}{{x}^{2}})\frac{{{\text{d}}f}}{{{\text{d}}x}}$ = ${{x}^{n}}(x - a).$ Интегрируя левую часть этого уравнения по частям и предполагая, что полученные интегралы существуют, получим

(57)
$\begin{gathered} \left[ {{{x}^{n}}({{b}_{0}} + {{b}_{1}}x + {{b}_{2}}{{x}^{2}})f} \right]_{{ - \infty }}^{\infty } - \\ - \,\,\int\limits_{ - \infty }^\infty {\left( {n{{b}_{0}}{{x}^{{n - 1}}} + (n + 1){{b}_{1}}{{x}^{n}} + (n + 2){{b}_{2}}{{x}^{{n + 1}}}} \right)} f{\text{d}}x = \\ = \int\limits_{ - \infty }^\infty {{{x}^{{n + 1}}}f{\text{d}}x - a} \int\limits_{ - \infty }^\infty {{{x}^{n}}f{\text{d}}x} . \\ \end{gathered} $

Предположим теперь, что выражение в квадратных скобках обращается в нуль на концах распределения $\mathop {{\text{lim}}}\limits_{f \to \pm \infty } {{x}^{{n + 2}}}f \to 0,$ или что распределение имеет бесконечный диапазон. Тогда, пользуясь обозначениями (11) для моментов относительно нуля, будем иметь

(58)
$\begin{gathered} n{{b}_{0}}{{m}_{{n - 1}}} + \{ (n + 1){{b}_{1}} - a\} {{m}_{n}} + \\ + \,\,\{ (n + 2){{b}_{2}} + 1\} {{m}_{{n + 1}}} = 0. \\ \end{gathered} $

Это уравнение позволяет выразить четыре константы $a,$ ${{b}_{0}},$ ${{b}_{1}}$ и ${{b}_{2}}$ через моменты ${{m}_{1}},$ ${{m}_{2}},$ ${{m}_{3}},$ ${{m}_{4}},$ или через три момента (13), взятых относительно среднего: ${{\mu }_{2}},$ ${{\mu }_{3}},$ ${{\mu }_{4}}.$ Полагая в (58) последовательно $n = 0,1,2,3,$ можно найти следующие соотношения:

(59)
$\begin{gathered} a = - \frac{{{{\mu }_{3}}\left( {{{\mu }_{4}} + 3\mu _{2}^{2}} \right)}}{A} = - \frac{{\sqrt {{{\mu }_{2}}} \sqrt {{{\beta }_{1}}} ({{\beta }_{2}} + 3)}}{{A{\kern 1pt} '}}, \\ {{b}_{0}} = - \frac{{{{\mu }_{2}}\left( {4{{\mu }_{2}}{{\mu }_{4}} - 3\mu _{3}^{2}} \right)}}{A} = - \frac{{{{\mu }_{2}}(4{{\beta }_{2}} - 3{{\beta }_{1}})}}{{A{\kern 1pt} {\kern 1pt} '}}, \\ {{b}_{1}} = - \frac{{{{\mu }_{3}}\left( {{{\mu }_{4}} + 3\mu _{2}^{2}} \right)}}{A} = - \frac{{\sqrt {{{\mu }_{2}}} \sqrt {{{\beta }_{1}}} ({{\beta }_{2}} + 3)}}{{A{\kern 1pt} {\kern 1pt} '}}, \\ {{b}_{2}} = - \frac{{\left( {4{{\mu }_{2}}{{\mu }_{4}} - 3\mu _{3}^{2} - 6\mu _{2}^{3}} \right)}}{A} = - \frac{{(2{{\beta }_{2}} - 3{{\beta }_{1}} - 6)}}{{A{\kern 1pt} '}}. \\ \end{gathered} $
где
(60)
$\begin{gathered} A \equiv 10{{\mu }_{2}}{{\mu }_{4}} - 18\mu _{2}^{3} - 12\mu _{3}^{2}, \\ A{\kern 1pt} ' \equiv 10{{\beta }_{2}} - 18 - 12{{\beta }_{1}}, \\ \end{gathered} $
а параметры ${{\beta }_{1}}$ и ${{\beta }_{2}}$ определяются формулами (14) и (16). В этих формулах нулевое значение принято в качестве среднего, так что распределения семейства (56) полностью определяются своими четырьмя первыми моментами.

Семейство кривых Пирсона, удовлетворяющее уравнению (56) приведено в табл. 2.

Таблица 2.  

Плотности распределений из семейства Пирсона

Тип уравнения Начало отчета для $x$ Область определения
I. Бета-распределение
$f(x) = {{f}_{0}}{{\left( {1 + \frac{x}{{{{a}_{1}}}}} \right)}^{{{{m}_{1}}}}}{{\left( {1 - \frac{x}{{{{a}_{2}}}}} \right)}^{{{{m}_{2}}}}}$
Мода $ - {{a}_{1}} < x < {{a}_{2}}$
II. $f = {{f}_{0}}{{\left( {1 - {{{{x}^{2}}} \mathord{\left/ {\vphantom {{{{x}^{2}}} {{{a}^{2}}}}} \right. \kern-0em} {{{a}^{2}}}}} \right)}^{m}}$ Мода (среднее) $ - a < x < a$
III. Гамма-распределение
$f = {{f}_{0}}{{e}^{{ - \gamma x}}}{{\left( {1 - {x \mathord{\left/ {\vphantom {x a}} \right. \kern-0em} a}} \right)}^{{\gamma a}}}$
Мода $ - a < x < \infty $
IV. Распределение Пирсона
четвертого типа
$f = {{f}_{0}}{{e}^{{{{ - {\gamma arctg}x} \mathord{\left/ {\vphantom {{ - {\gamma arctg}x} a}} \right. \kern-0em} a}}}}{{\left( {1 + \frac{{{{x}^{2}}}}{{{{a}^{2}}}}} \right)}^{{ - m}}}$
Среднее ${{ + \nu a} \mathord{\left/ {\vphantom {{ + \nu a} {(2m - 2)}}} \right. \kern-0em} {(2m - 2)}}$ $ - \infty < x < \infty $
V. $f = {{f}_{0}}{{e}^{{{{ - \gamma } \mathord{\left/ {\vphantom {{ - \gamma } x}} \right. \kern-0em} x}}}}{{x}^{{ - p}}}$ Начало кривой $0 < x < \infty $
VI. Распределение Пирсона
шестого типа $f = {{f}_{0}}{{(x - a)}^{{{{q}_{2}}}}}{{x}^{{ - {{q}_{1}}}}}$
Точка, отстоящая на $\left| a \right|$ от начала кривой $a < x < \infty $
VII. $f = {{f}_{0}}{{\left( {1 + {{{{x}^{2}}} \mathord{\left/ {\vphantom {{{{x}^{2}}} {{{a}^{2}}}}} \right. \kern-0em} {{{a}^{2}}}}} \right)}^{{ - m}}}$ Cреднее (мода) $ - \infty < x < \infty $

Из уравнения (56) следует, что мода равна $x = a.$ Из (59) следует, что для пирсоновской меры асимметрии (skewness) имеем

(61)
$Sk \equiv \frac{{ - a}}{{\sqrt {{{\mu }_{2}}} }} = \frac{{\sqrt {{{\beta }_{1}}} ({{\beta }_{2}} + 3)}}{{10{{\beta }_{2}} - 12{{\beta }_{1}} - 18}}.$

Далее, если $a = 0,$ то

${{{{{\text{d}}}^{2}}f} \mathord{\left/ {\vphantom {{{{{\text{d}}}^{2}}f} {{\text{d}}{{x}^{2}}}}} \right. \kern-0em} {{\text{d}}{{x}^{2}}}} = f\left( {{{b}_{0}} - {{b}_{2}}{{x}^{2}}} \right){{\left( {{{b}_{0}} + {{b}_{1}}x + {{b}_{2}}{{x}^{2}}} \right)}^{{ - 2}}},$
поэтому точки перегиба графика плотности распределения $f(x)$ определяются формулой ${{x}^{2}} = {{{{b}_{0}}} \mathord{\left/ {\vphantom {{{{b}_{0}}} {{{b}_{2}}}}} \right. \kern-0em} {{{b}_{2}}}}.$ Следовательно, у плотности распределений из семейства Пирсона существует не более чем две точки перегиба, и если их действительно две, то они отстоят от моды на одинаковом расстоянии. В общем случае может случиться, что одна из точек перегиба находится вне области, где сосредоточено распределение. Семейство кривых Пирсона, удовлетворяющее уравнению (56) приведено в табл. 2.

Дискриминант знаменателя в уравнении (56) равен: $D = b_{1}^{2}(1 - {{\kappa }^{{ - 1}}}),$ где

$\kappa = {{b_{1}^{2}} \mathord{\left/ {\vphantom {{b_{1}^{2}} {4{{b}_{0}}{{b}_{2}}}}} \right. \kern-0em} {4{{b}_{0}}{{b}_{2}}}} = \frac{{{{\beta }_{1}}{{{({{\beta }_{2}} + 3)}}^{2}}}}{{4(2{{\beta }_{2}} - 3{{\beta }_{1}} - 6)(4{{\beta }_{2}} - 3{{\beta }_{1}})}}$
– “каппа Пирсона” Общий интеграл уравнения (56) зависит от вида корней квадратного уравнения ${{b}_{0}} + {{b}_{1}}x$ + ${{b}_{2}}{{x}^{2}} = 0$ и определяется параметром $\kappa $-критерием Пирсона и дополнительными параметрами:

$\begin{gathered} a \equiv 2{{\beta }_{2}} - 3{{\beta }_{1}} - 6,\,\,\,\,c \equiv 8{{\beta }_{2}} - 15{{\beta }_{1}} - 36 = 0, \\ {\text{d}} \equiv {{\beta }_{2}} - {{\beta }_{1}} - 1 < 0. \\ \end{gathered} $

В табл. 3 приведены типы кривых Пирсона и соответствующие им критерии, а так же границы области кривых Пирсона. Граница 1 – это верхняя граница всех распределений, а граница 0 – граница кривых Пирсона.

Таблица 3.

Типы кривых Пирсона и соответствующие им критерии

Тип
кривой
Граница 0 I II III IV V VI VII Граница 1
Критерии     $\kappa < 0$ $\kappa = 0$ $\kappa = \pm \infty $ $0 < \kappa < 1$ $\kappa = 1$ $\kappa > 1$ $\kappa = 0$  
$d < 0$   ${{\beta }_{2}} < 3$ $a = 0$       ${{\beta }_{2}} > 3$ $c \geqslant 0$${{\mu }_{3}} = \infty $

Таким образом, если вдоль осей прямоугольной системы координат условиться откладывать отрезки, отвечающие величинам ${{\beta }_{1}}$ и ${{\beta }_{2}},$ то в плоскости ${{\beta }_{2}}0\,{{\beta }_{1}}$ различным типам кривых Пирсона будут соответствовать области, кривые и точки.

На рис. 3 указано такое разбиение плоскости ${{\beta }_{2}}0\,{{\beta }_{1}}$ для основных типов кривых Пирсона I–VII. Прямая линия с уравнением $d \equiv {{\beta }_{2}} - {{\beta }_{1}} - 1 = 0$ представляет собой верхнюю границу для допустимых точек $({{\beta }_{2}},{{\beta }_{1}}),$ так как не существует распределений, для которых $d < 0.$ Кроме этого, если кривая принадлежит семейству Пирсона, причем $c \equiv 8{{\beta }_{2}} - 15{{\beta }_{1}} - 36 \geqslant 0,$ то ${{\mu }_{3}} = \infty .$ На рис. 3 прямая с уравнением $c = 0$ служит нижней границей точек с координатами $({{\beta }_{2}},{{\beta }_{1}}).$

Заметим, что наиболее типичная колокообразная форма кривой типа I наблюдается тогда, когда ${{m}_{1}}$ и ${{m}_{2}}$ положительны (см. область I на рис. 3): для J-образного бета-распределения (см. область I(J)) один из этих показателей отрицателен. Если же ${{m}_{1}}$ и ${{m}_{2}}$ оба отрицательны, то бета-распределение имеет U-образную форму (см. на рис. 3 область I(U)).

– Граница области I(J) задается уравнением

$\begin{gathered} 4(4{{\beta }_{2}} - 3{{\beta }_{1}}){{(5{{\beta }_{2}} - 6{{\beta }_{1}} - 9)}^{2}} = \\ = {{\beta }_{1}}({{\beta }_{2}} + {{3}^{2}})(8{{\beta }_{2}} - 9{{\beta }_{1}} - 12). \\ \end{gathered} $

– Линия III типа: $2{{\beta }_{2}} - 3{{\beta }_{1}} - 6 = 0.$

– Линия V типа: ${{\beta }_{1}}{{({{\beta }_{2}} + 3)}^{2}}$ = $4(4{{\beta }_{2}} - 3{{\beta }_{1}})$ × $ \times \,\,(2{{\beta }_{2}} - 3{{\beta }_{1}} - 6).$

– Линия, ниже которой (а также на ней самой) для всех кривых Пирсона:

${{\mu }_{8}} = \infty ,\,\,\,\,8{{\beta }_{2}} - 15{{\beta }_{1}} - 36 = 0.$

Было установлено (Pearson, Hartley, 1954), что пирсоновские распределения нередко хорошо соответствуют результатам наблюдений. Другое достоинство этих распределений (в частности, распределений I и III типов) состоит в том, что ими можно с хорошей точностью приближать теоретические распределения, зная их моменты. Систематическое изложение техники подбора аппроксимирующих распределений было дано Elderton (1953). Все распределения Пирсона определяются своими четырьмя первыми моментами, за исключением распределений, задаваемых меньшим количеством моментов.

Пирсоновский метод подгонки состоит в следующем:

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

2. Вычисляются значения ${{\beta }_{1}}$ и ${{\beta }_{2}}$ и величина $\kappa $ (см. (61)), и, следовательно, определяется тип распределений.

3. Эмпирические моменты приравниваются моментам подходящего распределения, которые выражены в терминах его параметров.

4. Полученные уравнения разрешаются относительно неизвестных параметров и, следовательно, находится искомое распределение (см. Большев, Смирнов, 1983).

В следующем разделе рассмотрим пример иллюстрирующий этот процесс.

ВЫЧИСЛЕНИЕ МОМЕНТОВ РЕШЕНИЯ КИНЕТИЧЕСКОГО УРАВНЕНИЯ

Схему вычисления моментов решения рассмотрим на примере кинетического уравнения (1) с ядром коагуляции

(62)
$K({v},w) = G{{{v}}^{{\alpha }}}{{w}^{{\alpha }}}.$

Подставляя (62) в (21) и учитывая (44), получим систему уравнений для определения первых пяти моментов решения кинетического уравнения

(63)
$\begin{gathered} \frac{{\text{d}}}{{{\text{d}}\tau }}m_{0}^{'} = {{B}_{0}}{{\left( {m_{{\alpha }}^{'}} \right)}^{2}},\,\,\,\,\frac{{\text{d}}}{{{\text{d}}\tau }}m_{1}^{'} = 0, \\ \frac{{\text{d}}}{{{\text{d}}\tau }}m_{2}^{'} = {{B}_{2}}{{\left( {m_{{1 + {\alpha }}}^{'}} \right)}^{2}},\,\,\,\,\frac{{\text{d}}}{{{\text{d}}\tau }}m_{3}^{'} = {{B}_{3}}m_{{1 + {\alpha }}}^{'}m_{{2 + {\alpha }}}^{'}, \\ \frac{{\text{d}}}{{{\text{d}}\tau }}m_{4}^{'} = {{B}_{4}}m_{{1 + {\alpha }}}^{'}m_{{3 + {\alpha }}}^{'} + {{B}_{5}}{{\left( {m_{{2 + {\alpha }}}^{'}} \right)}^{2}}, \\ \end{gathered} $
где дифференцирование проводится по $\tau = Gm_{1}^{{2\alpha }}t.$ Здесь

(64)
$\begin{gathered} {{B}_{0}} = - \frac{1}{2}m_{1}^{{ - 2{\alpha }}}m_{0}^{{ - 1}}(0)m_{{1 + {\alpha }}}^{2}(0), \\ {{B}_{2}} = m_{1}^{{ - 2{\alpha }}}m_{2}^{{ - 1}}(0)m_{{1 + {\alpha }}}^{2}(0), \\ {{B}_{3}} = 3m_{1}^{{ - 2{\alpha }}}{{m}_{{1 + {\alpha }}}}(0){{m}_{{2 + \alpha }}}(0){{m}_{3}}(0), \\ {{B}_{4}} = 4m_{1}^{{ - 2{\alpha }}}{{m}_{{1 + {\alpha }}}}(0){{m}_{{3 + {\alpha }}}}(0)m_{4}^{{ - 1}}(0), \\ {{B}_{5}} = 3m_{1}^{{ - 2{\alpha }}}m_{{2 + {\alpha }}}^{2}(0)m_{4}^{{ - 1}}(0). \\ \end{gathered} $

Для случая двухточечной интерполяции дробных моментов через целые доопределяющие уравнения для системы (63) можно записать в виде:

(65)
$\begin{gathered} m_{{\alpha }}^{'} = {{\left( {m_{0}^{'}} \right)}^{{1 - {\alpha }}}},\,\,\,\,m_{{1 + {\alpha }}}^{'} = {{\left( {m_{2}^{'}} \right)}^{{\alpha }}}, \\ m_{{2 + {\alpha }}}^{'} = {{\left( {m_{2}^{'}} \right)}^{{1 - {\alpha }}}}{{\left( {m_{3}^{'}} \right)}^{{\alpha }}},\,\,\,\,m_{{3{\alpha }}}^{'} = {{\left( {m_{{3 + {\alpha }}}^{'}} \right)}^{{1 - {\alpha }}}}{{\left( {m_{4}^{'}} \right)}^{{\alpha }}}. \\ \end{gathered} $

Подставляя (65) в (63), получим

(66)

Первые четыре уравнения имеют аналитические решения, которые можно записать в виде:

при $\alpha = 1{\text{/}}2$

(67.1)
$\begin{gathered} m_{0}^{'} = {\text{exp}}({{B}_{0}}\tau ),\,\,\,\,m_{1}^{'} = 1,\,\,\,\,m_{2}^{'} = {\text{exp}}({{B}_{2}}\tau ), \\ m_{3}^{'} = {{\left[ {1 + ({{{{B}_{3}}} \mathord{\left/ {\vphantom {{{{B}_{3}}} {2{{B}_{2}}}}} \right. \kern-0em} {2{{B}_{2}}}}){\text{exp}}({{B}_{2}}\tau )} \right]}^{2}}; \\ \end{gathered} $
при $\alpha \ne 1{\text{/}}2$

(67.2)
$\begin{gathered} m_{0}^{'} = {{\left[ {1 + (2\alpha - 1){{B}_{0}}\tau } \right]}^{{{{ - 1} \mathord{\left/ {\vphantom {{ - 1} {(1 - 2{\alpha })}}} \right. \kern-0em} {(1 - 2{\alpha })}}}}},\,\,\,\,m_{1}^{'} = 1, \\ m_{2}^{'} = {{\left[ {1 + (2\alpha - 1){{B}_{2}}\tau } \right]}^{{{1 \mathord{\left/ {\vphantom {1 {(1 - 2{\alpha })}}} \right. \kern-0em} {(1 - 2{\alpha })}}}}},\,\,\,\,m_{3}^{'} = \\ = {{\left\{ {1 + \frac{{{{B}_{3}}}}{{2{{B}_{2}}}}\left[ {{{{\left( {1 + \frac{{{{B}_{2}}\tau }}{{1 - 2\alpha }}} \right)}}^{{{{2(1 - {\alpha })} \mathord{\left/ {\vphantom {{2(1 - {\alpha })} {(1 - 2{\alpha })}}} \right. \kern-0em} {(1 - 2{\alpha })}}}}} - 1} \right]} \right\}}^{{{1 \mathord{\left/ {\vphantom {1 {(1 - {\alpha })}}} \right. \kern-0em} {(1 - {\alpha })}}}}}. \\ \end{gathered} $

По первым пяти моментам решения можно сделать параметрическую оценку самого решения на основе диаграммы Пирсона, построенных в координатах квадрата асимметрии ${{\beta }_{1}}$ и эксцесса ${{\beta }_{2}}$оцениваемого распределения

(68)
$\begin{gathered} {{\beta }_{1}} = {{\left[ {{{\mu _{3}^{'}} \mathord{\left/ {\vphantom {{\mu _{3}^{'}} {{{{(\mu _{2}^{'})}}^{{{\text{3/2}}}}}}}} \right. \kern-0em} {{{{(\mu _{2}^{'})}}^{{{\text{3/2}}}}}}}} \right]}^{2}},\,\,\,\,{{\beta }_{2}} = {{\mu _{4}^{'}} \mathord{\left/ {\vphantom {{\mu _{4}^{'}} {{{{\left( {\mu _{2}^{'}} \right)}}^{2}}}}} \right. \kern-0em} {{{{\left( {\mu _{2}^{'}} \right)}}^{2}}}}, \\ \mu _{i}^{'} = \int\limits_0^\infty {{{{({v} - {\bar {v}})}}^{i}}p({v}){\text{d}}{v}} . \\ \end{gathered} $

Чтобы решить, какое из семейств распределений Пирсона аппроксимирует искомое распределение, нужно использовать годограф совокупности параметров ${{\beta }_{1}}(\tau )$ и ${{\beta }_{2}}(\tau ).$

На рис. 4 представлена диаграмма Пирсона для классических распределений, на которой нанесены годографы точек с координатами ${{\beta }_{1}}(\tau )$ и ${{\beta }_{2}}(\tau ),$ рассчитанные для ядра вида (62) при $\alpha = 0,$ $\alpha = 1{\text{/}}4$ и $\alpha = 1{\text{/}}2$ (Логунов, 1979). Годограф 4 соответствует уравнению с ядром $K({v},w)$ = $G({v} + w).$ При этом в качестве начальной плотности распределения частиц по размерам было выбрано гамма-распределение

(69)
$\begin{gathered} p({v},0) = \frac{{{{\lambda }^{{\beta }}}}}{{\Gamma (\beta )}}{{{v}}^{{{\beta } - 1}}}{\text{exp}}( - \lambda {v}), \\ {\text{при}}\,\,\,\,\lambda = 4 \times {{10}^{9}},\,\,\,\,\beta = 2, \\ \end{gathered} $
причем начальные значения моментов определялись как

(70)
${{m}_{{\gamma }}}(0) = {{\Gamma (\beta + \gamma )} \mathord{\left/ {\vphantom {{\Gamma (\beta + \gamma )} {{{\lambda }^{{\gamma }}}\Gamma (\beta )}}} \right. \kern-0em} {{{\lambda }^{{\gamma }}}\Gamma (\beta )}}.$
Рис. 4.

Диаграмма Пирсона с годографами решений кинетического уравнения соответствующих ядрам (62) при 1$\alpha = 0;$ 2$\alpha = 1/4;$ 3$\alpha = 1/2;$ 4$K({v},w)$ = $G({v} + w).$

В рассматриваемом случае точки ${{\beta }_{1}},$ ${{\beta }_{2}}$ лежат в области, ограниченной графиками кривых гамма-распределения и логнормального распределения (см. рис. 2). Следовательно, для описания полученных экспериментальных данных нужно выбрать семейство III (см. рис. 3) распределений Пирсона.

ЗАКЛЮЧЕНИЕ

Представленная работа является необходимым звеном в ряду исследований, посвященных математическим аспектам моделирования допланетного газопылевого облака, в процессе эволюции которого формируются и взаимодействуют друг с другом разномасштабные пылевые фрактальные агломераты, служащие, в конечном счете, основой зародышей рыхлых протопланетезималей. Поскольку строгое решение задачи образования и эволюции фрактальных кластеров включает одновременно с оценкой их скоростей также и определение функции распределения (спектра) кластеров по размерам (т.е. решение обобщенного нелинейного пространственно-неоднородного кинетического уравнения Смолуховского), то в ряде работ (Колесниченко, 2000; Колесниченко, Маров, 2014a; 2014б; Kolesnichenko, 2001; Marov, Kolesnichenko, 2001) были получены практически важные модели для ядер коагуляции фрактальных кластеров. При этом были рассмотрены две группы моделей образования кластеров в дисковой фрактальной среде. К первой группе моделей ядер коагуляции относятся модели, обусловленные прилипанием мономеров к кластеру; при этом учитывался связанный с характером движения отдельного мономера переход от кинетического режима к диффузионному. Ко второй группе моделей относятся модели, описывающие рост фрактальных агрегатов в результате ассоциации двух кластеров. Полученные алгебраические выражения для коэффициентов коагуляции во фрактальных средах существенно усложняются по сравнению с их видом в континуальной коагулирующей среде. Это обстоятельство дополнительнозатрудняет решение кинетического уравнения Смолуховского, моделирующего процессы аккумуляции допланетных тел. По этой причине требуется разработка подхода, позволяющего находить полуэмпирические распределения допланетных тел по размерам, когда имеются ограниченные данные о коэффициентах коагуляциив стохастической фрактальной системе.

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

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

Список литературы

  1. Большев Л.Н., Смирнов Н.В. Таблицы математической статистики. М.: Наука. Главная редакция физико-математической литературы, 1983. 416 с.

  2. Волощук В.М., Седунов Ю. Процессы коагуляции в дисперсных средах. М.: Гидрометиздат, 1975. 320 с.

  3. Волощук В.М. Кинетическая теория коагуляции. М.: Гидрометиздат, 1984. 283 с.

  4. Демидович Б.П., Марон И.А. Основы вычислительной математики. М.: Физматгиз, 1963. 664 с.

  5. Кендал М.Дж., Стьюарт А. Теория распределений. Пер. с англ. / Ред. Колмогоров А.Н. М.: Наука, 1966. 588 с.

  6. Колесниченко А.В. Проблемы моделирования процессов массопереноса и коагуляции в допланетном газопылевом облаке // Препр. ИПМ им. М.В. Келдыша РАН. 2000. № 41. 30 с.

  7. Колесниченко А.В. О синергетическом механизме возникновения когерентных структур в континуальной теории развитой турбулентности // Астрон. вестн. 2004. Т. 38. № 5. С. 405–427. (Kolesnichenko A. Synergetic mechanism of the development of coherent structures in the continual theory of developed turbulence // Sol. Syst. Res. 2004. V. 38. № 5. P. 351–371).

  8. Колесниченко А.В. О роли индуцированных шумом неравновесных фазовых переходов в структурировании гидродинамической турбулентности // Астрон. вестн. 2005. Т. 39. № 3. С. 243–262. (Kolesnichenko A. The role of noise-induced nonequilibrium phase transitions in the structurization of hydrodynamic turbulence // Sol. Syst. Res. 2005. V. 39. № 3. P. 214–230.)

  9. Колесниченко А.В., Маров М.Я. Моделирование процессов образования пылевых фрактальных кластеров как основы рыхлых прото-планетезималей в Солнечном допланетном облаке // Препр. ИПМ им. М.В. Келдыша РАН. 2014a. № 75. 44 с.

  10. Колесниченко А.В., Маров М.Я. Турбулентность и самоорганизация. Проблемы моделирования космических и природных сред. М.: БИНОМ. Лаборатория знаний, 2014б. 632 с.

  11. Крамер Г. Математические методы статистики. М.: ИЛ, 1948. 648 с.

  12. Левин Б.Р. Теоретические основы статистической радиотехники. М.: Советское радио, 1974. 656 с.

  13. Логунов В.И. Обезвоживание и обессоливание нефтей. М.: Химия, 1979. 216 с.

  14. Попель А.С., Регирер С.А., Шадрина Н.Х. Об уравнениях кинетики агрегационных процессов в суспензиях // Прикл. математика и механика. 1975. Т. 39. № 1. С. 131–143.

  15. Сафронов В.С. О гравитационной неустойчивости в плоских вращающихся системах с осевой симметрией // ДАН СССР. 1960. Т. 130. № 1. С. 53–56.

  16. Сафронов В.С. Эволюция допланетного облака и образование Земли и планет. М.: Наука, 1969. 244 с.

  17. Синайский Э.Г. Гидродинамика физико-химических процессов. М.: Недра, 1997. 339 с.

  18. Barrett J.C., Jheeta J.S. Improving the accuracy of the moments method for solving the Aerosol General Dynamic Equation // J. Aerosol Sci. 1996. V. 27. P. 1135–1142.

  19. Barrett J.C., Webb N.A. A comparison of some approximate methods for solving the aerosol general dynamic equation // J. Aerosol Sci. 1998. V. 29. P. 31–39.

  20. Blum J. Grain growth and coagulation // ASP Conf. Ser. V. 309. Astrophysics of Dust / Ed. Witt A.N., Clayton G.C., Draine B.T. San Francisco: ASP, 2004. P. 369–391.

  21. Chan T.L., Lin J.Z., Zhou K., Chan C.K. Simultaneous numerical simulation of nano and fine particle coagulation and dispersion in a round jet // J. Aerosol Sci. 2006. V. 37. P. 1545–1561.

  22. Dubrulle B., Morfill G., Sterzik M. The dust subdisk in the protoplanetary nebula// Icarus. 1995. V. 114. P. 237–246.

  23. Elderton W.P. Frequen curves and correlation. Cambridg Univ. Press, 1953. 290 p.

  24. Frenklach M., Harris S.J. Aerosol dynamics modeling using the method of moments // J. Coll. Interface Sci. 1987. V. 118. P. 252–261.

  25. Frenklach M. Method of moments with interpolative closure // Chem. Eng. Sci. 2002. V. 57. P. 2229–2239.

  26. Goldreich P., Ward W.R. The formation of planetesimals // Astrophys. J. 1973. V. 183. № 3. P. 1051–1061.

  27. Hahn G.J., Shapiro S.S. Statistical models in engineering. N.Y.: Wiley, 1967. 376 p.

  28. Kolesnichenko A.V. Hydrodynamic aspects of modeling of the mass transfer and coagulation processes in turbulent accretion disks // Sol. Syst. Res. 2001. V. 35. № 2. P. 125–140.

  29. Kolesnichenko A.V., Marov M.Ya. Fundamentals of the mechanics of heterogeneous media in the circumsolar protoplanetary cloud: The effects of solid particles on disk turbulence // Sol. Syst. Res. 2006. V. 40. № 1. P. 1–56.

  30. Lee K.W. Change of particle size distribution during brownian coagulation // J. Colloid Interface Sci. 1983. V. 92. № 2. P. 315–325.

  31. Marchisio D.L., Fox R.O. Solution of population balance equations using the direct quadrature method of moments // J. Aerosol Sci. 2005. V. 36. P. 43–73.

  32. Marov M.Ya., Kolesnichenko A.V. Mechanics of turbulence of multicomponent gases. Kluwer Acad. Publ, 2001. 375 p.

  33. Marov M.Ya., Kolesnichenko A.V. Turbulence and Self-Organization.Modeling Astrophysical Objects. Springer, 2013. 657 p.

  34. McGraw R. Description of aerosol dynamics by the quadrature method of moments // Aerosol Sci. Technol. 1997. V. 27. P. 255–265.

  35. Mingzhou Yu., Jianzhong Lin., Tatleung C.A. New moment method for solving the coagulation equation for particles in brownian motion // Aerosol Sci. and Technology. 2008. V. 42. P. 705–713.

  36. Nakagawa Y., Nakazawa K., Hayashi C. Growth and sedimentation of dust grains in the primordial solar nebula // Icarus. 1981. V. 45. P. 517–528.

  37. Nakagawa Y., Hayashi C., Nakazawa K. Accumulation of planetesimals in the solar nebula // Icarus. 1983. V. 54. P. 361–376.

  38. Nakagawa Y., Sekiya M., Hayashi C. Settling and growth of dust particles in a laminar phase of a low-mass solar nebula // Icarus. 1986. V. 67. P. 375–390.

  39. Nakamoto T., Nakagawa Y. Formation, early evolution, and gravitational stability of protoplanetary disks // Astrophys. J. 1994. V. 421. P. 640–651.

  40. Okuzumi S., Tanaka H., Takeuchu T., Sakagami M.-A. Electrostatic barrier against dust growth in protoplanetary disks. 1. Classifying the evolution of size distribution // Astrophys. J. 2011. V. 731. P. 95–115.

  41. Ormel C.W., Spaans M., Tielens A.G.G.M. Dust coagulation in protoplanetary disks: Porosity matters// Astron. and Astrophys. 2007. V. 461. P. 215–236.

  42. Pearson E.S., Hartley H.O. Biometrika Tables for Statisticians. V. I. Cambridge: Cambridge Univ. Press, 1954. 25 p.

  43. Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P. Numerical Recipes in Fortran: 2 nd edition. Cambridge Univ. Press, 1992. 994 p.

  44. Smoluchowski M.V. Experiments on a mathematical theory of kinetic coagulation of colloid solutions // Z. Phys. Chem.-Stochiometrie und Verwandtschaftslehre 1917. V. 92. P. 129–68.

  45. Suyama T., Wada K., Tanaka H. Numerical simulation of density evolution of dust aggregates in protoplanetary disks. I. Head-on collisions // Astrophys. J. 2008. V. 684. P. 1310–1322.

  46. Suyama T., Wada K., Tanaka H., Okuzumi S. Geometrical cross sections of dust aggregates and a compression model for aggregate collisions // arxiv:1205.1894vl [astro-ph. EP]. 2012. 28 p.

  47. Tarasov V.E. Fractional hydrodynamic equations for fractal media // Annals of Physics. 2005. V. 318. № 2. P. 286–307.

  48. Tarasov V.E. Fractional dynamics: Applicationsof fractional calculus to dynamics of particles, fields and media. Springer: Higher Education Press, 2010. 516 p.

  49. Wada K., Tanaka H., Suyama T., Kimura H., Yamamoto T. Simulation of dust aggregate collisions. II. Compression and disruption of three-dimensional aggregates in head-on collisions // Astrophys. J. 2008. V. 677. P. 1296–1308.

  50. Weidenschilling S.J. Dust to planetesimals: Settling and coagulation in the solar nebula// Icarus. 1980. V. 44. P. 172–189.

  51. Youdin A.N., Shu F. Planetesimal formation by gravitational instability // Astrophys. J. 2002. V. 580. P. 494–505.

  52. Yu M.Z., Lin J.Z., Chan T.L. A new moment method for solving the coagulation equation for particles in Brownian motion // Aerosol Sci. Technol. 2008. V. 42. P. 705–713.

  53. Yu M.Z., Lin J., Jin H., Jiang Y. The verification of the Taylor-expansion moment method for the nanoparticle coagulation in the entire size regime due to Brownian motion // J. Nanopart. Res. 2011.V. 13. P. 2007–2020.

Дополнительные материалы отсутствуют.