Линейный анализ устойчивости процесса затвердевания в ограниченной области при наличии конвекции в жидкости

ВВЕДЕНИЕ

Динамика межфазной границы (фронта кристаллизации), разделяющей твердую и жидкую фазы, определяет образование и рост элементов твердой фазы: например, рост дендритных кристаллов, перераспределение и поглощение примесей, а также рост и развитие различных микроструктур (это могут быть ячеистые, полосчатые, нерегулярные микроструктуры и микроструктуры смешанного типа) в кристаллизующемся материале [1$–$7]. Определив закономерности, которые отвечают за формирование различных типов микроструктур, а также за переходы между их различными типами, можно управлять процессом структурообразования с помощью физических и режимных параметров процесса затвердевания. Так, воздействовать на структурообразование материалов можно с помощью флуктуаций различных физических величин. Примерами флуктуаций могут быть: колебания температуры, постоянно присутствующие в природе; гидродинамические колебания скорости течения расплава, которые могут быть вызваны конвекцией; колебания концентрации примесей в жидкости, вызванные поступлением примесей извне или гидродинамическими колебаниями; механические колебания всей кристаллизующейся системы, вызванные сейсмическими процессами в природе или естественным фоном деятельности человека в лабораторной установке, и т.д. Такие флуктуации приводят к возмущениям температурных и градиентных полей, гидродинамическим возмущениям скорости движения жидкости, а это, в свою очередь, приводит к возмущению межфазной границы «твердое тело – жидкость» (фронта кристаллизации). Морфологические и динамические возмущения развиваются из-за неоднородности температурных и концентрационных полей, конвекции или концентрационного переохлаждения. В результате такие возмущения полностью меняют процесс затвердевания и приводят к формированию различных микроструктур. К примеру, морфологическая неустойчивость фронта кристаллизации влияет на развитие ячеистых структур, динамическая неустойчивость приводит к образованию полосчатых структур, а эволюция «дендритного леса» в целом порождает двухфазную область [8–13]. Линейный анализ морфологической устойчивости с приложениями к проблемам кристаллизации впервые был развит Маллинзом и Секеркой [14–16]. Затем их методика была расширена на описание различных особенностей процессов затвердевания [17–24].

Конвективные течения жидкости могут по-разному сказываться на процессе фазового превращения из жидкого состояния в твердое. Так, они могут выравнивать распределения температуры и концентрации в жидкости или, наоборот, создавать конвективные или холодные/горячие области в определенных местах у границы раздела «твердое тело – жидкость» [25, 26]. В результате этого в более холодных областях создаются благоприятные условия для роста твердой фазы и возникает морфологическая неустойчивость. На данный момент влияние конвекции на морфологическую/динамическую неустойчивость мало исследовано из-за сильной нелинейности математических моделей, разнообразия конвективных потоков и различных граничных условий на межфазной поверхности при наличии конвекции. Предыдущие исследования рассматривали либо определённые типы конвективных течений, либо течения в узком диапазоне параметров процесса. Например, для исследования влияния плоскопараллельного течения жидкости на устойчивость межфазной границы в работе [27], были применены упрощенные уравнения тепло- и массопереноса. В работах [28–30] авторы применили технику пограничного слоя для исключения малых нелинейных членов и упрощения уравнений модели. В работе [31] рассматривались пространственно-периодические течения жидкости для моделирования локализованных морфологических структур. Недавно для изучения влияния течений жидкости на устойчивость направленной кристаллизации была развита математическая модель в неограниченной области пространства, основанная на законах кондуктивного тепло- и массопереноса [32, 33]. Эта теория продемонстрировала, что дисперсионное соотношение и кривая нейтральной устойчивости при наличии конвекции существенно зависят от скорости расширения на границе раздела фаз.

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

Процесс затвердевания с плоской границей раздела твердой и жидкой фаз

Рассмотрим направленную кристаллизацию бинарной жидкости (расплава или раствора) с постоянной скоростью при наличии плоской межфазной границы раздела. Пусть ось $z$ совпадает с направлением затвердевания, которое вызвано определенными температурными градиентами в твердой и жидкой фазах. Рассматриваемый процесс начинается после охлаждения жидкости сверху и контролируется интенсивной конвекцией в жидкости, рис. 1а. Используя систему координат, которая движется с постоянной скоростью ${\upsilon }_s$ вместе с границей раздела фаз, имеем следующие уравнения теплопереноса в твердом теле и жидкости:

$\frac{\partial {\vartheta }_l}{\partial t}+u\cdot \nabla {\vartheta }_l=D_l{\nabla }^2{\vartheta }_l+{\upsilon }_s\frac{\partial {\vartheta }_l}{\partial z}, z>Z^{\prime }\left(t\right),$ (1)

 $\frac{\partial {\vartheta }_s}{\partial t}=D_s{\nabla }^2{\vartheta }_s+{\upsilon }_s\frac{\partial {\vartheta }_s}{\partial z},z<Z^{\prime }\left(t\right),$ (2)

где ${\vartheta }_l$ и ${\vartheta }_s$ $—$ температуры в жидкой и твердой фазах, $t$ - время, $u$ $—$ вектор скорости жидкости, $D_l$ и $D_s$ $—$ коэффициенты теплопроводности в жидкой и твердой фазах, ${\upsilon }_s$ $—$ скорость кристаллизации, $Z^{\prime }\left(t\right)$ $—$ координата границы раздела между твердым телом и жидкостью ( $Z^{\prime }\left(t\right)=0$ в случае невозмущенного стационарного режима затвердевания). При достаточно интенсивной конвекции в жидкости, которая выравнивает распределение концентрации, процесс диффузии примеси можно не рассматривать [34].

Рис. 1. а $–$ Схема морфологически неустойчивой кристаллизации с плоской границей раздела «твердое тело-жидкость»; б $–$ схема динамически неустойчивой кристаллизации с поверхностью разрыва.

На границе раздела фаз температура равна температуре фазового превращения, которая складывается из температуры фазового перехода чистого вещества ${\vartheta }_{*}$ и температурной поправки на кривизну границы раздела ГН. Также разность тепловых потоков определяет выделяющуюся теплоту кристаллизации $L_V$ [34, 35], то есть:

${\vartheta }_{int}={\vartheta }_l={\vartheta }_s={\vartheta }_{*}+ГН,z=Z^{\prime }\left(t\right),$ (3)

$\begin{array}{c}\left[c_l\left({\vartheta }_{\infty }-{\vartheta }_l\right)+L_V\right]u\cdot n=\\ =k_{s}n\cdot \nabla {\vartheta }_s-f\left({\vartheta }_l\right), z=Z^{\prime }\left(t\right),\end{array}$ (4)

Где Г $={\vartheta }_{*}\gamma /L_V$ $—$ коэффициент Гиббса, H$—$ кривизна границы раздела, $\gamma$ $—$ межфазная энергия, $L_V$ $—$ параметр скрытой теплоты, $n$ $—$ вектор нормали к границе раздела, $c_l$ $—$ удельная теплота, ${\vartheta }_{\infty }$ $—$ температура жидкости вдали от фронта кристаллизации, $k_s$ $—$ теплопроводность твердого тела, $f\left({\vartheta }_l\right)$ $—$ конвективный тепловой поток. Следует отметить, что в случае плоской границы раздела между твердым телом и жидкостью H = 0, а при малых морфологических возмущениях (линейная теория) $Н\approx {\nabla }^2{Z}^{\prime }\left(t\right)$. Конвективный тепловой поток $f\left({\vartheta }_l\right)$ задается в следующем виде [34, 35]:

$f\left({\vartheta }_l\right)=2^{4/3}\lambda k_l{\left(\frac{\alpha \left({\vartheta }_l\right)g}{D_l\varkappa }\right)}^{1/3}{\left({\vartheta }_{\infty }-{\vartheta }_l\right)}^{4/3},$ (5)

где $\lambda$ $—$ безразмерная постоянная, $k_l$ $—$ теплопроводность жидкости, $D_l=k_l/c_l$, $\varkappa$ $—$ кинематическая вязкость, $\alpha$ $—$ коэффициент теплового расширения, $g$ $—$ ускорение свободного падения. Отметим, что $\alpha$ может быть функцией от ${\vartheta }_l$: $\alpha \left({\vartheta }_l\right)={10}^{-4}\left(2.25+0.15{\vartheta }_l\right)$ $–$ для изопропанола [34].

Морфологическая неустойчивость

Предположим, что теплообмен преобладает в направлении пространственной оси $z$. В этом случае установившиеся профили температуры в жидкой и твердой фазах имеют вид ${\vartheta }_{l0}={\vartheta }_{l0}\left(z\right)$ и ${\vartheta }_{s0}={\vartheta }_{s0}\left(z\right)$ (подстрочный индекс «0» обозначает установившиеся решения).

Из граничного условия (4) находим стационарный градиент температуры в твердом теле при $z=0$:

$\begin{array}{c}\frac{d{\vartheta }_{s0}}{dz}=G_s=\frac{\left[c_l\left({\vartheta }_{\infty }-{\vartheta }_{int}\right)+L_V\right]{\upsilon }_s+f\left({\vartheta }_{int}\right)}{k_s},\\ z=\mathrm{0,}\end{array}$ (6)

где ${\vartheta }_{int}={\vartheta }_{*}$ для плоской границы раздела твердой и жидкой фаз.

Используя уравнения теплопроводности (1) и (2), находим вторые производные при $z=0$ (здесь используется условие прилипания жидкости к твердому телу, т.е. $u=0$ при $z=0$ ) в виде

$\begin{array}{c}\frac{d^2{\vartheta }_{s0}}{d{z}^2}=-\frac{{\upsilon }_s{G}_s}{D_s}, \\ \frac{d^2{\vartheta }_{l0}}{d{z}^2}=-\frac{{\upsilon }_s{G}_l}{D_l},\\ z=\mathrm{0,}\end{array}$ (7)

где $G_l=d{\vartheta }_{l0}/dz$ при $z=0$.

Далее, морфологически возмущаем плоскую границу раздела «твердое тело – жидкость» в виде $Z^{\prime }=Z_{А}E\left(x,y,t\right)$ с $E=\text{exp}\left(i{k}_{x}x+i{k}_{y}y+\varpi t\right)$. Здесь $x$ и $y$ $—$ декартовы координаты, направленные перпендикулярно оси затвердевания $z$, $k_x$ и $k_y$ $—$ волновые числа возмущений вдоль этих направлений, $i$ $—$ мнимая единица, $\varpi$ $—$ частота возмущений. Возмущение границы раздела фаз $Z^{\prime }$ соответствует температурным возмущениям в твердой и жидкой фазах: ${\vartheta }_s^{\text{'}}={\vartheta }_{sA}E\text{exp}\left({\beta }_{s}z\right)$, ${\vartheta }_l^{\text{'}}={\vartheta }_{lA}E\text{exp}\left({\beta }_{l}z\right)$, где $Z_{А}$, ${\vartheta }_{sA}$ и ${\vartheta }_{lA}$ - амплитуды возмущений, а ${\beta }_s$ и ${\beta }_l$ - коэффициенты нарастания/затухания возмущений, которые найдены ниже. Из рассматриваемой линейной теории следует, что $\left|{\vartheta }_s^{\text{'}}\right|\ll {\vartheta }_{s0}$ и $\left|{\vartheta }_l^{\text{'}}\right|\ll {\vartheta }_{l0}$.

Далее разлагаем граничные условия (3) и (4) в ряд Тейлора в окрестности невозмущенной границы раздела при $z=0$. Учитывая только линейные члены в возмущениях, получаем следующие выражения при $z=0$ (см. Приложение)

${\vartheta }_{sA}-{\vartheta }_{lA}+\left(G_s-G_l\right)Z_A=\mathrm{0,}$ (8)

${\vartheta }_{lA}+\left(\text{Г}{k}_h^2+G_l\right)Z_A=\mathrm{0,}$ (9)

$\begin{array}{c}{k}_s{\beta }_s{\vartheta }_{sA}+\left(c_l{\upsilon }_s-f^{\prime }\right){\vartheta }_{lA}-\\ -\left(\begin{array}{c}\frac{k_s{\upsilon }_s{G}_s}{D_s}+f^{\prime }{G}_l-c_l{\upsilon }_s{G}_l+\\ +c_l\varpi \left({\vartheta }_{\infty }-{\vartheta }_{int}\right)+L_V\varpi \end{array}\right)Z_A=\mathrm{0,}\end{array}$ (10)

где $k_h^2=k_x^2+k_y^2$ и $f^{\prime }=df/d{\vartheta }_l$ при ${\vartheta }_l={\vartheta }_{int},$ то есть

$\begin{array}{c}{f}^{\prime }=2^{4/3}\lambda k_l{\left[\frac{\alpha \left({\vartheta }_{int}\right)g}{D_l\varkappa }\left({\vartheta }_{\infty }-{\vartheta }_{int}\right)\right]}^{1/3}\times \\ \times \left(\frac{5\cdot {10}^{-6}\left({\vartheta }_{\infty }-{\vartheta }_{int}\right)}{\alpha \left({\vartheta }_{int}\right)}-\frac{4}{3}\right).\end{array}$ (11)

Уравнения теплопереноса (1) и (2) приводят к следующим соотношениям при $z=0$

$\left(\varpi -{\beta }_l{\upsilon }_s-\left({\beta }_l^2-k_h^2\right)D_l\right){\vartheta }_{lA}+U_z{G}_l{Z}_A=\mathrm{0,}$ (12)

${\beta }_s=-\frac{{\upsilon }_s}{2D_s}\left(1\pm \sqrt{1+\frac{4D_s}{{\upsilon }_s^2}\left(\varpi +k_h^2{D}_s\right)}\right),$ (13)

где $U_z=-\partial u_{z0}/\partial z$ при $z=0$ $—$ скорость расширения ( $u_{z0}$ $—$ установившаяся $z$ $—$ компонента скорости жидкости).

Объединяя уравнения (9) и (12), получаем ${\beta }_l$ в виде:

${\beta }_l=-\frac{{\upsilon }_s}{2D_l}\left(1\pm \sqrt{1+\frac{4D_l}{{\upsilon }_s^2}\left(\varpi +k_h^2{D}_l-\frac{U_z{G}_l}{\text{Г}{k}_h^2+G_l}\right)}\right).$ (14)

Также важно отметить, что знаки $\pm$ в выражениях для ${\beta }_s$ и ${\beta }_l$ определяют, будут ли возмущения затухать или же усиливаться. К примеру, если возмущения возникают в твердой фазе при $z<0$ на некотором расстоянии от межфазной границы $z=0$, то они будут затухать или усиливаться, в случае отрицательного или положительного значения ${\beta }_s$. Если же возмущения возникают в жидкости вдали от межфазной границы и распространяются в ее сторону, то они будут затухать или усиливаться в зависимости от знака ${\beta }_l$ (отрицательного или положительного). В реальных лабораторных установках или природных процессах это соответствует ограниченной области процесса кристаллизации: твердые стенки расположены на определенных расстояниях $z<0$ в твердом материале и $z>0$ в жидкой фазе.

Из уравнений (8) и (9) получаем дисперсионное соотношение:

$\varpi {\text{П}}_2\left(\varpi ,k_h\right)={\text{П}}_1\left(\varpi ,k_h\right),$ (15)

Где

$\begin{array}{c}{П}_1\left(\varpi ,k_h\right)=f^{\prime }Гk_h^2-c_l{\upsilon }_sГk_h^2-\\ -\frac{k_s{\upsilon }_s{G}_s}{D_s}-k_s{\beta }_s\left(Гk_h^2+G_s\right),\end{array}$ (16)

${П}_2\left(\varpi ,k_h\right)=c_l\left({\vartheta }_{\infty }-{\vartheta }_{int}\right)+L_V.$ (17)

На рис. 2a приведены дисперсионные кривые, рассчитанные в соответствии с соотношением (15) для расплавленного металла, параметры которого приведены в работе [34]. Анализ уравнения (15) показал наличие только вещественных решений. Из графиков видно, что процесс направленной кристаллизации при наличии интенсивной конвекции является морфологически неустойчивым ( $\varpi$ растет с увеличением волнового числа $k_h$ и уменьшением температуры жидкости вдали от фронта кристаллизации ${\vartheta }_{\infty }$, см. рис. 2a и рис. 3a). Такое поведение следует из $\varpi >0$, рассчитанного в широком диапазоне волновых чисел $k_h$. При этом для всех точек, показанных на рис.2a, значение ${\beta }_s$ отрицательно. Физически это означает, что возмущение температуры твердой фазы, возникающее на некотором расстоянии от фронта кристаллизации при $z<0$ (например, флуктуация температуры на охлаждаемой стенке), уменьшается с ростом $z$ и возмущает границу раздела «твердое тело $–$ жидкость» при $z=0$. Также важно отметить, что ${\beta }_s<0$ и ${\beta }_l<0$ описывают затухающие возмущения, распространяющиеся от охлаждаемой границы при $z<0$ в жидкую фазу при $z>0$.

Рис. 2. а $–$ Дисперсионные кривые, построенные в соответствии с выражением (15), $v_s={10}^{-6}$ м/с; б $–$ частота возмущений в зависимости от скорости кристаллизации при $k_h=0$ (динамические возмущения), показанная в соответствии с выражением (15). Для каждой температуры имеется два решения: одно решение устойчивое, при $\varpi =0$, второе $—$ неустойчивое, при $\varpi >0$. Физические параметры взяты для расплавленного металла из работы [34] (табл. 1).

Рис. 3. а $–$ Дисперсионные кривые, построенные в соответствии с выражением (15), $v_s={10}^{-6}$ м/с; б $–$ частота возмущения в зависимости от скорости кристаллизации при $k_h=0$ (динамические возмущения), показанная в соответствии с выражением (15). Для каждой температуры имеется два решения, одно решение устойчивое, при $\varpi =0$, второе $—$ неустойчивое, при $\varpi >0$. Физические параметры взяты для магмы [34] (табл.1).

На рис. 3а приведены дисперсионные кривые, рассчитанные аналогично металлу из уравнения (15), для магмы, данные также были взяты из работы [34] и представлены в табл. 1. На рис. 2a и рис. 3a видно, что поведение кривых подобно для металла и магмы. Значения ${\beta }_s<0$ и ${\beta }_l<0$, как и в случае металла, описывают затухающие возмущения, распространяющиеся от охлаждаемой твердой стенки к жидкости. Другими словами, при кристаллизации двух разных систем частота возмущений $\varpi$ в зависимости от волнового числа $k_h$ имеет подобное поведение.

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

Обнаруженная морфологическая неустойчивость является следствием возмущений, возникающих на охлаждаемой стенке при $z<0$. В случае бесконечной области кристаллизации такого решения не существует из-за требования затухания возмущений при $z\to -\infty$. Следствием обнаруженной неустойчивости является развитие морфологических возмущений в жидкости на границе раздела фаз и появление двухфазной области, заполненной дендритоподобными структурами, как это было показано в экспериментах [34].

Кроме морфологических возмущений с $k_h\ne 0$, в кристаллизующейся системе могут существовать и динамические возмущения с $k_h=0$, представляющие собой возмущения установившейся скорости кристаллизации ${\upsilon }_s$. На рис. 2б показана зависимость частоты возмущений от скорости ${\upsilon }_s$ при двух различных температурах ${\vartheta }_{\infty }$, для металла, а на рис. 3б представлена та же зависимость при двух разных температурах, для магмы. Были обнаружены две ветви решения, одна из которых является динамически неустойчивой ( $\varpi >0$, сплошная линия на рисунках 2б и 3б), а другая $—$ устойчивой ( $\varpi =0$, пунктирная линия на рисунках 2б и 3б). Отметим, что корень $\varpi =0$ из уравнения (15) находится аналитически для динамических возмущений при $k_h=0$. Кроме этого, $Z^{\prime }=Z_{\text{A}}$ и $d{Z}^{\prime }/dt=0$, т.е. система имеет нулевое возмущение скорости и кристаллизуется с установившейся скоростью ${\upsilon }_s$. Какая из этих двух ветвей решения будет реализована, зависит от того, усиливает ли конвекция (создавая неоднородности в полях температуры и концентрации) или ослабляет (выравнивая температуру и концентрацию) динамические возмущения вблизи границы раздела твердое тело/жидкость. Ответ на этот вопрос зависит от характера конвективных течений и расположения кристаллизующегося объекта.

Отметим, что двухфазная область часто возникает в различных геофизических явлениях замерзания льда и застывания магмы, а также в металлургических и химических процессах равновесной и неравновесной кристаллизации расплавов и растворов [36$–$44]. Этим объясняется необходимость разработки теории динамической устойчивости таких процессов под действием конвекции. При возникновении двухфазной области между чисто твердой и жидкой фазами малые температурные возмущения могут вызвать новый сценарий колебательной кристаллизации, когда двухфазный слой динамически неустойчив и колеблется вблизи своей установившейся скорости затвердевания. Такой процесс существенно изменяет распределение примеси в твердом материале и приводит к явлению слоистой ликвации примеси, возникающей в результате динамических колебаний двухфазной области. Для описания этого эффекта в рамках интенсивной конвекции в жидкости далее развита модель с поверхностью разрыва и выполнен анализ ее устойчивости к малым динамическим возмущениям. Эта поверхность разрыва отражает свойства реальной двухфазной области с помощью нового граничного условия. Данное условие представляет собой равенство градиента температуры фазового перехода и градиента температуры жидкой фазы на границе раздела. Кроме этого, это условие определяет отсутствие переохлаждения перед поверхностью разрыва (двухфазной областью) в жидкости. Такой анализ динамической устойчивости проводится ниже в духе ранее разработанных теорий без конвекции [42, 43].

Динамическая неустойчивость затвердевания с двухфазной областью

Рассмотрим узкую квазиравновесную двухфазную область, возникающую перед плоской межфазной границей между твердым телом и жидкостью. Переохлаждение в двухфазной зоне полностью компенсируется скрытой теплотой. Для простоты анализа мы далее заменяем двухфазную область поверхностью разрыва между чисто твердой и чисто жидкой фазами [44, 45] (рис. 1б). Отметим также, что конвекция заторможена в междендритных промежутках, что приводит к накоплению там примеси. Поэтому в математическую модель направленной кристаллизации с двухфазной областью (поверхностью разрыва) требуется включить концентрацию примеси $C_l$ бинарной системы. Таким образом, модель состоит из уравнений теплопроводности (1), (2) и уравнения диффузии примеси в жидкой фазе (диффузией примеси в твердом материале пренебрегаем):

$\begin{array}{c}\frac{\partial C_l}{\partial t}+u\cdot \nabla C_l=D_C{\nabla }^2{C}_l+{\upsilon }_s\frac{\partial C_l}{\partial z}, \\ z>Z^{\prime }\left(t\right),\end{array}$ (18)

где $D_C$ $—$ коэффициент диффузии [46, 47], а $Z^{\prime }\left(t\right)$ $—$ определяет динамические возмущения границы раздела, которая заменяет двухфазную область ( $Z^{\prime }=0$, если рассматривется стационарный сценарий кристаллизации). Поскольку температура фазового перехода зависит от концентрации примеси, граничное условие (3) следует переписать следующим образом:

${\vartheta }_{int}={\vartheta }_l={\vartheta }_s={\vartheta }_{*}+F\left(C_l\right), z=Z^{\prime }\left(t\right),$ (19)

где функция $F\left(C_l\right)$ определяется из фазовой диаграммы ( $F\left(C_l\right)=m{C}_l$ для линейной фазовой диаграммы, $m$ $—$ равновесный наклон линии ликвидус).

Условие теплового баланса (4) выполняется на границе раздела $Z^{\prime }\left(t\right)$. Для замыкания математической модели необходимо записать условие квазиравновесности границы раздела фаз [44, 45]:

$\frac{\partial {\vartheta }_l}{\partial z}=F^{\prime }\left(C_l\right)\frac{\partial C_l}{\partial z}, z=Z^{\prime }\left(t\right),$ (20)

где $F^{\prime }\left(C_l\right)=dF/d{C}_l$ и $F^{\prime }=m$ в случае линейного уравнения линии ликвидус.

В результате математическая модель направленной кристаллизации с двухфазной областью, замененной на поверхность разрыва фаз, состоит из уравнений (1), (2), (18), и граничных условий (4), (19) и (20).

В случае стационарного процесса затвердевания также выполняются выражения (6) и (7) при $Z^{\prime }=0$. Из уравнения диффузии примеси (18) при $z=0$, получаем

$\frac{d^2{C}_{l0}}{d{z}^2}=-\frac{{\upsilon }_s{G}_C}{D_C}, G_C=\frac{d{C}_{l0}}{dz}.$ (21)

Подстрочным индексом "0", как и ранее, обозначены решения, соответствующие установившемуся состоянию.

Двухфазная область может быть динамически неустойчивой, в то время как ее скорость кристаллизации колеблется около квазистационарной скорости ${\upsilon }_s$. Поэтому ниже анализируется динамическая неустойчивость кристаллизации с поверхностью разрыва (рис. 1б). В этом случае возмущения имеют тот же вид с $E=\text{exp}\left(\varpi t\right)$ и $C_l^{\text{'}}=C_{lA}E\text{exp}\left(\beta z\right)$, где $\beta$ $–$ коэффициент нарастания/затухания возмущений для концентрации примеси, а $C_{lA}$ $—$ амплитуда возмущений концентрации примеси.

Подставляя теперь возмущения в граничные условия (4), (19) и (20), получаем четыре уравнения для амплитуд возмущений при $z=0$. Как и прежде, сейчас выполняются уравнения (8) и (10). Два других уравнения имеют следующий вид:

${\vartheta }_{lA}-F^{\prime }{C}_{lA}+\left(G_l-F^{\prime }{G}_C\right)Z_A=\mathrm{0,}$ (22)

$\begin{array}{c}{\beta }_l{\vartheta }_{lA}-\left(\beta F^{\prime }+G_C{F}^{″}\right)C_{lA}+\\ +\left(\frac{{\upsilon }_s{G}_C{F}^{\prime }}{D_C}-\frac{{\upsilon }_s{G}_l}{D_l}-G_C^2{F}^{″}\right)Z_A=\mathrm{0,}\end{array}$ (23)

где $F^{″}=0$ в случае линейного уравнения линии ликвидус.

Уравнения (1), (2) и (18) позволяют выразить ${\beta }_l\left(\varpi \right)$, ${\beta }_s\left(\varpi \right)$ и $\beta \left(\varpi \right)$ следующим образом:

${\beta }_l=-\frac{{\upsilon }_s}{2D_l}\left(1\pm \sqrt{1+\frac{4D_l}{{\upsilon }_s^2}\left(\varpi +\frac{U_z{G}_C}{a}\right)}\right),$ (24)

${\beta }_s=-\frac{{\upsilon }_s}{2D_s}\left(1\pm \sqrt{1+\frac{4D_s\varpi }{{\upsilon }_s^2}}\right),$ (25)

$\beta =-\frac{{\upsilon }_s}{2D_C}\left(1\pm \sqrt{1+\frac{4D_C}{{\upsilon }_s^2}\left(\varpi +\frac{U_z{G}_C}{b}\right)}\right),$ (26)

где

$b\left(\varpi \right)=-\frac{k_s{\beta }_s\left(F^{\prime }{G}_C-G_s\right)+\left(c_l{\upsilon }_s-f^{\prime }\right)G_C{F}^{\prime }-\left(L_V+c_l\left({\vartheta }_{\infty }-{\vartheta }_{int}\right)\right)\varpi -k_s{\upsilon }_s{G}_s/D_s}{k_s{\beta }_s{F}^{\prime }+\left(c_l{\upsilon }_s-f^{\prime }\right)F^{\prime }},$

$a\left(\varpi \right)=F^{\prime }\left[b\left(\varpi \right)+G_C\right]-G_l.$

Как и ранее, знаки ${\beta }_l$, ${\beta }_s$ и $\beta$ определяют направление эволюции возмущений, распространяющихся к границе раздела из твердой или жидкой фаз.

Теперь, исключая амплитуды возмущений из уравнений (8), (10), (22) и (23), получаем следующее выражение для $\varpi$:

$\frac{{\upsilon }_s{G}_C{F}^{\prime }/D_C-{\upsilon }_s{G}_l/D_l-G_C^2{F}^{″}+{\beta }_l\left(F^{\prime }{G}_C-G_l\right)}{\left({\beta }_l-\beta \right)F^{\prime }-G_C{F}^{″}}+b\left(\varpi \right)=0.$ (27)

Знак $\varpi$ (действительной части $\varpi$ ) в уравнении (27) определяет устойчивость/неустойчивость процесса к малым динамическим возмущениям (устойчивость наступает при $\varpi <0$, а неустойчивость $–$ при $\varpi >0$ ).

Анализ выражения (27) для магмы показал наличие только вещественных корней. Они показаны на рис.4 в плоскости $\varpi \left({\upsilon }_s\right)$ для различных температур ${\vartheta }_{\infty }$ (значения расчетных параметров для магмы взяты из работ [34, 46, 47], табл.1).

Рис. 4. Зависимость частоты динамических возмущений от скорости кристаллизации в соответствии с выражением (27) для магмы. Физические параметры магмы приведены в табл. 1 по данным работ [34, 46].

На рисунке отражено, что частота возмущений положительна. Двигаясь слева направо вдоль горизонтальной оси ${\upsilon }_s$, мы видим, что $\varpi$ имеет два решения. Оба решения означают неустойчивость. Первое решение (положительные значения $\varpi$, сплошная линия), означает неустойчивость и увеличивается при увеличении скорости кристаллизации. Второе решение (также положительные значения $\varpi$, пунктирная линия), показывает неустойчивое решение, которое асимптотически приближается к нулю при увеличении скорости кристаллизации. Важным моментом является то, что ${\beta }_s<0$, ${\beta }_l<0$ и $\beta <0$ для всех точек на рис. 4. Как и ранее, это означает, что возмущение, возникающее на охлаждаемой стенке при $z<0$, распространяется через границу раздела фаз (поверхность разрыва) в жидкую фазу и приводит к динамическим возмущениям. В обоих случаях конвекция будет создавать локальные неоднородности полей температуры и концентрации, которые в свою очередь будут приводить к неустойчивости. Аналогичное поведение частоты возмущений было получено для расплавленного металла: были найдены две кривые, означающие неустойчивость. То есть возмущение, исходящее от охлаждаемой стенки, будет усиливаться ( $\varpi >0$, неустойчивость) за счет конвекции. В результате разрывная граница становится динамически неустойчивой и квазистационарный процесс кристаллизации с постоянной скоростью ${\upsilon }_s$ разрушается.

ЗАКЛЮЧЕНИЕ

В настоящей работе было изучено влияние малых морфологических/динамических возмущений на процесс направленной кристаллизации в конечной области с учетом интенсивной конвекции в жидкости. Прежде всего, исследован вопрос об устойчивости плоской границы раздела твердое тело/жидкость, движущейся в установившемся режиме, к малым морфологическим возмущениям поля температуры и скорости кристаллизации. Для этого выведено дисперсионное соотношение, определяющее частоту возмущений, как функцию волнового числа и других параметров системы. Анализ этого соотношения для расплавленного металла и магмы показал, что $\varpi$ положительна в широком диапазоне волновых чисел. Этому решению соответствуют отрицательные значения коэффициентов усиления/затухания возмущений ${\beta }_s$ и ${\beta }_l$. Отрицательные знаки этих коэффициентов свидетельствуют о том, что морфологические возмущения распространяются вдоль направления кристаллизации $z$. А именно, при колебаниях температуры на охлажденной стенке (при $z<0$ ) на некотором расстоянии от границы раздела фаз это возмущение распространяется на жидкую фазу и делает границу раздела твердое тело/жидкость морфологически неустойчивой. Кроме того, показано, что плоская граница раздела твердое тело/жидкость может быть неустойчива к динамическим возмущениям с нулевым волновым числом (возмущениям установившейся скорости кристаллизации ${\upsilon }_s$ ) в широком диапазоне ${\upsilon }_s$. Отметим, что здесь неустойчивая ветвь решений одновременно сосуществует с устойчивой. При этом кристаллизующийся расплав/раствор выбирает одну из них в зависимости от влияния конвекции, которая может либо усиливать, либо ослаблять динамические возмущения, исходящие от твердой фазы (охлаждаемой стенки). Эта морфологическая/динамическая неустойчивость развивается со временем и приводит к образованию двухфазного слоя между чисто твердым материалом и жидкостью.

Для изучения устойчивости такого процесса к динамическим возмущениям (возмущениям установившейся скорости кристаллизации с двухфазной зоной) проведен линейный анализ динамической устойчивости с конвекцией. Результатом этой теории, в которой двухфазная область заменена поверхностью разрыва, является уравнение для частоты динамических возмущений в зависимости от физических параметров системы. Анализ этого уравнения показал, что динамические возмущения могут эволюционировать от твердого тела к жидкому материалу. А именно, если возмущение возникает на охлаждаемой стенке при $z<0$, то оно уменьшается и распространяется в жидкую фазу при $z>0$ при ${\beta }_s<0$, ${\beta }_l<0$ и $\beta <0$. Это возмущение может эволюционировать со временем и приводить к динамической неустойчивости в широком диапазоне скоростей кристаллизации ${\upsilon }_s$. Система имеет две ветви решений (обе неустойчивые), которые сосуществуют одновременно. Конвекция в расплаве приводит к усилению динамических возмущений. В результате поверхность разрыва, имитирующая двухфазную область, оказывается неустойчивой к малым динамическим возмущениям скорости кристаллизации, т.е. устойчивый процесс с двухфазной областью разрушается. В системе образуется более сложный сценарий кристаллизации с неустановившейся скоростью и переменной толщиной двухфазного слоя.

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

Это означает, что ряд теорий устойчивой/неустойчивой кристаллизации с конвекцией (а возможно, и при ее отсутствии) необходимо пересмотреть с учетом ограниченности области кристаллизации. К этому, например, относится теория устойчивого роста вершины дендрита, позволяющая отбирать установившуюся скорость роста кристалла в зависимости от кривизны и переохлаждения [48$–$52]. Другой пример $—$ устойчивость границы раздела лед/океан к морфологическим возмущениям с учетом конечной глубины жидкости [53$–$55].

ФИНАНСИРОВАНИЕ

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

ПРИЛОЖЕНИЕ

В этом разделе мы поясняем математический вывод граничных условий (8)$–$(10), используя в качестве примера граничное условие (8): ${\vartheta }_l={\vartheta }_s$. Итак, разлагая это равенство в ряд Тейлора при $z=0$ и сохраняя только линейные члены, получаем

${\vartheta }_{s0}+{\left(\frac{d{\vartheta }_{s0}}{dz}\right)}_{z=0}{Z}^{\prime }+{\vartheta }_s^{\text{'}}={\vartheta }_{l0}+{\left(\frac{d{\vartheta }_{l0}}{dz}\right)}_{z=0}{Z}^{\prime }+{\vartheta }_l^{\text{'}}, z=0.$

Отметим, что возмущения предполагаются достаточно малыми, чтобы все граничные условия можно было перенести на невозмущенную границу $z=0$. Тогда, имея в виду невозмущенное равенство ${\vartheta }_{s0}={\vartheta }_{l0}$ при $z=0$, приходим к следующему выражению:

${\vartheta }_s^{\text{'}}-{\vartheta }_l^{\text{'}}+\left(G_s-G_l\right)Z^{\prime }=\mathrm{0,} z=0.$

Теперь, чтобы получить уравнение (8), учтем в этом равенстве возмущения ${\vartheta }_s^{\text{'}}$, ${\vartheta }_l^{\text{'}}$ и $Z^{\prime }$. Аналогичным образом можно получить выражения (9) и (10).

Теперь покажем, как получить ${\beta }_l$ из выражения (14). Для этого введем стационарные проекции скоростей жидкости $u_{x0}$, $u_{y0}$, $u_{z0}$ и их морфологические возмущения $u_x^{\text{'}}$, $u_y^{\text{'}}$ и $u_z^{\text{'}}$. Учитывая это, получаем $u_x=u_{x0}+u_x^{\text{'}}$, $u_y=u_{y0}+u_y^{\text{'}}$ и $u_z=u_{z0}+u_z^{\text{'}}$. Теперь, принимая во внимание ${\vartheta }_{l0}={\vartheta }_{l0}\left(z\right)$, получаем

$\begin{array}{l}\frac{\partial {\vartheta }_l^{\text{'}}}{\partial t}-{\upsilon }_s\frac{d{\vartheta }_{l0}}{dz}-{\upsilon }_s\frac{\partial {\vartheta }_l^{\text{'}}}{\partial z}+\left(u_{x0}+u_x^{\text{'}}\right)\frac{\partial {\vartheta }_l^{\text{'}}}{\partial x}+\left(u_{y0}+u_y^{\text{'}}\right)\frac{\partial {\vartheta }_l^{\text{'}}}{\partial y}+\left(u_{z0}+u_z^{\text{'}}\right)\left(\frac{d{\vartheta }_{l0}}{dz}+\frac{\partial {\vartheta }_l^{\text{'}}}{\partial z}\right)=\\ =D_l\left(\frac{d^2{\vartheta }_{l0}}{d{z}^2}+\frac{{\partial }^2{\vartheta }_l^{\text{'}}}{\partial z^2}+\frac{{\partial }^2{\vartheta }_l^{\text{'}}}{\partial x^2}+\frac{{\partial }^2{\vartheta }_l^{\text{'}}}{\partial y^2}\right).\end{array}$

Важными аспектами выполненного анализа являются следующие факторы:

(i) В теории учитываются только линейные члены по морфологическим и динамическим возмущениям.

(ii) Учитывается, что уравнение стационарной теплопроводности имеет вид:

$-{\upsilon }_s\frac{d{\vartheta }_{l0}}{dz}=D_l\frac{d^2{\vartheta }_{l0}}{d{z}^2}.$

(iii) Учитываются граничные условия прилипания жидкости к межфазной границе:

$u_{x0}=u_{y0}=u_{z0}, ïðè z=0.$

(iv) Учитывается линейный вид возмущения скорости жидкости:

$u_z=u_{z0}+\frac{\partial u_{z0}}{\partial z}{Z}^{\prime }+u_z^{\text{'}}=\mathrm{0,}$ при $z=0$ где $u_z^{\text{'}}=U_z{Z}^{\prime },$ при $z=0.$

Отметим, что подстановка в это уравнение возмущений, приводит к выражениям (12) и (14) при $z=0.$