Журнал вычислительной математики и математической физики, 2023, T. 63, № 4, стр. 629-638

Кинетика агрегации при седиментации. влияние диффузии частиц

Н. В. Бриллиантов 14*, Р. Р. Загидуллин 12**, С. А. Матвеев 23***, А. П. Смирнов 2****

1 Сколковский институт науки и технологий
121205 Москва, Большой бульвар, 30, стр. 1, Россия

2 Факультет вычислительной математики и кибернетики МГУ
119991 Москва, ул. Колмогорова, 1, стр. 52, Россия

3 Институт вычислительной математики РАН им. Г.И. Марчука
119333 Москва, ул. Губкина, 8, Россия

4 Универститет Лестера
LE1 7RH Лестер, Университетская ул., Великобритания

* E-mail: n.brilliantov@skoltech.ru
** E-mail: zagidullinrishat@gmail.com
*** E-mail: matseralex@cs.msu.ru
**** E-mail: sap@cs.msu.ru

Поступила в редакцию 28.05.2022
После доработки 27.09.2022
Принята к публикации 15.12.2022

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

Аннотация

В работе исследуется кинетика агрегации седиментирующих частиц аналитически и численно при использовании уравнения переноса-диффузии. Агломерация, вызванная этими механизмами (диффузией и переносом), важна как для мелких частиц (например, для первичных частиц пепла или сажи в атмосфере), так и для крупных частиц одинакового или близкого размера, где пространственная неоднородность менее выражена. Аналитические результаты можно получить для малых и больших чисел Пекле, определяющих относительное соотношение диффузии и переноса. При малых числах (пространственная неоднородность существует преимущественно из-за диффузии), мы получаем выражение для скорости агрегации через разложение чисел Пекле. При больших числах Пекле, когда перенос является основным источником пространственной неоднородности, мы получаем скорость агрегации из баллистических коэффициентов. Комбинируя эти результаты, предлагается аппроксимация рациональной функцией для всего диапазона чисел Пекле. Оцениваются скорости агрегации, численно решая уравнение переноса-диффузии. При этом результаты численного моделирования хорошо согласуются с аналитической теорией для большого диапазона рассматриваемых чисел Пекле (разница между минимальным и максимальным рассматриваемыми числами составляет 4 порядка). Библ. 26. Фиг. 2.

Ключевые слова: ядро коагуляции, пространственная неоднородность, число Пекле.

1. ВВЕДЕНИЕ

Осаждение коагулирующих частиц происходит в многочисленных природных и промышленных процессах, см., например, [15]. В качестве ярких примеров можно назвать агрегацию пыли, пепла или сажи, попадающих в воздух, коагулирующих органических частиц, оседающих в воде (озера, реки, моря). Частицы, падающие в жидкость, подвергаются действию движущей силы (силы тяжести и силе Архимеда) и силы трения. После короткого переходного периода эти силы уравновешиваются, что приводит к установившейся скорости падающей частицы. Различные установившиеся скорости приводят к столкновению частиц, если их траектории пересекаются. Если частица достаточно мала, становится важной стохастическая сила, возникающая из-за молекулярных флуктуаций в окружающей жидкости. Она порождает случайную, диффузионную составляющую движения частиц. Например, диффузионное движение сравнимо с баллистическим для частиц размером менее $2$ микрометра [68] в воздухе, что соответствует первичным частицам сажи [9].

Точно так же диффузия частиц может вызвать столкновения. Если частицы связаны, что всегда имеет место для мелких частиц из-за молекулярных поверхностных сил, они слипаются, образуя совместный агрегат. Эти агрегаты могут в дальнейшем сталкиваться с другими агрегатами или частицами. По мере эволюции системы появляются все более крупные кластеры, и распределение агрегатов по размерам в системе постоянно расширяется. Будем характеризовать агрегаты, состоящие из $k$ частиц (мономеров), их концентрацией ${{n}_{k}}(t)$ и радиусом ${{R}_{k}}$, считая их сферическими. Соответственно, концентрация первичных частиц (мономеров) радиуса ${{R}_{1}}$ равна ${{n}_{1}}(t)$.

Скорость процесса коагуляции, т.е. количество агрегатных столкновений в единицу времени и единицу объема агрегатов размера $i$ и $j$, можно записать как ${{K}_{{ij}}}{{n}_{i}}{{n}_{j}}$ [1012]. Она пропорциональна концентрации частиц ${{n}_{i}}$ и ${{n}_{j}}$ и константе скорости ${{K}_{{ij}}}$, которая зависит от природы частиц, их размера и конкретного процесса агрегации. В классическом случае Смолуховского для диффундирующих частиц, равномерно распределенных в пространстве и в отсутствие внешних сил, эти константы имеют вид ${{K}_{{ij}}} = 4\pi \left( {{{D}_{i}} + {{D}_{j}}} \right)\left( {{{R}_{i}} + {{R}_{j}}} \right)$, где ${{D}_{i}}$ и ${{D}_{j}}$ – соответственно коэффициенты диффузии агрегатов $i$ и $j$ [1012]. Следовательно, агрегация обусловлена взаимной диффузией частиц.

Если все агрегаты будут двигаться с одинаковой установившейся скоростью, то кинетические константы, естественно, не изменятся. Однако разница в установившихся скоростях меняет механизм агрегации, поскольку частицы могут сталкиваться, обгоняя друг друга. Если частицы достаточно велики, чтобы диффузионным движением можно было пренебречь, константы скорости следуют из чисто баллистических траекторий частиц, что дает $K_{{ij}}^{{{\text{bal}}}} = \pi {{\left( {{{R}_{i}} + {{R}_{j}}} \right)}^{2}}\left| {{{{v}}_{i}} - {{{v}}_{j}}} \right|$, где ${{{v}}_{i}}$ и ${{{v}}_{j}}$ – суммарные установившиеся скорости, соответствующие скоростям агрегации Смолуховского для сдвигового течения [13], [14]. Для учета гидродинамических взаимодействий между падающими частицами эти скорости модифицируются с помощью коэффициента эффективности столкновения $E({{R}_{i}},{{R}_{j}})$, из которых получаем: [6], [15], [16],

(1)
$K_{{ij}}^{{{\text{bal}}}} = \pi {{({{R}_{i}} + {{R}_{j}})}^{2}}E({{R}_{i}},{{R}_{j}})\left| {{{{v}}_{i}} - {{{v}}_{j}}} \right|.$
Множитель $E({{R}_{i}},{{R}_{j}})$ может быть получен в результате численного решения сложной гидродинамической задачи для двух обтекаемых сфер. К сожалению, в настоящее время отсутствует явное выражение для этого фактора; он доступен в виде таблиц со значениями $E$ для различных размеров частиц и давления воздуха [6].

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

Относительная важность адвективного движения и диффузионного движения количественно определяется числом Пекле ${\text{Pe}}$ (см. определение ниже). В [17], [18] аналитически исследована кинетика агрегации седиментирующих частиц с сильным преобладанием диффузионной составляющей (малые значения ${\text{Pe}}$), а в работе [19] рассмотрен другой предел сильного преобладания адвективного движения (большие значения ${\text{Pe}}$). В работе [4] было проведено обширное численное моделирование с использованием метода решеточных уравнений Больцмана в сочетании с методом больших вихрей. Однако во всех этих работах не было дано аналитического выражения для ядра скорости ${{K}_{{ij}}}$, адекватного для всего диапазона ${\text{Pe}}$.

Целью настоящей работы является построение явных выражений для ядра коагуляции ${{K}_{{ij}}}$, описывающих агрегацию седиментирующих частиц, которые можно использовать для произвольных значений числа Пекле. Такие коэффициенты скорости особенно важны для процесса агломерации, когда значения ${\text{Pe}}$ очень малы для первичных частиц и становятся большими для зрелых агрегатов. Для получения этих выражений рассмотрим отдельно два предельных случая малого и большого ${\text{Pe}}$. Затем мы строим простое рациональное приближение, которое точно описывает кинетику агрегации для всего диапазона числа Пекле. В этом подходе мы пренебрегаем гидродинамическими взаимодействиями между частицами; для учета последних взаимодействий мы предлагаем простое обобщение нашей аппроксимации.

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

2. МОДЕЛЬ АГРЕГАЦИИ СЕДИМЕНТИРУЮЩИХ ЧАСТИЦ

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

Рассмотрим первую осевшую одиночную частицу массы $m$, радиуса $R$ и плотности ${{\rho }_{p}}$. Если частица достаточно мала, она совершает случайное движение, которое подчиняется уравнению Ланжевена:

(2)
$m\frac{{d{\mathbf{v}}}}{{dt}} = - \gamma {\mathbf{v}} + {{{\mathbf{F}}}_{{\text{b}}}} + {{{\mathbf{F}}}_{{{\text{st}}}}}.$
Здесь сила вязкого трения равна $ - \gamma {\mathbf{v}}$ с $\gamma = 6\pi \eta R$, где $\eta $ – вязкость жидкости. Она действует вместе с силой ${{{\mathbf{F}}}_{b}}$, объединяющей гравитацию и плавучесть $(4{\text{/}}3)\pi {{R}^{3}}({{\rho }_{p}} - {{\rho }_{f}})g$, где ${{\rho }_{f}}$ – плотность жидкости, $g$ – ускорение свободного падения. Стохастическая сила ${{{\mathbf{F}}}_{{{\text{st}}}}}(t)$ имеет нулевое среднее значение и дисперсию в соответствии с теоремой диссипации-флуктуации (FDT) [10]:
(3)
$\left\langle {{{{\mathbf{F}}}_{{{\text{st}}}}}(t)} \right\rangle = 0,\quad \left\langle {F_{{{\text{st}}}}^{\alpha }(t)F_{{{\text{st}}}}^{\beta }(t{\kern 1pt} ')} \right\rangle = \Gamma {{\delta }_{{\alpha \beta }}}\delta \left( {t - t{\kern 1pt} '} \right),\quad \Gamma = 2{{k}_{{\text{B}}}}T\gamma ,$
где $\alpha ,\beta = x,y,z$, ${{k}_{{\text{B}}}}$ – постоянная Больцмана, $T$ – температура. Для малых частиц инерционные эффекты, связанные с членом $m{\mathbf{\dot {v}}}$, пренебрежимо малы, что приводит к передемпфированному уравнению движения:
(4)
${\mathbf{v}} = {\mathbf{\dot {r}}} = {{\gamma }^{{ - 1}}}{{{\mathbf{F}}}_{{\text{b}}}} + {{\gamma }^{{ - 1}}}{{{\mathbf{F}}}_{{{\text{st}}}}}.$
В приведенном выше уравнении Ланжевена рассматривается случайная величина – координата частицы ${\mathbf{r}}$, а роль стохастической силы играет ${{\gamma }^{{ - 1}}}{{{\mathbf{F}}}_{{{\text{st}}}}}$. Соответственно коэффициент $\Gamma $ в рамках флуктуационно-диссипационной теоремы следует заменить на $\Gamma {\kern 1pt} ' = \Gamma {\text{/}}{{\gamma }^{2}} = 2{{k}_{{\text{B}}}}T{\text{/}}\gamma $. Соответствующая функция распределения $P({\mathbf{r}},t\,{\text{|}}\,{{{\mathbf{r}}}_{0}},0)$ дает вероятность пребывания частицы в точке ${\mathbf{r}}$ в момент времени $t$ при условии, что в начальный момент времени $t = 0$ находился в точке ${{{\mathbf{r}}}_{0}}$. Она подчиняется уравнению Фоккера–Планка [10], [20]:
(5)
$\frac{{\partial P}}{{\partial t}} = - {{\nabla }_{{\mathbf{r}}}} \cdot {\mathbf{U}}P + D{{\Delta }_{{\mathbf{r}}}}P,$
где ${\mathbf{U}} = {{\gamma }^{{ - 1}}}{{{\mathbf{F}}}_{{\text{b}}}}$ – скорость установившегося движения частицы под действием силы ${{{\mathbf{F}}}_{{\;rmb}}}$ и
$D = \frac{{\Gamma {\kern 1pt} '}}{2} = \frac{{{{k}_{{\text{B}}}}T}}{\gamma } = \frac{{{{k}_{{\text{B}}}}T}}{{6\pi \eta R}}$
есть коэффициент диффузии.

Теперь применим приведенный выше анализ для двух частиц с массами ${{m}_{1}}$ и ${{m}_{2}}$ и соответствующими радиусами ${{R}_{1}}$ и ${{R}_{2}}$. Соответствующие уравнения Ланжевена с избыточным демпфированием получаются следующие:

(6)
${{{\mathbf{\dot {r}}}}_{1}} = \gamma _{1}^{{ - 1}}{{{\mathbf{F}}}_{{{\text{b}},1}}} + \gamma _{1}^{{ - 1}}{{{\mathbf{F}}}_{{{\text{st}},1}}},\quad {{{\mathbf{\dot {r}}}}_{2}} = \gamma _{2}^{{ - 1}}{{{\mathbf{F}}}_{{{\text{b}},2}}} + \gamma _{2}^{{ - 1}}{{{\mathbf{F}}}_{{{\text{st}},2}}},$
где, как и ранее, ${{\gamma }_{{1{\text{/}}2}}} = 6\pi \eta {{R}_{{1{\text{/}}2}}}$ и объединенные силы ${{{\mathbf{F}}}_{{{\text{b}},1{\text{/}}2}}}$ и стохастические силы ${{{\mathbf{F}}}_{{{\text{st}},1{\text{/}}2}}}$ подчиняются тем же соотношениям, что и выше, при соответствующих значениях параметров частиц. Мы предполагаем, что эти стохастические силы не скоррелированы, $\langle {{{\mathbf{F}}}_{{{\text{b}},1}}}(t){{{\mathbf{F}}}_{{{\text{b}},2}}}(t{\kern 1pt} ')\rangle = 0$.

Из уравнений (6) получаем уравнение Ланжевена для межчастичного расстояния ${{{\mathbf{r}}}_{{12}}}$:

(7)
${{{\mathbf{\dot {r}}}}_{{12}}} = {{{\mathbf{U}}}_{{12}}} + {{{\mathbf{F}}}_{{{\text{st}},12}}},$
где
(8)
${{{\mathbf{U}}}_{{12}}} = \gamma _{1}^{{ - 1}}{{{\mathbf{F}}}_{{{\text{b}},1}}} - \gamma _{2}^{{ - 1}}{{{\mathbf{F}}}_{{{\text{b}},2}}} = {{{\mathbf{U}}}_{1}} - {{{\mathbf{U}}}_{2}}$
это относительная установившаяся скорость со стохастической силой
(9)
$\langle {{{\mathbf{F}}}_{{{\text{st}},12}}}\rangle = 0,\quad \left\langle {F_{{{\text{st}},12}}^{\alpha }(t)F_{{{\text{st}},12}}^{\beta }(t{\kern 1pt} ')} \right\rangle = {{\Gamma }_{{12}}}{{\delta }_{{\alpha \beta }}}\delta (t - t{\kern 1pt} '),\quad {{\Gamma }_{{12}}} = \frac{{2{{k}_{{\text{B}}}}T}}{{{{\gamma }_{1}}}} + \frac{{2{{k}_{{\text{B}}}}T}}{{{{\gamma }_{2}}}} = 2({{D}_{1}} + {{D}_{2}}).$
Аналогично, как и для случая одной частицы, мы можем записать уравнение Фоккера-Планка для плотности вероятности ${{P}_{{12}}}({{{\mathbf{r}}}_{{12}}},t\,{\text{|}}\,{{{\mathbf{r}}}_{{12,0}}},0)$ межчастичного расстояния ${{{\mathbf{r}}}_{{12}}}$:
(10)
$\frac{{\partial {{P}_{{12}}}}}{{\partial t}} = - {{\nabla }_{{{{{\mathbf{r}}}_{{12}}}}}} \cdot {{{\mathbf{U}}}_{{12}}}{{P}_{{12}}} + {{D}_{{12}}}{{\Delta }_{{{{{\mathbf{r}}}_{{12}}}}}}{{P}_{{12}}},$
где ${{D}_{{12}}} = \frac{1}{2}{{\Gamma }_{{12}}} = {{D}_{1}} + {{D}_{2}}$ – относительный коэффициент диффузии для пары частиц. Отметим, что это уравнение получено без учета гидродинамических взаимодействий между частицами.

Теперь мы будем использовать уравнение (10), чтобы найти скорость агрегации. Сначала выберем направление силы тяжести по оси $OZ$ и запишем ${{{\mathbf{U}}}_{{12}}}$ как ${{{\mathbf{U}}}_{{12}}} = (2{\text{/}}9)\Delta \rho g\left( {R_{1}^{2} - R_{2}^{2}} \right){\mathbf{k}} = {{{v}}_{{12}}}{\mathbf{k}}$, где $\Delta \rho = {{\rho }_{p}} - {{\rho }_{f}}$ и ${\mathbf{k}}$ – единичный вектор в направлении $z$. Далее рассмотрим поток частиц радиуса ${{R}_{2}}$ (частиц сорта “2”) на поверхность частиц радиуса ${{R}_{1}}$ (частиц сорта “1”). Пусть частица сорта 1 находится в начале координат и $P({{{\mathbf{r}}}_{{12}}},t\,{\text{|}}\,{{{\mathbf{r}}}_{{12,0}}},0)$ дает вероятность нахождения частиц сорта 2 на расстоянии ${{{\mathbf{r}}}_{{12}}}$ от центра частиц сорта 1 в момент времени $t$. Удобно работать с концентрацией ${{n}_{2}}({\mathbf{r}},t)$ частиц сорта 2, которая просто пропорциональна плотности вероятности. Чтобы найти скорость агрегации, необходимо вычислить стационарное течение частиц сорта 2 на поверхности частицы сорта 1 в начале координат [10], [17], [18]. Мы предполагаем, что как только частицы касаются друг друга в точке ${\text{|}}{{{\mathbf{r}}}_{{12}}}{\kern 1pt} {\text{|}} = {{R}_{1}} + {{R}_{2}}$, мгновенно образуется совместный агрегат, и частицы сорта 2 исчезают. Это означает, что концентрация частиц 2 на поверхности ${{r}_{{12}}} = ({{R}_{1}} + {{R}_{2}})$ равна нулю. С другой стороны, на очень большом расстоянии от частицы сорта 1 концентрация частиц сорта 2 равна ${{n}_{{2,\infty }}}$. Тогда уравнение (10) может быть преобразовано в стационарном состоянии в форму

(11)
$D\Delta n - {v}\frac{n}{z} = 0,\quad n(r = R) = 0,\quad n(r \to \infty ) = {{n}_{\infty }}.$
Чтобы не загромождать обозначения, мы опускаем здесь все нижние индексы (т.е. $D = {{D}_{{12}}}$, $n = {{n}_{2}}$, ${v} = {{{v}}_{{12}}}$, ${{n}_{\infty }} = {{n}_{{2,\infty }}}$ и ${\mathbf{r}} = {{{\mathbf{r}}}_{{12}}}$) и используем сокращенную запись $R = {{R}_{1}} + {{R}_{2}}$. Далее введем безразмерные переменные ${\mathbf{r}} \to R{\mathbf{r}}$ и $n \to {{n}_{\infty }}n$, сохранив для простоты те же обозначения, что и для размерных величин. Тогда в сферических координатах с осью $OZ$ по силе тяжести (т.е. по скорости ${{{\mathbf{U}}}_{{12}}}$) безразмерное уравнение имеет вид
(12)
$\Delta n - 2\mu {{\partial }_{z}}n = 0,\quad {{\partial }_{z}}n = \cos \theta \frac{{\partial n}}{{\partial r}} - \frac{{\sin \theta }}{r}\frac{{\partial n}}{{\partial \theta }},$
где константа
(13)
$\mu = \frac{{{v}R}}{{2D}}$
соответствует числу Пекле. Также, используя замену
$n = u(r,\theta ){{e}^{{\mu z}}} + 1,$
мы преобразуем приведенное выше уравнение в форму,
(14)
$\Delta u - {{\mu }^{2}}u = 0,$
(15)
$u(r \to \infty ) = 0,$
(16)
$u(r = 1,\theta ) = - {{e}^{{ - \mu \cos \theta }}},$
где мы выбираем сферические координаты с осью $OZ$ вдоль скорости ${{U}_{{12}}}$ и используем симметрию системы относительно азимутального угла $\phi $.

Разделение переменных $u(r,\theta ) = \mathcal{R}(r)\Theta (\theta )$ дает два уравнения:

(17)
${{r}^{2}}\frac{{{{\partial }^{2}}\mathcal{R}}}{{\partial {{r}^{2}}}} + 2r\frac{{\partial{ \mathcal{R}}}}{{\partial r}} - \left( {{{\mu }^{2}}{{r}^{2}} + n(n + 1)} \right)\mathcal{R} = 0,$
(18)
${{\Delta }_{{\theta \phi }}}\Theta + n(n + 1)\Theta = 0,$
где ${{\Delta }_{{\theta \phi }}}$ – угловая часть лапласиана и $n = 0,1, \ldots $. Решение первого уравнения (17) можно представить в виде $\mathcal{R}(r) = {{K}_{{n + \frac{1}{2}}}}(\mu r){\text{/}}\sqrt r $, что является модифицированной функцией Бесселя II рода разделенной на квадратный корень из $r$. Решение второго уравнения (18) – многочлены Лежандра, $\Theta (\theta ) = {{P}_{n}}(\cos \theta )$. Следовательно, общее решение уравнения (14) имеет вид
(19)
$u(r,\theta ) = \frac{1}{{\sqrt r }}\sum\limits_{n = 0}^\infty \,{{A}_{n}}{{P}_{n}}(\cos \theta ){{K}_{{n + \frac{1}{2}}}}(\mu r),$
где коэффициенты ${{A}_{n}}$ следуют из граничных условий уравнений (15) и (16) [21]. Разлагая граничное условие (16) по полиномам Лежандра и используя $r = 1$ в уравнении (19), находим
(20)
${{A}_{n}} = - \frac{{(2n + 1)}}{{2{\kern 1pt} {{K}_{{n + \frac{1}{2}}}}(\mu )}}{{F}_{n}},$
(21)
${{F}_{n}} = \int\limits_{ - 1}^1 \,{{e}^{{ - \mu x}}}{{P}_{n}}(x)dx.$
Радиальная часть потока
(22)
${{J}_{r}} = - D\frac{n}{r} + {v}{\mathbf{k}} \cdot {{{\mathbf{e}}}_{r}}n = - D\frac{n}{r} + {v}\cos \theta n,$
где ${{{\mathbf{e}}}_{r}}$ – единичный вектор в радиальном направлении. На границе $r = R$ получаем поток
(23)
${{J}_{r}} = - \frac{{{{n}_{\infty }}D}}{R}\left[ {\frac{1}{2} - \mu \cos (\theta ) - \mu {{e}^{{\mu \cos (\theta )}}}{\kern 1pt} \sum\limits_{n = 0}^\infty \frac{{(2n + 1)K_{{n + \frac{1}{2}}}^{'}(\mu )}}{{2{{K}_{{n + \frac{1}{2}}}}(\mu )}}{{F}_{n}}{{P}_{n}}(\cos \theta )} \right],$
где $K_{{n + \frac{1}{2}}}^{'}$ обозначает производную функции. Коэффициент коагуляции можно получить, интегрируя радиальный поток по поверхности контакта двух частиц. Отсюда получаем
(24)
$K = \frac{1}{{{{n}_{\infty }}}}\int\limits_0^{2\pi } \,d\phi \int\limits_0^\pi \,{{R}^{2}}\sin \theta d\theta ( - {{J}_{r}})H( - {{J}_{r}}) = 2\pi DR\left[ {1 - \mu \sum\limits_{n = 0}^\infty \frac{{(2n + 1)K_{{n + \frac{1}{2}}}^{'}(\mu )}}{{2{{K}_{{n + \frac{1}{2}}}}(\mu )}}{\kern 1pt} {{F}_{n}}{{C}_{n}}} \right],$
где
(25)
${{C}_{n}} = \int\limits_0^\pi \,{{P}_{n}}(\cos \theta ){{e}^{{\mu \cos \theta }}}H\left[ { - {{J}_{r}}(\cos \theta )} \right]\sin \theta d\theta .$
Единичная ступенчатая функция Хевисайда $H(x)$ в подынтегральном выражении уравнений (24) и (25) гарантирует, что только радиальная составляющая потока, направленная к центру частицы, дает вклад в скорость коагуляции $K$. Используя соотношение $xK_{{1{\text{/}}2}}^{'}(x){\text{/}}{{K}_{{1{\text{/}}2}}}(x) = - \left( {x + 1{\text{/}}2} \right)$, легко проверить, что в случае, когда $\mu = 0$ получается стандартный результат Смолуховского $K = 4\pi RD$.

Трудности вычисления коэффициентов ${{C}_{n}}$ связаны с множителем $H( - {{J}_{r}})$ под интегралом. Наш численный анализ показал, что при малых числах Пекле $\mu < 0.5$, когда преобладает диффузия, поток отсутствует, т.е. ${{J}_{r}} < 0$ для всех углов $\theta $. Это позволяет убрать множитель $H( - {{J}_{r}})$ в уравнении (25) для малых $\mu $. Более того, при малых $\mu $ можно также использовать в уравне-ниях (21) и (25) для ${{F}_{n}}$ и ${{C}_{n}}$ разложения Тейлора, ${{e}^{{ \pm \mu x}}} = 1 \pm \mu x + {{(\mu x)}^{2}}{\text{/}}2 + \ldots $. Это дает следующее аналитическое выражение для $K$ при малых $\mu $:

(26)
$K = 4\pi DR\left( {1 + \mu - \frac{1}{3}{{\mu }^{2}} + \ldots } \right).$
Для очень больших $\mu $, когда преобладает перенос, можно пренебречь диффузионным движением и аппроксимировать поток выражением $J = {v}{\mathbf{k}}{{n}_{\infty }}$ с ${{J}_{r}} = {v}{{n}_{\infty }}\cos \theta $. Тогда из уравнения (24) получаем
(27)
$K = 2\pi \int\limits_0^\pi \,{{R}^{2}}{v}{\kern 1pt} \sin \theta \cos \theta H( - \cos \theta )d\theta = \pi {{R}^{2}}{v} = 2\pi RD\mu .$
Стоит найти экстраполяционное выражение, которое могло бы описать не только предельные случаи $\mu \to 0$ и $\mu \to \infty $, но и весь диапазон $\mu $. Это можно сделать, применив приближение, которое правильно воспроизводит предельные случаи и, как мы надеемся, обеспечивает приемлемую точность для промежуточных значений $\mu $.

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

(28)
$\mathcal{F}(\mu ) = \frac{{{{A}_{0}} + {{A}_{1}}\mu + {{A}_{2}}{{\mu }^{2}}}}{{1 + B\mu }}.$
Если поделить числитель и знаменатель на $\mu $ и допустить, что $\mu \to \infty $. Тогда для выполнения условия (27) надо, чтобы
(29)
$\frac{{{{A}_{2}}}}{B} = \frac{1}{2}.$
Для получения других уравнений, которые позволили бы однозначно определить коэффициенты, разложим знаменатель в ряд Тейлора до квадратичного члена. Тогда при соответствующих степенях и в соответствии с (26) получаем следующие соотношения:

(30)
${{A}_{0}} = 1,$
(31)
${{A}_{1}} - {{A}_{0}}B = 1,$
(32)
${{A}_{2}} - B{{A}_{1}} + {{A}_{0}}{{B}^{2}} = - \frac{1}{3}.$

Тогда мы находим соответствующие коэффициенты и подставляем их в ядро коагуляции:

(33)
$K = 4\pi DR\mathcal{F}(\mu ),\quad \mathcal{F}(\mu ) = \left( {\frac{{1 + \frac{5}{3}\mu + \frac{1}{3}{{\mu }^{2}}}}{{1 + \frac{2}{3}\mu }}} \right).$

Также стоит записать ядро коагуляции в начальных размерных обозначениях:

(34)
${{K}_{{ij}}} = 4\pi {{D}_{{ij}}}({{R}_{i}} + {{R}_{j}})\mathcal{F}({{\mu }_{{ij}}}),$
(35)
${{\mu }_{{ij}}} = \frac{{\left| {{{{v}}_{i}} - {{{v}}_{j}}} \right|\left( {{{R}_{i}} + {{R}_{j}}} \right)}}{{2{{D}_{{ij}}}}},$
где
(36)
${{D}_{{ij}}} = ({{D}_{i}} + {{D}_{j}}),$
(37)
$\left| {{{{v}}_{i}} - {{{v}}_{j}}} \right| = \frac{2}{9}\left| {\Delta \rho } \right|g\left( {R_{i}^{2} - R_{j}^{2}} \right)$
суть коэффициент взаимной диффузии и стоксова относительная скорость седиментирующих частиц, при этом ${{D}_{l}} = {{k}_{{\text{B}}}}T{\text{/}}(6\pi \eta {{R}_{l}})$, $l = i,j$. Функция $\mathcal{F}(x)$ определена в уравнении (33).

При выводе аналитического выражения (34) мы пренебрегаем гидродинамическими взаимодействиями между агрегирующими частицами. Для количественной оценки гидродинамической коррекции необходимо решение уравнений гидродинамики для двух сфер в жидкости. Это сложно сделать даже численно. В то время как поправочный коэффициент $E({{R}_{i}},{{R}_{j}})$ для баллистических столкновений представлен в виде нескольких таблиц, анализ коэффициента взаимной диффузии был выполнен аналитически, см., например, [17], [18], [24]. К сожалению, все эти исследования основаны на разложениях гидродинамических сил при большом межчастичном расстоянии; аналитические выражения для сил для малых расстояний в настоящее время отсутствуют. Однако наиболее важным является воздействие гидродинамических сил в условиях близкого расположения частиц, непосредственно перед столкновением. Поэтому получение надежных аналитических оценок представляется маловероятным. Хотя вывод соответствующих гидродинамических поправок является сложной задачей, возможно, стоит получить хотя бы грубую аналитическую оценку этого влияния на кинетику агломерации.

Мы предлагаем следующую феноменологическую модификацию (34) с гидродинамическими поправками:

(38)
$K_{{ij}}^{{{\text{(h)}}}} = 4\pi {{D}_{{ij}}}({{R}_{i}} + {{R}_{j}})\left( {\frac{{1 + {{\mu }_{{ij}}} + \mu _{{ij}}^{3}{{E}_{{ij}}}}}{{1 + \mu _{{ij}}^{2}}}} \right),$
где ${{E}_{{ij}}} = E({{R}_{i}},{{R}_{j}})$ – коэффициент эффективности столкновений [6]. Для больших ${{\mu }_{{ij}}} \gg 1$, когда диффузией можно пренебречь, приведенное выше уравнение даст кинетические скорости (1). В то же время при ${{\mu }_{{ij}}} \ll 1$ скорости (38) будут совпадать со скоростями (34) с точностью до членов порядка $\mathcal{O}\left( {mu_{{ij}}^{2}} \right)$. Обратите внимание, что уравнение (38) не учитывает гидродинамические поправки для коэффициентов взаимной диффузии из-за отсутствия надежных аналитических выражений для гидродинамических сил ближнего поля.

Мы ожидаем, что приведенное выше выражение (38) будет особенно полезно для описания кинетики агрегации для частиц близкого размера, когда ${{R}_{i}} \approx {{R}_{j}}$, так как стандартное соотношение (1) дает нулевые скорости агрегации, что не может быть верно. Потому что при малой разности плотностей $\Delta \rho $ диапазон разности радиусов $R_{i}^{2} - R_{j}^{2}$, приводящий к малым относительным скоростям $\left| {{{{v}}_{i}} - {{{v}}_{j}}} \right| \to 0$, может быть достаточно большим.

3. РЕЗУЛЬТАТЫ ЧИСЛЕННЫХ ЭКСПЕРИМЕНТОВ

Чтобы проверить точность нашей теории для коэффициентов коагуляции, когда присутствуют как диффузия, так и перенос, мы численно решаем уравнение (11) для различных чисел Пекле $\mu = {v}R{\text{/}}(2D)$. Мы используем неявный метод конечных разностей для видоизмененных радиальной и угловой операторов. А именно:

(39)
${{\Delta }_{r}}n = \frac{1}{{{{r}^{2}}}}\frac{\partial }{{\partial r}}\left( {\frac{{{{r}^{2}}}}{g}\frac{{\partial (gn)}}{{\partial r}}} \right),$
(40)
${{\Delta }_{\theta }}n = \frac{1}{{{{r}^{2}}\sin \theta }}\frac{\partial }{{\partial \theta }}\left( {\frac{{\sin \theta }}{g}\frac{{\partial (gn)}}{{\partial \theta }}} \right),$
(41)
$g = {{e}^{{ - 2\mu r\cos \theta }}}.$

Это преобразование необходимо для стабилизации численной ошибки, возникающей при больших числах Пекле $\mu $. Включение экспоненциального множителя $g$ приводит к членам $\mu {\kern 1pt} (\cos {{\theta }_{j}} - \cos {{\theta }_{{j \pm 1/2}}})$ и $\mu {\kern 1pt} dr$ (нижние индексы относятся к узлам сетки), которые могут уменьшить дестабилизирующий эффект больших $\mu $ за счет применения меньших шагов дискретизации [25].

Таким образом получается следующая численная схема:

(42)
$\begin{gathered} {{\Lambda }_{r}} = \frac{1}{{r_{i}^{2}}}\frac{1}{{{{h}_{r}}}}\left( {{{r}_{i}}{{r}_{{i + 1}}}\frac{{{{n}_{{i + 1,j}}}{{e}^{{ - 2\mu \cos {{\theta }_{j}}\left( {{{r}_{{i + 1}}} - {{r}_{{i + 0.5}}}} \right)}}} - {{n}_{{i,j}}}{{e}^{{ - 2\mu \cos {{\theta }_{j}}\left( {{{r}_{i}} - {{r}_{{i + 0.5}}}} \right)}}}}}{{{{h}_{r}}}}} \right. - \\ \, - \left. {{{r}_{i}}{{r}_{{i - 1}}}\frac{{{{n}_{{i,j}}}{{e}^{{ - 2\mu \cos {{\theta }_{j}}\left( {{{r}_{i}} - {{r}_{{i - 0.5}}}} \right)}}} - {{n}_{{i - 1,j}}}{{e}^{{ - 2\mu \cos {{\theta }_{j}}\left( {{{r}_{{i - 1}}} - {{r}_{{i - 0.5}}}} \right)}}}}}{{{{h}_{r}}}}} \right), \\ \end{gathered} $
(43)
$\begin{gathered} {{\Lambda }_{\theta }} = \frac{1}{{r_{i}^{2}\sin {{\theta }_{j}}}}\frac{1}{{{{h}_{\theta }}}}\left( {\sin {{\theta }_{{j + 0.5}}}\frac{{{{e}^{{ - 2\mu {{r}_{i}}\left( {\cos {{\theta }_{{j + 1}}} - \cos {{\theta }_{{j + 0.5}}}} \right)}}}{{n}_{{i,j + 1}}} - {{e}^{{ - 2\mu {{r}_{i}}\left( {\cos {{\theta }_{j}} - \cos {{\theta }_{{j + 0.5}}}} \right)}}}{{n}_{{i,j}}}}}{{{{h}_{\theta }}}}} \right. - \\ \;\left. { - \sin {{\theta }_{{j - 0.5}}}\frac{{{{e}^{{ - 2\mu {{r}_{i}}\left( {\cos {{\theta }_{j}} - \cos {{\theta }_{{j - 0.5}}}} \right)}}}{{n}_{{i,j}}} - {{e}^{{ - 2\mu {{r}_{i}}\left( {\cos {{\theta }_{{j - 1}}} - \cos {{\theta }_{{j - 0.5}}}} \right)}}}{{n}_{{i,j - 1}}}}}{{{{h}_{\theta }}}}} \right), \\ \end{gathered} $
(44)
${{\Lambda }_{r}} + {{\Lambda }_{\theta }} = 0,$
где $i = \{ 1, \ldots ,N - 2\} $, $j = \{ 0, \ldots ,M - 1\} $, ${{h}_{r}}$ – шаг по радиальной переменной, ${{h}_{\theta }}$ – шаг по угловой переменной. При $i = 0$, $i = N - 1$ схема заменяется граничными условиями для уравнения (11). Если разбить сетку по $\theta $ таким образом, что ${{\theta }_{{ - 0.5}}} = 0$ ($j = 0$) и ${{\theta }_{{M - 1 + 0.5}}} = \pi $ ($j = M - 1$), благодаря множетилям $\sin {{\theta }_{{j - 0.5}}}$ и $\sin {{\theta }_{{j + 0.5}}}$, угловой оператор остается в границах расчетной области по $n$ без необходимости явно определять граничные условия для численной схемы по угловой переменной.

Наконец, из полученной схемы мы строим систему линейных уравнений с разреженной матрицей размера $NM$ на $NM$, где $N$ – количество точек по $r$, $M$ – количество точек по $\theta $, решение которой дает нам численное решение уравнения (11). Для решения систем линейных уравнений мы использовали стандартный прямой решатель из пакета umfpack [26] через интерфейс scipy.sparse.linalg на языке программирования Python3.

Мы численно исследуем следующий диапазон числа Пекле, $\mu \in [0,01,100,0]$, используя сетку из 100 точек для угловой координаты $\theta $ и 500 точек для радиальной координаты $r$. Погрешность численной схемы по норме Чебышёва во всех расчетах не превышала $0,01$. Результат для распределения концентрации $n$ показан на фиг. 1 для трех различных чисел Пекле. Мы видим, что диффузионное “гало” исчезает при больших $\mu $ и существует при меньших $\mu $. Данный результат соотносится с интуитивными ожиданиями, что при больших числах Пекле диффузия не дает никакого вклада в коагуляцию.

Фиг. 1.

Численное решение стационарного уравнения переноса-диффузии (11) для различных чисел Пекле $\mu = {v}R{\text{/}}(2D)$.

На фиг. 2 представлены результаты для скорости агрегации $K$, полученные прямым численным решением уравнения (11), сравнивается с аналитическим результатом (33). Фигура демонстрирует близость аналитического приближения для всего диапазона чисел Пекле, охватывающего четыре порядка величины $\mu $.

Фиг. 2.

Приведенная скорость агрегации $K{\text{/}}(4\pi RD)$ как функция числа Пекле $\mu = {v}R{\text{/}}(2D)$. Результат, полученный прямым численным решением уравнения переноса-диффузии (11) (сплошная зеленая линия), очень близок к аналитической аппроксимации (33) для всего диапазона $\mu $. Также показаны асимптотики (26) для $\mu \to 0$ (сплошная синяя линия) и (27) для $\mu \to \infty $ (сплошная коричневая линия). Аналитическое выражение близко аппроксимирует численные результаты.

4. ЗАКЛЮЧЕНИЕ

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

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

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

  1. Vowinckel B., Withers J., Luzzatto-Fegiz P., Meiburg E. Settling of cohesive sediment: particle-resolved simulations // J. Fluid Mech. 2019. V. 858. P. 5–44.

  2. Fischer A., Chatterjee A., Speck T. Aggregation and sedimentation of active Brownian particles at constant affinity // J. Chem. Phys. 2019. V. 150. № 6. P. 064910.

  3. Yang Y.-J., Kelkar A.V., Corti D.S., Franses E.I. Effect of Interparticle Interactions on Agglomeration and Sedimentation Rates of Colloidal Silica Microspheres // Langmuir 2016. V. 32. № 20 P. 5111–5123.

  4. Hongsheng Ch., Wenwei L., Zhiwei Ch., Zhong Zh. A numerical study on the sedimentation of adhesive particles in viscous fluids using LBM-LES-DEM // Powder Technol. 2021. V. 391. P. 467–478.

  5. Whitmer J.K., Luijten E. Sedimentation of aggregating colloids // J. Chem. Phys. 2011. V. 134. P. 034510.

  6. Pinsky M., Khain A., Shapiro M. Collision Efficiency of Drops in a Wide Range of Reynolds Numbers: Effects of Pressure on Spectrum Evolution // J. Atmos. Sci. 2001. V. 58. P. 742–766.

  7. Khodzher T.V., et al. Study of Aerosol Nano-and Submicron Particle Compositions in the Atmosphere of Lake Baikal During Natural Fire Events and Their Interaction with Water Surface // Wat. Air and Soil Poll. 2021. V. 232. P. 266.

  8. Zhamsueva G., et al. Studies of the Dispersed Composition of Atmospheric Aerosol and Its Relationship with Small Gas Impurities in the Near-Water Layer of Lake Baikal Based on the Results of Ship Measurements in the Summer of 2020 // Atmosphere. 2022. V. 13. P. 139.

  9. Shahad H.A.K. An experimental investigation of soot particle size inside the combustion chamber of a diesel engine // Energy Convers. Manag. 1989. V. 29. P. 141–149.

  10. Krapivsky P.L., Redner A., Ben-Naim E. A Kinetic View of Statistical Physics. Cambridge: Cambridge University Press, 2010.

  11. Leyvraz F. Scaling theory and exactly solved models in the kinetics of irreversible aggregation // Phys. Rep. 2003. V. 383. P. 95–212.

  12. Smoluchowski M.V. Attempt for a mathematical theory of kinetic coagulation of colloid solutions // Z. Phys. Chem. 1917. V. 92. P. 265–271.

  13. Higashitani K., Ogawa R., Hosokawa G., Matsuno Y. Kinetic Theory of Shear Coagulation for Particles in a Viscous Fluid // J. Chem. Eng. Japan. 1982. V. 15. № 4. P. 299–304.

  14. Saffman P.F., Turner N.F. On the collision of drops in turbulent clouds // J. Fluid Mech. 1956. V. 1. № 1. P. 16–30.

  15. Falkovich G., Fouxon A., Stepanov, M. Acceleration of rain initiation by cloud turbulence // Nature. 2002. V. 419. P. 151.

  16. Falkovich G., Stepanov M.G., Vucelja M. Rain Initiation Time in Turbulent Warm Clouds // J. Appl. Meteor. Climatol. 2006. V. 45. P. 591.

  17. van de Ven T.G.M., Mason S.G. The microrheology of colloidal dispersions VIII. Effect of shear on perikinetic doublet formation // Colloid Polym. Sci. 1977. V. 255. № 8. P. 794–804.

  18. Melik D.H., Fogler H.S. Effect of gravity on Brownian flocculation // J. Colloid Interface Sci. 1984. V. 101. № 1. P. 84–97.

  19. Feke D.L., Scjowalter W.R. The effect of Brownian diffusion on shear-induced coagulation of colloidal dispersions // J. Fluid Mech. 1983. V. 133. P. 17–35.

  20. van Kampen N. Stochastic Processes in Physics and Chemistry. Amsterdam: Elsevier, 1992.

  21. Tikhonov A.N., Samarsky A.A. Equations of mathematical physics. New York: Dover Publications, 2013.

  22. Andrianov I., Shatrov A. Padé Approximants, Their Properties, and Applications to Hydrodynamic Problems // Symmetry. 2021. V. 13. P. 1869.

  23. Brezinski C. History of Continued Fractions and PadГ© Approximants. Berlin: Springer, 1991.

  24. Reed C.C., Anderson J.L. Hindered settling of a suspension at low Reynolds number // AIChE J. 1980. V. 26. № 5. P. 816–827.

  25. Samarskii A.A., Vabishchevich P.N Computational heat transfer. New York: John Wiley and Sons, 1995.

  26. Davis T.A. Algorithm 832 // ACM Trans. Math. Softw. 2004 V. 30. № 2 P. 196–199.

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