К теории роста системы кристаллов в переохлаждeнных / пересыщенных жидкостях

Полный текст

ВВЕДЕНИЕ

Объемные фазовые превращения из метастабильного и неравновесного состояния широко распространены в природе (замерзание воды, застывание лавы), а также часто используются в лабораторных и промышленных установках для кристаллизации металлических расплавов и пересыщенных растворов (например, затвердевание жидких капель в электромагнитных левитаторах) [1–7]. Важным примером практического применения являются кристаллизаторы с непрерывным перемешиванием для производства фармацевтических препаратов и пищевых добавок с особыми свойствами [8–11].

Все эти процессы основаны на эволюции полидисперсного ансамбля частиц на начальных и промежуточных стадиях объемной кристаллизации, когда размеры кристаллов малы, а среднее расстояние между ними намного больше их характерного размера. Кинетические и балансовые уравнения, описывающие эти стадии процесса объемной кристаллизации, зависят от скорости роста отдельных кристаллов, которые в свою очередь зависят от текущего размера частиц (радиуса или объема частиц) и переохлаждения/пересыщения жидкости [12–16].

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

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

Для того, чтобы учесть нестационарность теплового/концентрационного поля вокруг сферического зародыша, была разработана аналитическая теория на основе приближенного решения задачи Стефана с использованием метода дифференциальных рядов [17–19]. Таким образом была определена скорость роста отдельной частицы в тепловом поле [17], концентрационном поле [18] и термоконцентрационном поле [19].

Найденные законы для скорости роста зародыша зависят от времени и переохлаждения/пересыщения жидкости. Устранив явную зависимость от времени с помощью соответствующего выражения для радиуса частицы, ученые нашли скорость роста сферических кристаллов в зависимости от их радиуса и текущего значения переохлаждения/пересыщения жидкости [20, 21]. Эти законы существенно отличаются от ранее известных квазистационарных формул и имеют предельные переходы к ним.

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

УРАВНЕНИЯ МОДЕЛИ

Рассмотрим однокомпонентную переохлажденную или пересыщенную систему, в которой происходят нуклеация и рост частиц. Частота (темп) нуклеации является основным кинетическим параметром уменьшения степени метастабильности жидкости.

Методы расчета частоты нуклеации частиц, а также сама молекулярно-кинетическая теория зарождения рассматриваются во многих работах (см., например, [24, 25]). Молекулярно-кинетическая теория нуклеации основана на предположении, что в результате флуктуаций плотности в метастабильных расплавах или растворах могут образовываться микроядра новой фазы. Впоследствии эти ядра становятся способными к дальнейшему росту, когда их размер превышает критический радиус r_*.

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

$W\left(r\right)=4\pi r^2{\gamma }_i-\frac{4}{3}\pi r^3{\rho }_s\left({\mu }_l-{\mu }_s\right),$                                                                (1)

где: γ_i - поверхностное натяжение, ρ_s – плотность стабильной твердой фазы, и μ_l > μ_s – химические потенциалы исходной жидкой и стабильной твердой фаз.

Для метастабильной жидкой системы μ_l > μ_s , следовательно,W(r) достигает максимума, соответствующего некоторому критическому размеру частиц r_*. Поскольку W(r) равно изменению термодинамического потенциала Гиббса, который минимален в стабильном состоянии, частицы с радиусом r < r_* неустойчивы и исчезают, а частицы с радиусом r > r_* неустойчивы и растут. Ядра с радиусом r = r_* находятся в неустойчивом равновесии с матричной средой. Работа их образования W_* следует из уравнения (1) при r = r_*:

$W_{*}=\frac{16{\text{πγ}}_i^3}{3{\text{ρ}}_s^2{\left({\text{μ}}_l-{\text{μ}}_s\right)}^2}.$                                                                                   (2)

Разность между химическими потенциалами может быть рассчитана из уравнения Гиббса–Гельмгольца [26]:

${\text{μ}}_l-{\text{μ}}_s=\frac{L\Delta \vartheta }{{\vartheta }_p}, \Delta \vartheta ={\vartheta }_p-{\vartheta }_l,$                                                                   (3)

где L — скрытая теплота фазового перехода на единицу массы вещества, ϑ_l и ϑ_p – температуры расплава и фазового перехода.

Возникновение ядер новой фазы в метастабильной системе можно рассматривать как прохождение через энергетический барьер высотой W_*, который препятствует зарождению. В этом случае скорость нуклеации I может быть представлена как экспоненциальная функция высоты энергетического барьера [25–27]. Следовательно, используя уравнения (2) и (3), имеем

$I=I_{*}\exp \left(\frac{-W_{*}}{k_B{\vartheta }_l}\right)=I_{*}\text{exp}\left(\frac{-16{\text{πγ}}_i^3{\vartheta }_p^2}{3{\text{ρ}}_s^2{L}^2\Delta {\vartheta }^2{k}_B{\vartheta }_l}\right),$                                                    (4)

где k_B — постоянная Больцмана.

Можно также отметить, что предэкспоненциальный фактор I_* не сильно зависит от переохлаждения/пересыщения [24]. Предэкспоненциальный коэффициент может быть функцией r [28], но в рассматриваемой теории, где скорость нуклеации является частью граничного условия, I_*

 можно считать постоянным. Более того, поскольку $\Delta \vartheta \ll {\vartheta }_p$, уравнение (4) может быть переписано как

$I=I_{*}\exp \left(\frac{-16{\text{πγ}}_i^3{\vartheta }_p}{3{\text{ρ}}_s^2{L}^2\Delta {\vartheta }^2{k}_B}\right).$                                                                               (5)

Упомянутая модель справедлива для однокомпонентных расплавов и может быть модифицирована в случае бинарных расплавов, поскольку изменение свободной энергии в результате зарождения зависит от состава расплава. Как показано в работах [29–31], поверхностной примесной диффузией растущих зародышей можно пренебречь, и скорость нуклеации может быть записана аналогично уравнению (5). В этом случае температура фазового перехода ϑ_p, входящая в Δϑ, зависит от концентрации C_l примеси, и поэтому ϑ_p следует заменить функцией ϑ_p (C_l) . Для удобства введем начальное переохлаждение Δϑ₀ и запишем уравнение (5) в виде

$I=I_{*}\exp \left[\frac{-p}{{\left(\Delta \vartheta /\Delta {\vartheta }_0\right)}^2}\right], p=\frac{16 {\text{π γ}}_i^3{\vartheta }_p}{3 {\text{ρ}}_s^2{L}^2\Delta {\vartheta }^2 k_B},$                                                       (6)

где p — безразмерное число Гиббса, соответствующее начальному переохлаждению Δϑ₀, а уравнение (6) дает скорость нуклеации как функцию относительного переохлаждения ϖ = Δϑ /Δϑ₀.

В случае пересыщенных растворов выражение, аналогичное уравнению (6), имеет вид [32],

$I=I_{*}\exp \left[\frac{-p}{{\ln }^2\left(C_l/C_p\right)}\right], p=\frac{16 {\text{π γ}}_i^3{M}_s^2}{3 {\text{ρ}}_s^2{R}_g^2 {\vartheta }_s^3 k_B}$,                                           (7)

где C_p — концентрация при насыщении, M_s — молекулярная масса, R_g — универсальная газовая постоянная, ϑ_s — температура раствора.

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

Далее будем использовать формулу, полученную в работах [20, 21] для скорости роста сферического кристалла в однокомпонентном переохлажденном расплаве. Это уравнение является обобщением, включающим в себя нестационарность теплового поля вокруг сферической частицы и зависимость температуры кристаллизации от кривизны межфазной границы:

$\frac{dr}{d\text{τ}}=\widetilde{{\text{β}}_{\text{*}}}\left(\Delta \vartheta -\frac{\text{χ}}{r}\right)\sqrt{1-2\widetilde{{\text{β}}_{\text{*}}}qr}, q=\frac{{\text{ρ}}_{s}L}{{\text{λ}}_l},$                                                (8)

где τ — время роста частицы, $\widetilde{{\text{β}}_{\text{*}}}={\text{β}}_{*}/\left(1+{\text{β}}_{\text{*}}/{\text{μ}}_k\right)$, $\text{χ}={\vartheta }_p\text{α}/\left({\text{ρ}}_{s}L\right)$, ${\text{β}}_{*}$ и ${\text{μ}}_k$ — кинетические коэффициенты, а α — коэффициент поверхностного натяжения, λ_l — коэффициент теплопроводности (см. подробнее ссылки в работах [20, 21]).

Уравнение (8) показывает, что скорость роста является нелинейной функцией переохлаждения расплава Δϑ. Для очень маленьких кристаллов, когда $r\ll 1/\left(2\widetilde{{\text{β}}_{\text{*}}}q\right)$, скорость роста не зависит от их размера для χ = 0.

В случае пересыщенных растворов выражение (8) преобразуется таким образом:

$\frac{dr}{d\text{τ}}=\widetilde{{\text{β}}_{\text{*}}}\left(\Delta C-\frac{\text{χ}}{r}\right)\sqrt{1-2\widetilde{{\text{β}}_{\text{*}}}qr}, q=\frac{C_p}{D},$                                                (9)

где C_p обозначает концентрацию насыщения.

Уравнение (9) показывает ту же зависимость скорости роста, но для пересыщенных систем, где ΔC и D — соответственно пересыщение и коэффициент диффузии.

Обратим внимание, что существуют теории роста, которые используют выражения, отличные от уравнений (8) и (9) [8–11]. Использование этих выражений позволяет изучить взаимодействие между зарождением и эволюцией ансамблей частиц в метастабильных средах. Теория, разработанная далее, может быть легко скорректирована для использования различных выражений скорости роста частиц.

Рассмотрим суспензию твердых сферических частиц в макроскопически однородной переохлажденной однокомпонентной или пересыщенной жидкости. Эволюция такой двухфазной среды описывается кинетическим уравнением типа Фоккера–Планка для функции распределения частиц по радиусу $\text{ψ}\left(r,\text{τ}\right)$ и уравнением теплового баланса для переохлаждения расплава Δϑ. Пренебрегая случайными флуктуациями в скоростях роста частиц и предполагая, что распределение температуры не зависит от пространственных координат, имеем

$\frac{\partial \text{ψ}}{\partial \text{τ}}+\frac{\partial }{\partial r}\left(\frac{dr}{d\text{τ}}\text{ψ}\right)=\mathrm{0,} r>r_{*},$                                                            (10)

${\text{ρ}}_m{C}_m\frac{d\Delta \vartheta }{d\text{τ}}=-4{\text{πρ}}_{s}L\underset{r_{*}}{\overset{\infty }{\int }}{r}^2\frac{dr}{d\text{τ}}\text{ψ}dr,$                                                      (11)

где ρ_m и C_m — постоянные плотность и теплоемкость, а τ обозначает время процесса эволюции системы частиц. Здесь мы предполагаем, что жидкость, которая ранее не содержала кристаллов, почти мгновенно охлаждается в нулевой момент времени до температуры ниже температуры кристаллизации на величину Δϑ₀. Тогда граничные условия, накладываемые на решения уравнений (10) и (11), имеют вид

$\text{ψ}=\mathrm{0,} \Delta \vartheta =\Delta {\vartheta }_0, \text{τ}=\mathrm{0,}$                                                                    (12)

$\frac{dr}{d\text{τ}}\text{ψ}=I\left(\Delta \vartheta \right), r=r_{*}.$                                                                          (13)

Последнее условие определяет поток кристаллов, который проходит критический барьер нуклеации. Выражения (4)–(13) представляют собой общую систему уравнений для описания начальной и промежуточной стадий объемной кристаллизации.

КИНЕТИКА СНЯТИЯ ПЕРЕОХЛАЖДЕНИЯ/ПЕРЕСЫЩЕНИЯ

Введем, а также следующие безразмерные функции и параметры

$t=\frac{\text{τ}}{{\text{τ}}_0}, s=\frac{r}{l_0}, \varpi =\frac{\Delta \vartheta }{\Delta {\vartheta }_0}, \text{Ψ}=l_0^4\text{ψ}, {\text{τ}}_0={\left(\widetilde{{\text{β}}_{*}}\Delta {\vartheta }_0^3{I}_0\right)}^{-1/4},$            (14)

$l_0={\left(\widetilde{{\text{β}}_{\text{*}}}\Delta {\vartheta }_0\right)}^{1/4}{I}_0^{-1/4}, {\text{γ}}_{*}=\frac{2\widetilde{{\text{β}}_{*}}{\text{ρ}}_{s}L{l}_0}{{\text{λ}}_l}, b_1=\frac{4{\text{πρ}}_{s}L}{{\text{ρ}}_m{C}_m\Delta {\vartheta }_0}, {\text{α}}_{*}=\frac{\text{χ}}{l_0\Delta {\vartheta }_0}.$

Здесь t и s — безразмерные временная и пространственная переменные, ϖ — безразмерное переохлаждение, Ψ — безразмерная функция распределения частиц по радиусу, τ и l₀ — характерные масштабы времени и длины, а безразмерные параметры b₁ и ${\text{α}}_{*}$ описывают соответственно эффекты скрытой теплоты кристаллизации и поверхностного натяжения.

Упомянутая математическая модель может быть переписана в безразмерной форме с помощью выражений (14). Для простоты математических преобразований рассмотрим случай $r_{*}\to 0$. Тогда получим

$\frac{\partial \text{Ψ}}{\partial t}+\frac{\partial }{\partial s}\left(\text{Ψ}\left(\varpi -\frac{{\text{α}}_{\text{*}}}{s}\right)\sqrt{1-{\text{γ}}_{*}s}\right)=\mathrm{0,}$                                                (15)

$\frac{d\varpi }{dt}=-b_1\underset{0}{\overset{\infty }{\int }}\text{Ψ}{s}^2\left(\varpi -\frac{{\text{α}}_{\text{*}}}{s}\right)\sqrt{1-{\text{γ}}_{*}s} ds,$                                              (16)

$\text{Ψ}=\mathrm{0,} \varpi =\mathrm{1,} t=\mathrm{0,}$                                                                      (17)

$\text{Ψ}=\frac{1}{\varpi -\frac{{\text{α}}_{*}}{s}}\exp \left[pg\left(\varpi \right)\right], s=0.$                                                         (18)

Уравнения модели (15)–(18) имеют место для переохлажденных и пересыщенных жидкостей, а функция g(ϖ) в граничном условии (18) определяет кинетику нуклеации. Например, рассматривая кинетику Вебера–Вольмера–Френкеля–Зельдовича (Weber–Volmer–Frenkel–Zel’dovich, WVFZ) для переохлажденных расплавов (подстрочный индекс sm), получаем из уравнения (6)

$g\left(\varpi \right)=g_{sm}=1-{\varpi }^{-2}.$                                                                     (19)

При анализе различных экспериментальных данных часто используется кинетика нуклеации Майера (Meirs) [33]. В этом случае, $I=I_{*}{\left(\Delta \vartheta \right)}^p$, и

$g\left(\varpi \right)=g_{sm}=\ln \varpi .$                                                                          (20)

Отметим, что параметр p для кинетики Майера должен быть найден из экспериментальных данных. Также следует учитывать, что ϖ — относительное переохлаждение, когда речь идет о переохлажденных расплавах.

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

$t=\frac{\text{τ}}{{\text{τ}}_0}, s=\frac{r}{l_0}, \varpi =\frac{\Delta C}{\Delta C_0}, \text{Ψ}=l_0^4\text{ψ}, {\text{τ}}_0={\left(\widetilde{{\text{β}}_{\text{*}}}\Delta C_0^3{I}_0\right)}^{-1/4},$            (21)

$l_0={\left(\widetilde{{\text{β}}_{\text{*}}}\Delta C_0\right)}^{1/4}{I}_0^{-1/4}, {\text{γ}}_{*}=\frac{2\widetilde{{\text{β}}_{*}}{C}_p{l}_0}{D}, b_1=\frac{4\text{π}{C}_p}{\Delta C_0}, {\text{α}}_{*}=\frac{\text{χ}}{l_0\Delta C_0}.$

Здесь ϖ — безразмерное пересыщение, ΔC₀— начальное пересыщение раствора,  описывает эффект концентрации при насыщении, а остальные безразмерные переменные и параметры описаны ранее.

Аналогично уравнению (19) для переохлажденных расплавов из уравнения (7) получаем следующее выражение для пересыщенных растворов, используя кинетику WVFZ [33]:

$g\left(\varpi \right)=g_{ss}={\ln }^{-2}\left(1+\frac{1}{{\varpi }_p}\right)-{\ln }^{-2}\left(1+\frac{\varpi }{{\varpi }_p}\right),$                                       (22)

где ${\varpi }_p=C_p/\Delta C_0$. В случае кинетики Майера для пересыщенных растворов $I=I_{*}{\left(\Delta C\right)}^p$, что дает [33]

$g\left(\varpi \right)=g_{ss}=\ln \varpi .$                                                                            (23)

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

Особо подчеркнем, что выражения (15) и (16) являются сильно нелинейными интегро-дифференциальными уравнениями, которые управляют кинетикой переохлаждения/пересыщения. Решение этой модели основано на методе седловой точки для интеграла типа Лапласа [15, 33]. Используя эту технику, мы продемонстрируем аналитический метод решения рассматриваемой модели.

Для простоты рассмотрим случай пренебрежимо малого α_*. Теория, учитывающая этот эффект, была разработана в статье [34]. Точное аналитическое решение уравнения (15), с условиями (17) и (18), имеет вид

$\text{Ψ}\left(t,s\right)=\frac{1}{\sqrt{1-{\text{γ}}_{*}s}}\text{φ}\left(x_1\left(t\right)-s_1\left(s\right)\right)\text{η}\left(x_1\left(t\right)-s_1\left(s\right)\right),$                                  (24)

где

$x_1\left(t\right)=\underset{0}{\overset{t}{\int }}\varpi \left(t_2\right)d{t}_2, s_1\left(s\right)=\underset{0}{\overset{s}{\int }}\frac{d{s}_2}{\sqrt{1-{\text{γ}}_{*}{s}_2}}=-\frac{2\sqrt{1-{\text{γ}}_{*}s}}{{\text{γ}}_{*}}+\frac{2}{{\text{γ}}_{*}},$

$\text{φ}\left(x_1\left(t\right)\right)=\frac{1}{\varpi \left(x_1\left(t\right)\right)}\exp \left[pg\left(\varpi \left(x_1\left(t\right)\right)\right)\right],$

где $x_1$ — модифицированное время, а $\text{η}\left(x_1\right)$ — функция Хевисайда.

Интегрирование безразмерной скорости роста кристаллов дает

$\frac{ds}{dt}=\varpi \sqrt{1-{\text{γ}}_{*}s}.$

Этот закон роста получен для переохлажденных расплавов с помощью выражений (8) и (14). Имея дело с пересыщенными растворами, приходим к тому же результату, используя выражения (9) и (21). В обоих случаях получаем безразмерный радиус частицы s как функцию t и v с учетом начального условия s = 0 при t = v:

${\text{γ}}_{*}s=1-{\left[1-\frac{{\text{γ}}_{*}}{2}\left(x_1\left(t\right)-x_1\left(\nu \right)\right)\right]}^2.$                                                      (25)

При v = 0 выражение (25) определяет радиус $s_m\left(t\right)$ самых крупных частиц, существующих в метастабильной жидкости.

Заменяя теперь  при любом фиксированном t на новую переменную v, получаем

$x_1\left(\nu \right)=x_1\left(t\right)-s_1\left(s\right).$

Это выражение показывает, что пределы интегрирования $s = 0$ и $s=s_m\left(t\right)$ относительно t соответствуют пределам интегрирования v = t и v = 0 относительно v. Теперь, комбинируя (24), (25) и (16), имеем:

$\frac{d\varpi }{dt}=-b_1\varpi \underset{0}{\overset{t}{\int }}h\left(\nu ,t\right)\exp \left[pg\left(\nu \right)\right]d\nu ,$                                                   (26)

где

$h\left(\nu ,t\right)={\left[\frac{1}{{\text{γ}}_{*}}\left(1-{\left[1-\frac{{\text{γ}}_{\text{*}}}{2}\left(x_1\left(t\right)-x_1\left(\nu \right)\right)\right]}^2\right)\right]}^2\left(1-\frac{{\text{γ}}_{\text{*}}}{2}\left(x_1\left(t\right)-x_1\left(\nu \right)\right)\right),$       (27)

а $g\left(\nu \right)=g\left(\varpi \left(\nu \right)\right)$ определяется формулами (19), (20), (22) и (23).

Важным моментом является то, что $dg/d\nu <0$ для всех рассматриваемых кинетик нуклеации. Этот вывод следует из выражений (19), (20), (22) и (23). Например, рассматривая упомянутые кинетики для переохлажденных расплавов, получаем $dg/d\nu =2{\varpi }^{-3}d\varpi /d\nu <0$ (WVFZ) и $dg/d\nu ={\varpi }^{-1}d\varpi /d\nu <0$ (Meirs). Это означает, что g(v) имеет точку максимума при v = 0 (в граничной точке).

Оценивая теперь производные по v с помощью выражения (26), получаем, что три первых из них равны нулю при v = 0 для всех кинетик нуклеации. Однако n-я производная $g^{\left(n\right)}{|}_{\left(\nu =0\right)}\ne 0$ при n ≥ 4. Первые ненулевые производные перечислены в табл. 1, где

$\varkappa ={\left(\frac{4}{b_1}\right)}^{1/4}, k=\frac{1}{\left(1+{\varpi }_p\right){\ln }^3\left(1+{\varpi }_p^{-1}\right)}.$

 

Таблица 1. Параметры и функции аналитического решения для различных кинетических механизмов

Примечание: z = z(x₁ (t)) = $\frac{{\gamma }_{*}}{2}$x_i(t), а Ψ₀, Ψ₁₅, Ψ₂ и Ψ₃ приведены в прил. A.

Оценивая интеграл типа Лапласа (26) с помощью метода седловой точки [22, 23], получаем

$\underset{0}{\overset{t}{}}h\left(\nu ,t\right)\exp \left[pg\left(\nu \right)\right]d\nu \approx \sum _{k=0}^{\infty }{p}^{-\left(k+1\right)/4)}{a}_k\left(x_1\left(t\right)\right),$                                    (28)

$a_k\left(x_1\left(t\right)\right)={\left(-1\right)}^{k+1}\frac{4^k}{k!}\mathit{\text{Γ}}\left(\frac{k+1}{4}\right)H^k\left(0\right)\frac{d^k}{d{\nu }^k}{\left[h\left(\nu ,t\right)H\left(\nu \right)\right]}_{\nu =0},$                  (29)

где Г(‧) обозначает гамма-функцию Эйлера, а

$H\left(\nu \right)=\frac{{\left[g\left(\nu \right)\right]}^{3/4}}{g\text{'}\left(\nu \right)}.$

Первые четыре коэффициента  приведены в прил. А.

Подставляя

$\frac{d\varpi }{dt}=\frac{d\varpi }{d{x}_1}\frac{d{x}_1}{dt}=\varpi \frac{d\varpi }{d{x}_1}$

в левую часть выражения (26) с учетом (28), получаем дифференциальное уравнение для безразмерной степени метастабильности ϖ (переохлаждения или пересыщения)

$\frac{d\varpi }{d{x}_1}=-b_1\sum _{k=0}^{\infty }{p}^{-\left(k+1\right)/4}{a}_k\left(x_1\right).$

Интегрируя это уравнение, имеем

$\varpi \left(x_1\right)=1-b_1{p}^{-1/4}\underset{0}{\overset{x_1}{\int }}\sum _{k=0}^{\infty }{p}^{-k/4}{a}_k\left(x_2\right)d{x}_2.$                                       (30)

Теперь, комбинируя (30) и $\varpi =d{x}_1/dt$, находим

$t\left(x_1\right)=\underset{0}{\overset{x_1}{\int }}\frac{d{x}_2}{\varpi \left(x_2\right)}.$                                                                           (31)

Выражения (30) и (31) определяют полное аналитическое решение для ϖ в параметрической форме (x₁ представляет собой параметр решения). Заметим, что ϖ₀ — основной (фундаментальный) вклад, а ϖ = ϖ_i (i  ≥ 1) дает i-е приближение, учитывающее поправки к ϖ₀ (эти функции выписаны в прил. A и табл. 1). Для i-го приближения имеем

$\varpi ={\varpi }_i=1-b_1{p}^{-1/4}\underset{0}{\overset{x_1}{\int }}\sum _{k=0}^i{p}^{-k/4}{a}_k\left(x_2\right)d{x}_2.$

На рис. 1–5 показано аналитическое решение для степени метастабильности ϖ и функции распределения Ψ частиц по радиусам (теплофизические параметры приведены в табл. 2).

 

Рис. 1. Безразмерная степень метастабильности ϖ = Δϑ/Δϑ₀ в зависимости от безразмерного времени t = τ/τ₀ : а – кинетика нуклеации WVFZ для Δϑ₀ = 373K, р = 66.9 и γ_* = 1.66·10⁻⁴; б – кинетика нуклеации WVFZ для Δϑ₀ = 573K, р = 7.4 и  γ_* = 2.18·10⁻⁴; в – кинетика нуклеации Meirs для Δϑ₀= 573K, р = 4 и γ_* = 2.18·10⁻⁴.

 

Рис. 2. Безразмерная степень метастабильности ϖ = Δϑ/Δϑ₀ в зависимости от безразмерного времени t = τ/τ₀ : а – кинетика нуклеации WVFZ для Δϑ₀ = 373K, р = 66.9 и γ_* = 1.66·10⁻¹; б – кинетика нуклеации WVFZ для Δϑ₀ = 573K, р = 7.4 и γ_* = 2.18·10⁻¹; в – кинетика нуклеации Meirs для Δϑ₀= 573K, р = 4 и γ_* = 2.18·10⁻¹.

   

Табл. 2. Теплофизические параметры

Параметр	Значение
Поверхностное натяжение, γi	0.3 Дж/м2
Температура фазового перехода, ϑp	1 273 K
Плотность твердой фазы, ρs	7.8103 кг/м3
Скрытая теплота, L	0.9106 Дж/кг
Плотность, ρ1, ρm	7103 кг/м3
Коэффициент теплопроводности,  λl	63 Дж/м·К·с
Кинетический коэффициент, β	0.510–7 м/с·K
Теплоемкость, Cm	840 Дж/кг·К
Частота нуклеации, I0	1011 1/м3·с

 

Расчеты показывают, что основной вклад ϖ₀ имеет правильную тенденцию, но численно отличается от полного решения уравнения (30) с поправками (отличается от решений ϖ₂ и ϖ₃). Принимая во внимание ϖ₁, ϖ₂ и ϖ₃, приходим к выводу, что учет третьей (ϖ₂) и четвертой (ϖ₃) поправок дает почти совпадающее решение. Поэтому с практической точки зрения мы можем ограничиться третьим приближением ϖ₂.

Наши расчеты показывают, что этот вывод справедлив и для достаточно малых значений p ∼10^-1. В этом случае ряд (30) быстро сходится из-за малых значений γ_* (которые порядка 10^–2 – 10^–4 для реальных метастабильных расплавов) при любых фиксированных x₁ и t.

Аналитические решения (24) и (30) позволяют найти функцию распределения частиц по радиусам помня, что максимальный размер частиц задается формулой

$s_m\left(x_1\left(t\right)\right)=\frac{1-{\left[1- \frac{{\text{γ}}_{*}}{2} x_1 \left(t\right) \right]}^2}{{\text{γ}}_{*}}.$                                                    (32)

Таким образом, функция распределения имеет вид:

$\text{Ψ}\left(x_1\left(t\right),s\right)=\frac{1}{\sqrt{1-{\text{γ}}_{*}s}}{E}_1^{-1} \left(x_1 \left(t\right),s\right) \mathrm{ex}p\left[-\frac{p \left(1- E_1^2\right)}{E_1^2}\right],$                 (33)

где

$E_1\left(x_1\left(t\right),s\right)=1-b_1$$\underset{0}{\overset{x_1\left(t\right)+\frac{2\sqrt{1-{\text{γ}}_{*}s}}{{\text{γ}}_{*}}-\frac{2}{{\text{γ}}_{*}}}{}}$$\sum _{k=0}^{\infty }{p}^{-\left(k+1\right)/4}{a}_k\left(x_2\right)d{x}_2.$

Важным моментом является то, что переменная s всегда меньше или равна s_m, а x₁(t) определяется (31) как обратная функция. i-я аппроксимация функции распределения (33) приведена в прил. B.

 

Рис. 3. Функция распределения Ψ(s,t) в зависимости от безразмерного радиуса кристаллов s = r/l₀ в различные элементы t = τ/τ₀ для кинетики нуклеации WVFZ, где Ψᵢ(s,t) = Ψ₃(s,t): a – Δϑ₀ = 373K, p = 66.9, (1) – x₁ = 0.05, t = 0.05, (2) – x₁ = 0.25, t = 0.254, (3) – x₁ = 0.45, t = 0.493, (4) – x₁ = 0.625, t = 1.489; б – Δϑ₀, (1) – x₁ = 0.05, t = 0.05, (2) – x₁ = 0.25, t = 0.255, (3) – x₁ = 0.45, t = 0.475, (4) – x₁ = 0.742, t = 2.038. Вертикальные линии показывают максимальный размер sm частиц для различных моментов времени. Линии (4), соответствуют полностью переохлажденной жидкости (ϖ = 0).

 

На рис. 3 показана функция распределения частиц по радиусам, построенная с учетом четырех коэффициентов a_k (x₁)  (см. прил. A). Эта функция эволюционирует со временем и ограничена максимальным радиусом s_m (t) сферических кристаллов.

В моменты, близкие к начальному времени, функция распределения содержит кристаллы размером 0 < s < s_m (t). Кристаллы растут, выделяют скрытую теплоту кристаллизации и частично компенсируют переохлаждение жидкости. Поэтому процесс их зарождения практически прекращается, и происходит практически только укрупнение существующих частиц в менее переохлажденной жидкости (функция распределения смещается в сторону больших размеров кристаллов).

Полученные аналитические решения позволяют найти общее число N(x₁(t))  частиц в метастабильной жидкости и их средний размер L̅(x₁(t))  в момент :

$N\left(x_1\left(t\right)\right)=\frac{1}{l_0^3}\underset{0}{\overset{s_m\left(x_1\left(t\right)\right)}{\int }}\text{Ψ}\left(x_1\left(t\right),s\right)ds,$                                                     (34)

$\overline{L}\left(x_1\left(t\right)\right)=l_0\underset{0}{\overset{s_m\left(x_1\left(t\right)\right)}{\int }}s\text{Ψ}ds{\left(\underset{0}{\overset{s_m\left(x_1\left(t\right)\right)}{\int }}\text{Ψ}ds\right)}^{-1}.$                                          (35)

i-е приближения этих функций выписаны в прил. B.

Особо отметим влияние нестационарности в скоростях роста частиц (параметр γ_*) на эволюцию ансамбля кристаллов. Во-первых, при увеличении γ_* скорость роста частиц $ds/dt=\varpi \sqrt{1-{\text{γ}}_{*}s}$ уменьшается, кристаллы выделяют меньше теплоты кристаллизации, а степень метастабильности (переохлаждение или пересыщение) уменьшается медленнее.

 

Рис. 4. Функция распределения Ψ(s,t) в зависимости от безразмерного радиуса кристаллов s = r/l₀ в различные элементы t = τ/τ₀ для кинетики нуклеации WVFZ, где Ψᵢ(s,t) = Ψ₃(s,t): Δϑ₀ = 373K, p = 66.9 (1) – x₁ = 0.625, t = 1.489, (2) – x₁ = 0.639, t = 1.489. Вертикальные линии показывают максимальный размер sm частиц для различных времен.

 

i-е приближения для безразмерного числа кристаллов и их среднего радиуса показаны на рис. 5. 

 

Рис. 5. Относительное количество кристаллов nᵢ = Nᵢ(x₁)/N₃(x_*) (a) и их средний радиус uᵢ = L̅ᵢ(x₁)/L̅₃(x*) (б) в зависимости от безразмерного времени t = τ/τ₀ для кинетики нуклеации WVFZ, ∆i₀ = 573K , p = 7.4 . Здесь нижний индекс “_*” обозначает полностью переохлажденное состояние жидкости (ϖ = 0, x_* = 0.742, t_* = 2.038).

 

Как легко видеть, третье и четвертое приближения совпадают, т. е. n₂ ≈ n₃, и u₂ ≈ u₃. Это снова означает быструю сходимость решений и возможность использования третьего приближения для практических расчетов.

Выпишем теперь переохлаждение расплава Δϑ и функцию распределения частиц по радиусу Ψ в размерной форме:

$\Delta \vartheta =\Delta {\vartheta }_i\left(x_1\left(\text{τ}\right)\right)=\Delta {\vartheta }_0-\frac{4{\text{πρ}}_{s}L{p}^{-1/4}}{{\rho }_m{C}_m}\underset{0}{\overset{x_1\left(\frac{\text{τ}}{{\text{τ}}_0}\right)}{\int }}\sum _{k=0}^i{p}^{-k/4}{a}_k\left(x_2\right)d{x}_2,$            (36)

$\text{ψ}={\text{ψ}}_i\left(x_1\left(\text{τ}\right),r\right)=\frac{1}{l_0^4\sqrt{1-\frac{2\widetilde{{\text{β}}_{*}}{\text{ρ}}_{s}Lr}{{\text{λ}}_l}}}{E}_1^{-1}\left(x_1\left(\frac{\text{τ}}{{\text{τ}}_0}\right),\left(\frac{r}{l_0}\right)\right)\exp \left[-\frac{p\left(1-E_1^2\right)}{E_1^2}\right],$  (37)

$\text{τ}={\text{τ}}_0\underset{0}{\overset{x_1\left(\frac{\text{τ}}{{\text{τ}}_0}\right)}{\int }}\frac{\Delta {\vartheta }_{0}d{x}_2}{\Delta \vartheta \left(x_2\right)}.$                                                                        (38)

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

ЗАКЛЮЧЕНИЕ

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

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

Наши расчеты показывают, что первые четыре приближения по методу седловой точки дают сходящееся решение для рассматриваемых кинетик нуклеации (WVFZ и Meirs). Функция распределения увеличивается вплоть до максимального размера кристаллов и смещается в сторону больших радиусов с увеличением времени. Нестационарность роста отдельных частиц существенно влияет на динамическое поведение всего кристаллического ансамбля.

Довольно часто зарождение и рост кристаллов сопровождаются направленным фазовым превращением вдоль выбранного пространственного направления. Движущей силой для направленного движения является градиент температуры и/или концентрации. В таких процессах зарождение и рост кристаллов могут проиодить в двухфазном слое, заполненном дендритоподобными структурами [35–41]. Поэтому имеет смысл расширить настоящую теорию на такие более сложные системы с подвижными границами области кристаллизации. Это представляет собой важное направление для будущих исследований, которое следует развивать путем обобщения настоящей теории и теории двухфазного слоя [42–51].

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

Развиваемая теория определяет начальное состояние расплавов и растворов на заключительной стадии фазового превращения, когда такие процессы, как оствальдово созревание, коагуляция и фрагментация кристаллов, являются преобладающими [25, 57–60].

Работа выполнена при финансовой поддержке Российского научного фонда (грант 23–19–00337).

ПРИЛОЖЕНИЕ A

Коэффициенты и функции, входящие в аналитическое решение:

$a_0\left(z\right)=-\text{Γ}\left(\frac{1}{4}\right)h\left(\mathrm{0,}t\right)H\left(0\right),$

$a_1\left(z\right)=4\text{Γ}\left(\frac{1}{2}\right)H\left(0\right)\left[H\left(0\right){\left(\frac{\partial h}{\partial \nu }\right)}_{\nu =0}+h\left(\mathrm{0,}t\right)H^{\text{'}}\left(0\right)\right],$

$a_2\left(z\right)=-8\text{Γ}\left(\frac{3}{4}\right)H^2\left(0\right)\left[2H^{\text{'}}\left(0\right){\left(\frac{\partial h}{\partial \nu }\right)}_{\nu =0}+H\left(0\right){\left(\frac{{\partial }^{2}h}{\partial {\nu }^2}\right)}_{\nu =0}+h\left(\mathrm{0,}t\right)H^{\text{'}\text{'}}\left(0\right)\right],$

$a_3\left(z\right)=\frac{32}{3}{H}^3\left(0\right)\left[3H^{\text{'}\text{'}}\left(0\right){\left(\frac{\partial h}{\partial \nu }\right)}_{\nu =0}+h\left(\mathrm{0,}t\right)H^{\text{'}\text{'}\text{'}}\left(0\right)+3H^{\text{'}}\left(0\right){\left(\frac{{\partial }^{2}h}{\partial {\nu }^2}\right)}_{\nu =0}+H\left(0\right){\left(\frac{{\partial }^{3}h}{\partial {\nu }^3}\right)}_{\nu =0}\right],$

$h\left(\mathrm{0,}t\right)=\frac{z}{{\text{γ}}_{*}^2}{\left(1-z^2\right)}^2, {\left(\frac{\partial h}{\partial \nu }\right)}_{\nu =0}=-\frac{2z^2}{{\text{γ}}_{*}}\left(1-z^2\right)+\frac{{\left(1-z^2\right)}^2}{2{\text{γ}}_{*}},$

${\left(\frac{{\partial }^{2}h}{\partial {\nu }^2}\right)}_{\nu =0}=5z^3-3z, {\left(\frac{{\partial }^{3}h}{\partial {\nu }^3}\right)}_{\nu =0}=\frac{15}{2}{\text{γ}}_{*}{z}^2-\frac{3}{2}{\text{γ}}_{*},$

${\text{ψ}}_0=\text{Γ}\left(\frac{1}{4}\right)\frac{\varkappa b_1}{4{\text{γ}}_{*}^3{p}^{-1/4}}\left(z^2-z^4+\frac{z^6-1}{3}\right),$

${\text{ψ}}_1=2\text{Γ}\left(\frac{1}{2}\right)\frac{b_1{\varkappa }^2}{{\text{γ}}_{*}^2{p}^{1/2}}\left(\frac{1}{6}+\frac{2z}{19}-\frac{z^2}{2}-\frac{4z^3}{19}+\frac{2z^5}{19}+\frac{z^4}{2}-\frac{z^6}{6}\right),$

${\text{ψ}}_2=\text{Γ}\left(\frac{3}{4}\right)\frac{b_1{\varkappa }^3}{{\text{γ}}_{\text{*}}{p}^{3/4}}\left(\frac{4}{100}+\frac{z}{15}-\frac{25z^2}{100}-\frac{2z^3}{15}+\frac{56z^4}{100}-\frac{z^5}{15}-\frac{17z^6}{100}\right),$

${\text{ψ}}_3=\frac{4}{3p}\left(-\frac{15}{100}-\frac{39z}{100}-\frac{6z^2}{10}-\frac{65z^3}{100}+\frac{13z^4}{10}+\frac{z^5}{100}-\frac{4z^6}{10}\right),$

ПРИЛОЖЕНИЕ B

i-е приближения для функции распределения частиц по радиусам, общего числа частиц и их среднего размера определяются следующими выражениями:

${\text{Ψ}}_i\left(x_1\left(t\right),s\right)=\frac{1}{\sqrt{1-{\text{γ}}_{*}s}}{E}_{1i}^{-1}\left(t,s\right)\exp \left[-\frac{p\left(1-E_{1i}^2\right)}{E_{1i}^2}\right],$

$E_{1i}\left(x_1\left(t\right),s\right)=1-b_1\underset{0}{\overset{x_1\left(t\right)+\frac{2\sqrt{1-{\text{γ}}_{*}s}}{{\text{γ}}_{\text{*}}}-\frac{2}{{\text{γ}}_{\text{*}}}}{\int }}\sum _{k=0}^i{p}^{-\left(k+1\right)/4}{a}_k\left(x_2\right)d{x}_2,$

$N_i\left(x_1\left(t\right)\right)=\frac{1}{l_0^3}\underset{0}{\overset{s_m\left(x_1\left(t\right)\right)}{\int }}{\text{Ψ}}_i\left(x_1\left(t\right),s\right)ds,$

$\overline{L_i}\left(x_1\left(t\right)\right)=l_0\underset{0}{\overset{s_m\left(x_1\left(t\right)\right)}{\int }}s{\text{Ψ}}_{i}ds{\left(\underset{0}{\overset{s_m\left(x_1\left(t\right)\right)}{\int }}{\text{Ψ}}_{i}ds\right)}^{-1},$

$n_i\left(x_1\left(t\right)\right)=\frac{N_i\left(x_1\left(t\right)\right)}{N_3\left(x_{*}\right)}, u_i\left(x_1\left(t\right)\right)=\frac{\overline{L_i}\left(x_1\left(t\right)\right)}{\overline{L_3}\left(x_{*}\right)},$

где i -0, 1, 2, ...

Об авторах

Е. В. Маковеева

Уральский федеральный университет им. Б.Н. Ельцина

Email: dmitri.alexandrov@urfu.ru
Россия, Екатеринбург

И. Е. Корозникова

Уральский федеральный университет им. Б.Н. Ельцина

Email: dmitri.alexandrov@urfu.ru
Россия, Екатеринбург

А. Е. Глебова

Уральский федеральный университет им. Б.Н. Ельцина

Email: dmitri.alexandrov@urfu.ru
Россия, Екатеринбург

А. А. Иванов

Уральский федеральный университет им. Б.Н. Ельцина

Email: dmitri.alexandrov@urfu.ru
Россия, Екатеринбург

М. А. Никишина

Уральский федеральный университет им. Б.Н. Ельцина

Email: dmitri.alexandrov@urfu.ru
Россия, Екатеринбург

Л. В. Торопова

Уральский федеральный университет им. Б.Н. Ельцина

Email: dmitri.alexandrov@urfu.ru
Россия, Екатеринбург

Д. В. Александров

Уральский федеральный университет им. Б.Н. Ельцина

Автор, ответственный за переписку.
Email: dmitri.alexandrov@urfu.ru
Россия, Екатеринбург

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

Дополнительные файлы