Цель работы

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

Задание

В рамках этапа «Алгоритмы» требуется:

• Описать общий подход к моделированию агрегации, ограниченной диффузией (DLA).

• Привести математические формулы для вычисления фрактальной размерности двумя методами.

• Составить пошаговые алгоритмы для базовой модели DLA.

• Описать модификации алгоритма для случаев с переменной вероятностью прилипания.

• Изложить алгоритмы для бессеточной модели, баллистической агрегации и кластер-кластерной агрегации.

Общий подход к моделированию

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

Процесс моделируется пошагово:

• Задается начальное состояние системы (затравочная частица или подложка).

• Генерируется новая блуждающая частица на границе области.

• Частица совершает случайные блуждания.

• При контакте с кластером частица присоединяется к нему с заданной вероятностью.

• Процесс повторяется до достижения нужного размера кластера.

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

Математический аппарат

Случайные блуждания

Движение частицы описывается дискретным случайным процессом. На каждом шаге $t$ координаты частицы изменяются на вектор $(dx, dy)$, где компоненты выбираются случайно из множества разрешенных направлений.

Для сеточной модели:

Направления выбираются из четырех возможностей: вверх, вниз, влево, вправо. Вероятность каждого направления равна $1/4$:

$$ P(\text{вверх}) = P(\text{вниз}) = P(\text{влево}) = P(\text{вправо}) = \frac{1}{4} $$

Для бессеточной модели:

Угол движения $\alpha$ выбирается равномерно из интервала $[0, 2\pi)$. Смещение на шаге $h$ задается как:

$$ \Delta x = h \cdot \cos(\alpha), \quad \Delta y = h \cdot \sin(\alpha) $$

Фрактальная размерность

Для количественной характеристики получаемых разветвленных структур используется понятие фрактальной размерности $D$.

Метод: зависимость массы от радиуса гирации

Масса кластера (число частиц) $N$ связана с его характерным размером степенным законом:

$$ N \sim R_g^D $$

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

$$ R_g = \sqrt{ \frac{1}{N} \sum_{i=1}^{N} \left[ (x_i - x_c)^2 + (y_i - y_c)^2 \right] } $$

где $(x_c, y_c)$ — координаты центра масс.

В логарифмических координатах эта зависимость становится линейной:

$$ \ln N = D \cdot \ln R_g + C $$

Тангенс угла наклона прямой, построенной методом наименьших квадратов, дает значение размерности $D$.

Метод: подсчет клеток (Box counting)

Пространство, содержащее кластер, покрывается сеткой с квадратными ячейками размера $\varepsilon$. Подсчитывается число ячеек $N(\varepsilon)$, содержащих хотя бы одну точку кластера. При уменьшении размера ячейки число непустых ячеек растет по закону:

$$ N(\varepsilon) \sim \varepsilon^{-D} $$

В логарифмических координатах:

$$ \ln N(\varepsilon) = -D \cdot \ln \varepsilon + C $$

По наклону графика $\ln N(\varepsilon)$ от $-\ln \varepsilon$ определяется размерность $D$.

Вероятностные правила прилипания

В базовой модели частица прилипает при первом же касании кластера (вероятность $p = 1$). В более сложных модификациях вводится параметр вероятности прилипания $p < 1$. При проверке контакта генерируется случайное число $\xi$, равномерно распределенное на $[0, 1]$. Прилипание происходит при выполнении условия:

$$ \xi \leq p $$

В химически-ограниченной агрегации вероятность $p$ зависит от локального окружения — числа уже занятых соседних узлов $k$:

$$ p = p(k), \quad k = 1, 2, 3, 4 $$

Диффузия кластеров

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

$$ D(M) = D_0 \cdot M^{-\gamma} $$

где $\gamma = 1/2$ (для свободно-сочлененных цепочек) или $\gamma = 1/d_f$ (для фрактальных агрегатов).

Алгоритмы

Базовый алгоритм DLA на сетке

Входные данные:

• Размер сетки $L \times L$
• Число частиц $N_{total}$
• Вероятность прилипания $p$

Шаг 1. Инициализация. Создать пустую сетку $grid[L][L]$. Поместить затравочную частицу в центр $(L/2, L/2)$. Установить счетчик $N = 1$.

Шаг 2. Основной цикл. Пока $N < N_{total}$:

• Вычислить текущий максимальный радиус кластера $R_{max}$.
• Сгенерировать стартовую позицию на окружности радиусом $R_{start} = R_{max} + \Delta$: $\alpha = random(0, 2\pi)$, $x = L/2 + R_{start} \cdot \cos(\alpha)$, $y = L/2 + R_{start} \cdot \sin(\alpha)$.
• Установить флаг прилипания $stuck = False$.

Шаг 3. Блуждание. Пока не $stuck$ и расстояние от центра $< R_{kill}$ ( $\approx 2 \cdot R_{max} + \delta$ ):

• Выбрать случайное направление из ${\uparrow, \downarrow, \leftarrow, \rightarrow}$.
• Сделать шаг: $(x, y) \to (x+dx, y+dy)$.
• Если вышли за границы сетки — прервать блуждание.
• Проверить 4 соседние клетки.
• Если хотя бы одна занята и $random() < p$: $grid[x][y] = 1$, $N = N + 1$, $stuck = True$.

Шаг 4. Если частица не прилипла (ушла далеко) — запустить новую (вернуться к шагу 2).

Алгоритм расчета размерности методом радиуса гирации

Входные данные:

• Массив координат всех точек кластера
• История роста (пары $N_i, R_{g,i}$)

Алгоритм:

Для каждого размера кластера $N$ (с шагом $\Delta N$) вычислить:

• Центр масс: $x_c = \frac{1}{N} \sum x_i$, $y_c = \frac{1}{N} \sum y_i$
• Квадрат радиуса гирации: $R_g^2 = \frac{1}{N} \sum \left[ (x_i - x_c)^2 + (y_i - y_c)^2 \right]$
• Радиус гирации: $R_g = \sqrt{R_g^2}$
• Сохранить пару $(N, R_g)$

Построить график в координатах $X = \ln(R_g)$, $Y = \ln(N)$. Выполнить линейную регрессию $Y = D \cdot X + C$. Наклон $D$ — искомая фрактальная размерность.

Алгоритм расчета размерности методом ящиков

Входные данные: бинарное изображение кластера (матрица занятости).

Алгоритм:

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

Для размера ячейки $\varepsilon = L/2, L/4, L/8, \dots, \varepsilon_{min}$:

• Разбить область на ячейки $\varepsilon \times \varepsilon$
• Подсчитать число непустых ячеек $N_{box}$
• Сохранить пару $(\varepsilon, N_{box})$

Построить график в координатах $X = -\ln(\varepsilon)$, $Y = \ln(N_{box})$. Выполнить линейную регрессию $Y = D \cdot X + C$. Наклон $D$ — фрактальная размерность.

Алгоритм с переменной вероятностью прилипания

Модификация шага проверки прилипания в базовом алгоритме:

При обнаружении контакта с кластером подсчитать число занятых соседей $k$. Определить вероятность $p_k$ по таблице:

Если $random() < p_k$ — прилипание. Иначе — отступить на предыдущую позицию и продолжить блуждание.

Алгоритм бессеточной модели DLA

Координаты частиц — вещественные числа $(x, y \in \mathbb{R})$.

Начальная позиция: $R_{start} = R_{max} + 5 \cdot r$, $\alpha = random(0, 2\pi)$, $x = R_{start} \cdot \cos(\alpha)$, $y = R_{start} \cdot \sin(\alpha)$, где $r$ — радиус частицы.

Шаг блуждания: $\alpha = random(0, 2\pi)$, $x = x + h \cdot \cos(\alpha)$, $y = y + h \cdot \sin(\alpha)$, где $h$ — длина шага.

Проверка прилипания: для каждой точки кластера $(x_c, y_c)$ вычислить расстояние $d = \sqrt{(x - x_c)^2 + (y - y_c)^2}$. Если $d \leq d_{stick}$ ( $\approx 2 \cdot r$ ) — прилипание.

Уничтожение: если расстояние от центра $> 2 \cdot R_{max} + 20 \cdot r$ — частица уничтожается.

Алгоритм баллистической агрегации

Инициализация: создать подложку — линию занятых узлов при $y = 0$.

Для каждой новой частицы:

• $x = random(0, W)$ ($W$ — ширина области), $y = H_{start}$ (запуск с высоты)
• Пока $y > 0$:
  • Если узел $(x, y-1)$ занят: занять узел $(x, y)$, выход из цикла
  • $y = y - 1$ (падение вниз)
  • С вероятностью $p_{walk}$ ($\approx 0.2$): $x = x \pm 1$
• Если $y = 0$ (достигнута подложка): занять узел $(x, 0)$

Граница кластера определяется как множество точек с максимальной $y$-координатой для каждого $x$.

Алгоритм кластер-кластерной агрегации

Инициализация: создать $N$ одиночных частиц в случайных позициях. Каждая частица — отдельный кластер размера 1.

Эволюция (пока число кластеров $> 1$):

• Выбрать случайный кластер $i$
• Вычислить коэффициент диффузии: $D_i = D_0 / \sqrt{size_i}$
• Сгенерировать смещение: $dx = Gauss(0, D_i)$, $dy = Gauss(0, D_i)$
• Сместить все точки кластера на $(dx, dy)$
• Проверить пересечение с другими кластерами. Если расстояние между точками разных кластеров $< 2 \cdot r$:
  • Слить кластеры $i$ и $j$ в один
  • Удалить $i$ и $j$ из списка
  • Добавить объединенный кластер
• Если пересечений нет — обновить позицию кластера $i$

Результат — один большой фрактальный агрегат.

Обоснование выбора подходов

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

Метод стартовой окружности позволяет сократить время расчета в десятки раз, исключая бесполезное блуждание частиц на периферии.

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

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

Выводы

В рамках второго этапа проекта:

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

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