Журнал вычислительной математики и математической физики, 2020, T. 60, № 10, стр. 1787-1794

ВВЕДЕНИЕ

Работа посвящена математическому моделированию процесса образования и эволюции движения космических облаков – скоплений большого числа космических частиц – на базе уравнений движения классической задачи $n$-тел (см., например, [1]–[5]). В литературе имеются и более сложные модели, основанные на гидродинамических, физико-химических, релятивистских и других аспектах космогонических процессов (см., например, [6]–[8]). В Cолнечной системе такие скопления обнаружены в устойчивых треугольных точках либрации систем Солнце–Юпитер, Солнце–Сатурн, Солнце–Земля и др. (см. [9]).

Существование скопления частиц в окрестности треугольных точек либрации системы Земля–Луна впервые было зафиксировано в 1961 г. (см. [10]) польским астрономом К. Кордылевским. Затем последовала череда противоречащих одна другой публикаций, подтверждающих и ставящих под сомнение существование этих облаков, их называли “неуловимые облака Кордылевского”. Дискуссии длились более полувека. Объяснение этого феномена было дано в работах [11]–[20] с учетом гравитационного и светового возмущения от Солнца и электростатических взаимодействий частиц.

В отличие от указанных систем, в формировании планетных систем основную роль играют начальный кинетический момент, внутренние взаимодействия частиц и диссипация механической энергии, которая происходит, главным образом, за счет столкновения частиц (см. [21]–[23]).

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

Другая трудность интегрирования задачи $n$-тел связана с квадратичным ростом числа парных взаимодействий $n(n - 1){\text{/}}2$ частиц с ростом $n$. Эта трудность преодолевается введением упрощенных алгоритмов вычисления: древовидного упорядочения взаимодействий с отбрасыванием несущественных дальнодействий и близкодействий (см. [24], [25]), метод частиц в ячейке (particles in cell – PiC) (см. [5]) и др. Эти алгоритмы реализованы для последовательного и параллельного вычислений и для применения графических процессоров (см. [26]–[28]).

В настоящей работе исследование проводится в рамках классической задачи $n$-тел Ньютона и гипотезы Ньютона о неупругом ударе твердых тел. За основу взята модель, разрабатывавшаяся с начала 80-х годов прошлого века Т.М. Энеевым и Н.Н. Козловым (см. [3]). В этой работе обосновывались дополнения к гипотезе Шмидта, связанные с механизмом взаимодействия сталкивающихся частиц. В 2001 г. Т.М. Энеев (см. [4]) выдвинул гипотезу о том, что из отделившегося пояса частиц могут образоваться только маленькие спутники, а происхождение больших спутников объясняется захватом других небесных тел с близкими орбитами. В [29] с применением модифицированного потенциала была показана возможность образования на ранней стадии не только мелких, но и крупных спутников планеты, а также систем типа двойных звезд. В настоящей работе основное внимание уделено алгоритмической реализации модели, основанной на теории неупругого удара твердых тел.

1. ОПИСАНИЕ МОДЕЛИ

В статье моделирование проводится в рамках классической задачи $n$-тел, притягивающихся по закону всемирного тяготения Ньютона:

где ${{m}_{i}},\;{{r}_{i}},\;{{v}_{i}}$ – массы, радиус-векторы и скорости центров масс частиц в инерциальной системе отсчета Кёнига с началом в центре масс системы соответственно. В случае контактного взаимодействия используется теория удара твердых тел, базирующаяся на введенном Ньютоном понятии коэффициента восстановления относительной скорости частиц после удара $(0 \leqslant e \leqslant 1)$.

При рассмотрении ударов поверхности частиц предполагаются гладкими сферами с диаметрами ${{d}_{i}}$ и с центрами в центрах масс частиц. Это соответствует предположению о малости кинетической энергии вращения частицы вокруг центра масс по сравнению с кинетической энергией движения центра масс. Сила трения в плоскости, касательной к контактирующим поверхностям, не учитывается. Считаем, что при ударе скорости меняются мгновенно, сохраняется количество движения соударяющихся частиц, и, согласно гипотезе Ньютона, нормальная составляющая ${{u}_{{ij}}}$ разности скоростей ${{v}_{{ij}}} = {{v}_{j}} - {{v}_{i}}$ частиц $i$ и $j$ после удара меняет знак и умножается на коэффициент восстановления $0 \leqslant e \leqslant 1$ (штрихами отмечены скорости после удара):

При $e = 0$ имеем абсолютно неупругий удар, при $e = 1$ – абсолютно упругий удар. Из (2) следуют формулы для скоростей после удара и для величины потери кинетической энергии при ударе В частном случае (${{m}_{i}} = {{m}_{j}} = m$) имеем

Такой подход сопряжен со следующими трудностями. Во-первых, с увеличением числа частиц n объем вычислений парных силовых взаимодействий возрастает пропорционально $n(n - 1){\text{/}}2$. Эта трудность преодолевается применением PiC-метода (Particles in Cell) (см. [17]) для вычисления ньютоновского потенциала как решения уравнения Пуассона. Во-вторых, в процессе уплотнения системы частота столкновений увеличивается, и при достижении некоторого предела необходимо переходить на другую модель взаимодействия, например, с модифицированным потенциалом и рассмотрением взаимного проникновения частиц. В этом случае PiC-метод неприменим, и, кроме того, при сближении хотя бы одной пары частиц приходится уменьшать шаг интегрирования. Указанные особенности необходимо отслеживать на каждом шаге интегрирования.

2. АЛГОРИТМ ИНТЕГРИРОВАНИЯ

В силу необходимости сокращения шага интегрирования уравнений (1) при сближении частиц с нарушением условия

и необходимости пересчета начальных условий при каждом ударе, интегрирование проводится одношаговым методом Кутты–Мерсона с автоматическим определением шага.

На шаге интегрирования $[{{t}_{k}},{{t}_{{k + 1}}}]$, на котором нарушается условие (6) для пары $ij$, уточняется момент времени первого касания сферических поверхностей, связанных с частицами. Для этого строится кубическая интерполяция функции $t(\rho )$, $\rho = {\text{|}}{{r}_{{ij}}}{\kern 1pt} {\text{|}}$, по известным значениям этой функции и ее производной $t{\kern 1pt} ' = dt{\text{/}}d\rho = {\text{|}}{{r}_{{ij}}}{\kern 1pt} {\text{|/}}({{r}_{{ij}}}{{{v}}_{{ij}}})$ на концах отрезка $[{{\rho }_{k}},{{\rho }_{{k + 1}}}]$, соответствующего отрезку $[{{t}_{k}},{{t}_{{k + 1}}}]$:

где При $\rho = {{d}_{{ij}}}$ получаем момент времени ${{t}_{{ij}}}$ касания частиц $i$ и $j$.

Алгоритм интегрирования включает следующие действия.

1. На каждом шаге интегрирования определяется множество $\alpha $ пар индексов $ij$, для которых нарушается условие (6).

2. Для каждой пары $ij$ из $\alpha $ по значениям ${{t}_{k}},\;{{r}_{{ij}}}({{t}_{k}})$, ${{v}_{{ij}}}({{t}_{k}})$ и ${{t}_{{k + 1}}},\;{{r}_{{ij}}}({{t}_{{k + 1}}})$, ${{{v}}_{{ij}}}({{t}_{{k + 1}}})$ вычисляем значения $dt{\text{/}}d\rho $ в точках ${{t}_{k}}$ и ${{t}_{{k + 1}}}$, строим кубическую интерполяцию функции $t(\rho )$ на отрезке $[{{t}_{k}},{{t}_{{k + 1}}}]$ и вычисляем значение ${{t}_{{ij}}} = t({{d}_{{ij}}})$.

3. Среди всех значений ${{t}_{{ij}}}$ выбираем наименьшее $t{\kern 1pt} * = min{{t}_{{ij}}}$ и определяем подмножество $\beta \subset \alpha $, для которого этот минимум достигается с заданной точностью.

4. Делаем шаг интегрирования от точки ${{t}_{k}}$ до точки $t{\kern 1pt} *$, при этом для всех контактирующих пар из множества $\beta $ будут выполняться с заданной точностью равенства ${\text{|}}{{r}_{{ij}}}{\kern 1pt} {\text{|}} = {{d}_{{ij}}}$.

5. Если множество $\beta $ состоит из одного элемента или если элементов несколько, но ни одна частица не входит более, чем в одну пару из $\beta $, то это – независимые парные столкновения и проводится независимый пересчет скоростей по формулам (3) для каждой пары из $\beta $.

6. Если в множестве $\beta $ имеется одно или несколько подмножеств $\sigma $, для которых в каждой паре хотя бы одна из частиц входит хотя бы в одну другую пару частиц из $\sigma $, то мы имеем кратный удар, пересчет скоростей в котором будет рассмотрен ниже.

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

В дальнейшем будем различать напряженный и ненапряженный контакт частиц. Будем говорить, что пара соприкасающихся частиц $ij$ находится в напряженном контакте, если, кроме условия ${\text{|}}{{r}_{{ij}}}{\kern 1pt} {\text{|}} = {{d}_{{ij}}}$, еще выполнено условие ${{r}_{{ij}}}{{v}_{{ij}}} > 0$, означающее, что при продолжении непрерывного движения условие ${\text{|}}{{r}_{{ij}}}{\kern 1pt} {\text{|}} \geqslant {{d}_{{ij}}}$ нарушится.

3. АЛГОРИТМ РАЗРЕШЕНИЯ КРАТНОГО УДАРА

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

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

Алгоритм пересчета скоростей состоит из следующих действий.

1. Упорядочиваем последовательность пар в множестве $\sigma $ по возрастанию времени наступления контакт и в заданном порядке производим пересчет скоростей пар, находящихся в напряженном контакте (с учетом ранее произведенных пересчетов).

2. Если при переборе всех пар из $\sigma $ была обнаружена хотя бы одна напряженная пара, то перебор повторяется.

3. Если при переборе всех пар из $\sigma $ не было ни одной пары в напряженном состоянии, пересчет заканчиваем, и возобновляется интегрирование с новыми начальными условиями.

4. ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ

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

Для наглядной и простой иллюстрации приведем результаты двумерного моделирования 1024 одинаковых шариков с массой $m$ и диаметром $d$. В начальный момент времени шарики равномерно с шагом $10d$ заполняют два прямоугольника, содержащие $32 \times 19$ и $32 \times 13$ шариков. Величина постоянной всемирного тяготения принята равной $\gamma = 6.6743 \times {{10}^{{ - 11}}}$ м³ с^–2 кг^–1. В качестве единицы длины принят диаметр шариков $L = d$, единицы массы – масса шариков $M = m$, а единица времени $T$ выбрана так, чтобы в новых единицах выполнялось равенство $\gamma = 1{{L}^{3}}{{T}^{{ - 2}}}{{M}^{{ - 1}}}$. Если принять плотность материала шариков равной $\rho = 2208$ кг м^–3, примерно соответствующую плотности хондритов, из которых состоят метеорные тела солнечной системы, то единица времени составит T = 3600 c = 1 ч. При таком выборе единиц измерения уравнения движения и результаты расчета не меняются при изменении размера частиц. На фиг. 1 показано начальное распределение шариков. Большая часть шариков не видна за левой и правой границами рисунка. Скорости прямоугольников равны $1.54$ и $ - 1.2$ (на фигурах скорости показаны короткими штрихами), расстояние между передними кромками прямоугольников равно $90$. Коэффициент восстановления равен $e = 0.97$.

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

5. УПРОЩЕНИЕ АЛГОРИТМА

При интегрировании уравнений (1) существенную долю машинного времени занимает вычисление гравитационного потенциала и гравитационных сил, действующих на каждую частицу со стороны остальных частиц. При значительном количестве частиц и, особенно, при рассмотрении пространственных движений используется ускоренный расчет потенциала с помощью PIC-метода (Particles In Cell) (см. [5]), в котором система частиц рассматривается как гравитирующая сплошная среда, погруженная в прямоугольную неподвижную сетку. Потенциал Ньютона $U$ распределенной гравитирующей массы обладает замечательным свойством: он удовлетворяет уравнению Пуассона $\Delta U(r) = \rho $ с массовой плотностью $\rho $, пропорциональной количеству частиц, попавших в ячейку. Значения функции $U$ и плотности $\rho $ вычисляются в центрах ячеек сетки. Уравнение Пуассона интегрируется методом быстрого преобразования Фурье FFT (Fast Fourier Transform).

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

При ${\text{|}}{{r}_{j}} - r{\kern 1pt} {\text{|}} = d$ функции ${{V}_{j}}$ имеют минимум, и определяемые ими силы обращаются в нуль. При ${\text{|}}{{r}_{j}} - r{\kern 1pt} {\text{|}} > d$ эти силы близки к ньютоновским силам тяготения, а при ${\text{|}}{{r}_{j}} - r{\kern 1pt} {\text{|}} < d$ имеем отталкивающие силы, стремящиеся к бесконечности при ${\text{|}}{{r}_{j}} - r{\kern 1pt} {\text{|}} \to 0$. При таком подходе возможно сквозное интегрирование уравнений движения без выделения моментов соударения частиц, и устраняется неоднозначность решения, однако возникают трудности, связанные с жесткостью системы.

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