Integral formulation of solutions of boundary-value problems of heat and mass transfer in domains with moving boundaries

Full Text

Введение

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

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

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

Структура материалов, формирующаяся в процессе протекания в них фазовых превращений, зависит от конечного распределения в них центров новой фазы по размерам. Это распределение определяется следующими двумя процессами: зарождением центров новой фазы и процессом их последующего роста. Первый из указанных процессов моделируется начально-краевыми задачами для уравнения типа Фоккера – Планка, второй – краевыми задачами для уравнений теплопроводности или/и диффузии в областях с движущимися границами или моделями фазового поля. Второй из указанных процессов играет не менее важную роль в формировании конечной структуры материала, чем процесс зарождения центров новой фазы, ввиду известной проблемы устойчивости межфазной границы раздела фаз, приводящей к формированию дендритной структуры материала. Имеется огромное число работ, посвященных исследованию устойчивости межфазной границы, разработке аналитических и численных методов решения соответствующих краевых задач [8–18].

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

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

Однако в общем случае процессы тепло- и/или массопереноса, обусловленные фазовым превращением, происходят, как правило, в обеих областях двухфазной области и их необходимо учитывать в математических моделях кинетики фазовых превращений, особенно в тех случаях, когда характерные времена процессов переноса – одного порядка. Насколько нам известно, для случая двухфазных областей решения в интегральной форме получены только в отдельных частных случаях [19, 20].

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

Формулировка модели

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

$\frac{\partial {\Theta }_s\left(M,t\right)}{\partial t}=a_s\Delta {\Theta }_s\left(M,t\right),M\in G_s\left(t\right),t>\mathrm{0,}$ (1)

${\Theta }_s\left(M\mathrm{,0}\right)=\mathrm{0,}M\in G_s\left(0\right),$ (2)

$v_n=k\cdot \left({\Theta }_{eq}\left(M,t\right)-{\Theta }_l\left(M,t\right)\right),M\in S\left(t\right),t>\mathrm{0,}$ (3)

${\Theta }_s\left(M,t\right)={\Theta }_l\left(M,t\right)+\Delta {\Theta }_{ls},M\in S\left(t\right),t>\mathrm{0,}$ (4)

$\frac{\partial {\Theta }_l\left(M,t\right)}{\partial t}=a_l\Delta {\Theta }_l\left(M,t\right),M\in G_l\left(t\right),t>\mathrm{0,}$ (5)

${\Theta }_l\left(M\mathrm{,0}\right)=\mathrm{0,}M\in G_l\left(0\right),$ (6)

$q{v}_n={\lambda }_s\frac{\partial {\Theta }_s\left(M,t\right)}{\partial n}-{\lambda }_l\frac{\partial {\Theta }_l\left(M,t\right)}{\partial n},M\in S\left(t\right),t>\mathrm{0,}$ (7)

$\left|{\Theta }_s\left(M,t\right)\right|<\infty ,$ (8)

$\underset{M\to \infty }{\lim }{\Theta }_l\left(M,t\right)=0.$ (9)

Здесь ${\Theta }_s\left(M,t\right)=T_s\left(M,t\right)-T_{si}$, $T_s\left(M,t\right)$ - температура в точке $M\in G_s\left(t\right)$ области $G_s\left(t\right)$, занятой твердой фазой в момент времени t;

$T_{si}$ - начальная температура твердой фазы;

${\Theta }_l\left(M,t\right)=T_l\left(M,t\right)-T_{li}$, $T_l\left(M,t\right)$ - температура в точке $M\in G_l\left(t\right)$ области $G_l\left(t\right)$, занятой жидкой фазой (расплавом) в момент времени t;

$T_{lн}$ - начальная температура жидкой фазы;

$\Delta {\Theta }_{ls}=T_{li}-T_{si}$; ${\Theta }_{eq}\left(M,t\right)=T_{eq}\left(M,t\right)-T_{li},$$M\in S\left(t\right)$ есть переохлаждение жидкой фазы у межфазной поверхности;

$T_{eq}\left(M,t\right)-$ температура равновесия между жидкой и твердой фазами, зависящая в общем случае от кривизны межфазной границы;

$a_s,a_l-$ температуропроводности твердой и жидкой фаз соответственно;

${\lambda }_s,{\lambda }_l-$ теплопроводности твердой и жидкой фаз соответственно;

$v_n-$ нормальная скорость движения межфазной границы;

$\rho$– плотность твердой фазы;

$q-$удельная теплота плавления на единицу объема;

$S\left(t\right)-$ поверхность межфазной границы.

Начальное условие (2) сформулировано в предположении, что в начальный момент времени твердая фаза отсутствовала.

Уравнение (3) есть граничное условие, определяющее температуру на межфазной поверхности. В общем случае это кинетическое условие, которое формулируется в зависимости от механизма роста кристалла. В данном случае граничное условие (3) записано в предположении нормального механизма роста кристалла [1, 2], который реализуется при достаточно больших переохлаждениях жидкой фазы. Условие (4) есть условие равенства температур твердой и жидкой фаз на межфазной поверхности.

При заданной скорости роста условий (1)–(6), (8), (9) достаточно для однозначного определения соответствующих температурных полей. Однако если скорость роста не задана, то необходимо задавать дополнительное условие для ее определения. Таким условием является условие теплового баланса на межфазной поверхности, которое в общем виде записывается в виде, аналогичном уравнению (7).

Условие теплового баланса существенно усложняет возможность получения аналитического решения соответствующей краевой задачи, так как необходимо согласованным образом найти и температурное поле, и зависимость $S\left(t\right)$, связанные друг с другом. В результате для определения зависимости $S\left(t\right)$ необходимо в общем случае решить интегро-дифференциальное уравнение, получающееся после подстановки выражений температурных полей в условие массового баланса [1]. Приближенные аналитические методы решения самосогласованных краевых задач для уравнений теплопроводности и диффузии представлены в [1].

Метод решения

Решение краевых задач (1)–(9) ищем с помощью преобразования Фурье, полагая, что функции ${\Theta }_s\left(M,t\right)$ и ${\Theta }_l\left(M,t\right)$ вне областей определения $G_s\left(t\right)$ и $G_l\left(t\right)$ в своих краевых задачах соответственно равны нулю, а начало координат находится в области $G_s\left(t\right)$.Тогда изображения этих функций определяются следующими формулами:

 ${\overline{\Theta }}_s\left({\overline{k}}_s,t\right)=\frac{1}{{\left(2\pi \right)}^3}\underset{G_s\left(t\right)}{\iiint }{\Theta }_s\left(M,t\right)\exp \left(i{\overline{k}}_s\cdot \overline{r}\right)dV,$ (10)

${\overline{\Theta }}_l\left({\overline{k}}_l,t\right)=\frac{1}{{\left(2\pi \right)}^3}\underset{G_l\left(t\right)}{\iiint }{\Theta }_l\left(M,t\right)\exp \left(i{\overline{k}}_l\cdot \overline{r}\right)dV.$ (11)

Преобразуем далее уравнения краевых задач (1)–(9) по формулам (10), (11) преобразования Фурье с учетом формул Грина и того, что

$\frac{d}{dt}\underset{G\left(t\right)}{\iiint }f\left(M,t\right)dV=\underset{G\left(t\right)}{\iiint }\frac{\partial }{\partial t}f\left(M,t\right)dV+\underset{G\left(t\right)}{\iiint }div\left(\overline{v}f\left(M,t\right)\right)dV\Rightarrow$,

$\frac{d}{dt}\underset{G\left(t\right)}{\iiint }f\left(M,t\right)dV=\underset{G\left(t\right)}{\iiint }\frac{\partial }{\partial t}f\left(M,t\right)dV+\underset{S\left(t\right)}{\iint }{v}_{n}f\left(M,t\right)dS$

где $\overline{v}=d\overline{r}\left(M,t\right)/dt,M\in S\left(t\right),S\left(t\right)=\partial G\left(t\right)-$ граница области $G\left(t\right)$;

$\overline{r}\left(M,t\right),M\in S\left(t\right)-$радиус-вектор точки М, принадлежащей поверхности $S\left(t\right)$;

$v_n-$ нормальная скорость движения границы области $G\left(t\right)$;

$n-$внешняя нормаль к $S\left(t\right)$.

Получим в результате уравнения

$\frac{d{\overline{\Theta }}_s\left({\overline{k}}_s,t\right)}{dt}+a_s{\overline{k}}_s^2{\overline{\Theta }}_s\left({\overline{k}}_s,t\right)=\frac{1}{{\left(2\pi \right)}^3}\underset{S\left(t\right)}{\iint }{v}_n\left(M,t\right){\Theta }_s\left(M,t\right)\exp \left(i{\overline{k}}_s\cdot {\overline{r}}_M\right)d{S}_M+$

$+\frac{a_s}{{\left(2\pi \right)}^3}\underset{S\left(t\right)}{\iint }\left(\exp \left(i{\overline{k}}_s\cdot {\overline{r}}_M\right)\frac{\partial {\Theta }_s\left(M,t\right)}{\partial n_M}-{\Theta }_s\left(M,t\right)\frac{\partial \exp \left(i{\overline{k}}_s\cdot {\overline{r}}_M\right)}{\partial n_M}\right)d{S}_M=g_s\left(t\right),$ (12)

$\frac{d{\overline{\Theta }}_l\left({\overline{k}}_l,t\right)}{dt}+a_l{\overline{k}}_l^2{\overline{\Theta }}_l\left({\overline{k}}_l,t\right)=-\frac{1}{{\left(2\pi \right)}^3}\underset{S\left(t\right)}{\iint }{v}_n\left(M,t\right){\Theta }_l\left(M,t\right)\exp \left(i{\overline{k}}_l\cdot {\overline{r}}_M\right)d{S}_M-$

$-\frac{a_l}{{\left(2\pi \right)}^3}\underset{S\left(t\right)}{\iint }\left(\exp \left(i{\overline{k}}_l\cdot {\overline{r}}_M\right)\frac{\partial {\Theta }_l\left(M,t\right)}{\partial n_M}-{\Theta }_l\left(M,t\right)\frac{\partial \exp \left(i{\overline{k}}_l\cdot {\overline{r}}_M\right)}{\partial n_M}\right)d{S}_M=g_l\left(t\right).$ (13)

Здесь $n_M-$внешняя нормаль к $S\left(t\right)$ относительно области $G_s\left(t\right)$, $M\in S\left(t\right),$

${\overline{r}}_M-$ радиус-вектор точки М.

Начальные условия для этих задач Коши имеют вид, в рамках предположения о постоянстве начальных температур, ${\overline{\Theta }}_s\left({\overline{k}}_s\mathrm{,0}\right)=0$, ${\overline{\Theta }}_l\left({\overline{k}}_l\mathrm{,0}\right)=0$. Тогда решения уравнений (12), (13) можно представить в виде

${\overline{\Theta }}_s\left({\overline{k}}_s,t\right)=\underset{0}{\overset{t}{\int }}{g}_s\left(\eta \right)\exp \left(-a_s{\overline{k}}_s^2\left(t-\eta \right)\right)d\eta$,${\overline{\Theta }}_l\left({\overline{k}}_l,t\right)=\underset{0}{\overset{t}{\int }}{g}_l\left(\eta \right)\exp \left(-a_l{\overline{k}}_l^2\left(t-\eta \right)\right)d\eta$.

Подставляя в эти формулы соответствующие выражения $g_s\left(t\right)$ и $g_l\left(t\right)$ из формул (12) и (13), и далее ${\overline{\Theta }}_s\left({\overline{k}}_s,t\right)$ и ${\overline{\Theta }}_l\left({\overline{k}}_l,t\right)$ в формулы обращения преобразования Фурье

${\Theta }_s\left(P,t\right)=\int {\overline{\Theta }}_s\left({\overline{k}}_s,t\right)\exp \left(-i{\overline{k}}_s{\overline{r}}_P\right)d{\overline{k}}_s,{\Theta }_l\left(P,t\right)=\int {\overline{\Theta }}_l\left({\overline{k}}_l,t\right)\exp \left(-i{\overline{k}}_l{\overline{r}}_P\right)d{\overline{k}}_l$,

получим следующее интегральное представление решений краевых задач (1)-(9):

$\begin{array}{l}{\Theta }_s\left(P,t\right)=\\ =\frac{1}{{\left(2\pi \right)}^3}\underset{0}{\overset{t}{\int }}d\eta \underset{S\left(\eta \right)}{\iint }\left(a_s\frac{\partial {\Theta }_s\left(M,t\right)}{\partial n_M}+v_n\cdot {\Theta }_s\right)\int \exp \left(-a_s{\overline{k}}_s^2\left(t-\eta \right)-i{\overline{k}}_s\cdot \Delta \overline{r}\right)d{\overline{k}}_{s}d{S}_M-\end{array}$ (14)

$\begin{array}{l}-\frac{a_s}{{\left(2\pi \right)}^3}\underset{0}{\overset{t}{\int }}d\eta \underset{S\left(\eta \right)}{\iint }{\Theta }_s\left(M,\eta \right)\frac{\partial }{\partial n_M}\int \exp \left(-a_s{\overline{k}}_s^2\left(t-\eta \right)-i{\overline{k}}_s\cdot \Delta \overline{r}\right)d{\overline{k}}_{s}d{S}_M,M\in S\left(\eta \right),\\ P\in G_s\left(t\right),\end{array}$

$\begin{array}{l}{\Theta }_l\left(P,t\right)=\\ =-\frac{1}{{\left(2\pi \right)}^3}\underset{0}{\overset{t}{\int }}d\eta \underset{S\left(\eta \right)}{\iint }\left(a_l\frac{\partial {\Theta }_l\left(M,t\right)}{\partial n_M}+v_n\cdot {\Theta }_l\right)\int \exp \left(-a_l{\overline{k}}_l^2\left(t-\eta \right)-i{\overline{k}}_l\cdot \Delta \overline{r}\right)d{\overline{k}}_{l}d{S}_M+\end{array}$(15)

$\begin{array}{l}+\frac{a_l}{{\left(2\pi \right)}^3}\underset{0}{\overset{t}{\int }}d\eta \underset{S\left(\eta \right)}{\iint }{\Theta }_l\left(M,\eta \right)\frac{\partial }{\partial n_M}\int \exp \left(-a_l{\overline{k}}_l^2\left(t-\eta \right)-i{\overline{k}}_l\cdot \Delta \overline{r}\right)d{\overline{k}}_{l}d{S}_M,M\in S\left(\eta \right),\\ P\in G_l\left(t\right),\end{array}$

где $d\overline{k}=d{k}_{x}d{k}_{y}d{k}_z$,$\Delta \overline{r}={\overline{r}}_P-{\overline{r}}_M={\overline{r}}_P-\overline{\rho }\left(M,\eta \right)$,$\overline{\rho }\left(M,\eta \right)\in S\left(\eta \right)$ – радиус-вектор точки M $S\left(\eta \right)$;

${\overline{r}}_P-$радиус-вектор точки Р, а интегрирование по $\overline{k}=\left(k_x,k_y,k_z\right)$ выполняется по каждой компоненте вектора $\overline{k}$ по промежутку $\left(-\infty ,+\infty \right)$.

Выполняя далее в (14), (15) стандартное интегрирование по $\overline{k}$ переходом к сферической системе координат в $\overline{k}$ -пространстве, получим

$\begin{array}{l}\underset{}{\overset{}{\int }}\exp \left(-a{\overline{k}}^2\left(t-\eta \right)-i\overline{k}\cdot \Delta \overline{r}\right)d\overline{k}=\\ =\underset{0}{\overset{\infty }{\int }}\exp \left(-a{\left|\overline{k}\right|}^2\left(t-\eta \right)\right){\left|\overline{k}\right|}^{2}d\left|\overline{k}\right|\underset{0}{\overset{\pi }{\int }}\exp \left(-i\left|\overline{k}\right|\left|\Delta \overline{r}\right|\cos \theta \right)\sin \theta d\theta \underset{0}{\overset{2\pi }{\int }}d\phi =\end{array}$

$=2\pi \underset{0}{\overset{\infty }{\int }}\exp \left(-a{\left|\overline{k}\right|}^2\left(t-\eta \right)\right){\left|\overline{k}\right|}^{2}d\left|\overline{k}\right|\underset{-1}{\overset{1}{\int }}\exp \left(-i\left|\overline{k}\right|\left|\Delta \overline{r}\right|\eta \right)d\eta =$

$=4\pi \underset{0}{\overset{\infty }{\int }}\exp \left(-a{\left|\overline{k}\right|}^2\left(t-\eta \right)\right){\left|\overline{k}\right|}^{2}d\left|\overline{k}\right|\underset{0}{\overset{1}{\int }}\cos \left(\left|\overline{k}\right|\left|\Delta \overline{r}\right|\eta \right)d\eta =$

$=4\pi \underset{0}{\overset{\infty }{\int }}\exp \left(-a{\left|\overline{k}\right|}^2\left(t-\eta \right)\right)\left|\overline{k}\right|\sin \left(\left|\overline{k}\right|\left|\Delta \overline{r}\right|\right)d\left|\overline{k}\right|=\frac{\pi \sqrt{\pi }}{\sqrt{a{\left(t-\eta \right)}^3}}\exp \left(-a\frac{{\left|\Delta \overline{r}\right|}^2}{4a\left(t-\eta \right)}\right)$.

С учетом этого результата выражения (14) и (15) преобразуются к виду

${\Theta }_s\left(P,t\right)=$$\frac{1}{8\sqrt{{\left(\pi a_s\right)}^3}}\times$

$\times \underset{0}{\overset{t}{\int }}\frac{d\eta }{\sqrt{{\left(t-\eta \right)}^3}}\underset{S\left(\eta \right)}{\iint }\left(a_s\frac{\partial {\Theta }_s\left(M,\eta \right)}{\partial n_M}+v_n\cdot {\Theta }_s\left(M,\eta \right)\right)\exp \left(-\frac{{\left({\overline{r}}_P-\overline{\rho }\left(M,\eta \right)\right)}^2}{4a_s\left(t-\eta \right)}\right)d{S}_M-$ (16)

$\begin{array}{l}-\frac{a_s}{8\sqrt{{\left(\pi a_s\right)}^3}}\underset{0}{\overset{t}{\int }}\frac{d\eta }{\sqrt{{\left(t-\eta \right)}^3}}\underset{S\left(\eta \right)}{\iint }{\Theta }_s\left(M,\eta \right)\frac{\partial }{\partial n_M}\exp \left(-\frac{{\left({\overline{r}}_P-\overline{\rho }\left(M,\eta \right)\right)}^2}{4a_s\left(t-\eta \right)}\right)d{S}_M,\\ P\in G_s\left(t\right),\end{array}$

${\Theta }_l\left(P,t\right)=$

$=-\frac{1}{8\sqrt{{\left(\pi a_l\right)}^3}}\underset{0}{\overset{t}{\int }}\frac{d\eta }{\sqrt{{\left(t-\eta \right)}^3}}\underset{S\left(\eta \right)}{\iint }\left(a_l\frac{\partial {\Theta }_l\left(M,\eta \right)}{\partial n_M}+v_n\cdot {\Theta }_l\right)\exp \left(-\frac{{\left({\overline{r}}_P-\overline{\rho }\left(M,\eta \right)\right)}^2}{4a_l\left(t-\eta \right)}\right)d{S}_M+$ (17)

$\begin{array}{l}+\frac{a_l}{8\sqrt{{\left(\pi a_l\right)}^3}}\underset{0}{\overset{t}{\int }}\frac{d\eta }{\sqrt{{\left(t-\eta \right)}^3}}\underset{S\left(\eta \right)}{\iint }{\Theta }_l\left(M,\eta \right)\frac{\partial }{\partial n_M}\exp \left(-\frac{{\left({\overline{r}}_P-\overline{\rho }\left(M,\eta \right)\right)}^2}{4a_l\left(t-\eta \right)}\right)d{S}_M,\\ P\in G_l\left(t\right).\end{array}$

Формулы (16), (17) являются искомыми интегральными представлениями температурных полей в каждой из фаз двухфазной области. Они особенно удобны для расчета температурных полей, когда заданы температура на межфазной границе и нормальная скорость ее движения. В случае процесса роста центра новой фазы, когда температура на межфазной границе и скорость ее роста не заданы, а подлежат самосогласованному определению из условий равенства температур на межфазной границе (4) и теплового баланса (7), мы приходим к системе двух интегро-дифференциальных уравнений, из которой необходимо найти $\partial {\Theta }_s\left(M,t\right)/\partial n$,$\partial {\Theta }_l\left(M,t\right)/\partial n$ и нормальную скорость движения межфазной границы $v_n$.

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

Иллюстрация метода

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

$\begin{array}{l}{\Theta }_l\left(P,t\right)=\\ =\frac{1}{8\sqrt{{\left(\pi a_l\right)}^3}}\underset{0}{\overset{t}{\int }}\frac{d\eta }{\sqrt{{\left(t-\eta \right)}^3}}\underset{S\left(\eta \right)}{\iint }{v}_n\left(\frac{q}{{С}_l{\rho }_l}-{\Theta }_l\left(M,\eta \right)\right)\exp \left(-\frac{{\left({\overline{r}}_P-\overline{\rho }\left(M,\eta \right)\right)}^2}{4a_l\left(t-\eta \right)}\right)d{S}_M+\end{array}$ (18)

$\begin{array}{l}+\frac{a_l}{8\sqrt{{\left(\pi a_l\right)}^3}}\underset{0}{\overset{t}{\int }}\frac{d\eta }{\sqrt{{\left(t-\eta \right)}^3}}\underset{S\left(\eta \right)}{\iint }{\Theta }_l\left(M,\eta \right)\frac{\partial }{\partial n_M}\exp \left(-\frac{{\left({\overline{r}}_P-\overline{\rho }\left(M,\eta \right)\right)}^2}{4a_l\left(t-\eta \right)}\right)d{S}_M,\\ P\in G_l\left(t\right),\end{array}$

где ${С}_l$ – удельная теплоемкость расплава;

${\rho }_l-$ плотность расплава.

В частности, в случае роста сферического кристалла в переохлажденном расплаве из формулы (18), полагая $d{S}_M={\rho }^2\sin \theta d\theta d\phi$ и интегрируя по $\phi$ от 0 до $2\pi$, получим:

${\Theta }_l\left(r,t\right)=$$\frac{2\pi }{8\sqrt{{\left(\pi a_l\right)}^3}}\underset{0}{\overset{t}{\int }}{\rho }^2\left(\eta \right)\frac{d\rho \left(\eta \right)}{d\eta }\left(\frac{q}{{С}_l{\rho }_l}-{\Theta }_l\left(\eta \right)\right)\times$

$\times \exp \left(-\frac{r^2+{\rho }^2\left(\eta \right)}{4a_l\left(t-\eta \right)}\right)\frac{d\eta }{\sqrt{{\left(t-\eta \right)}^3}}\underset{0}{\overset{\pi }{\int }}\exp \left(\frac{r\rho \left(\eta \right)\cos \theta }{2a_l\left(t-\eta \right)}\right)\sin \theta d\theta +$ (19)

$+\frac{2\pi a_l}{8\sqrt{{\left(\pi a_l\right)}^3}}\underset{0}{\overset{t}{\int }}\frac{{\rho }^2\left(\eta \right){\Theta }_l\left(\eta \right)d\eta }{\sqrt{{\left(t-\eta \right)}^3}}\underset{0}{\overset{\pi }{\int }}\frac{\partial }{\partial \rho }\exp \left(-\frac{r^2+{\rho }^2\left(\eta \right)-2r\rho \left(\eta \right)\cos \theta }{4a_l\left(t-\eta \right)}\right)\sin \theta d\theta .$

Здесь $\theta$ есть угол между ${\overline{r}}_P$ и $\overline{\rho }\left(M,\eta \right)$; $r=\left|{\overline{r}}_P\right|,\rho \left(\eta \right)=\left|\overline{\rho }\left(\eta \right)\right|$,${\Theta }_l\left(t\right)$ есть значение ${\Theta }_l\left(r,t\right)$ на межфазной границе, то есть ${\Theta }_l\left(t\right)={\Theta }_l\left(M,t\right),M\in S\left(t\right)$. Проинтегрируем далее по $\theta$:

$\exp \left(-\frac{r^2+{\rho }^2\left(\eta \right)}{4a_l\left(t-\eta \right)}\right)\underset{0}{\overset{\pi }{\int }}\exp \left(\frac{r\rho \left(\eta \right)\cos \theta }{2a_l\left(t-\eta \right)}\right)\sin \theta d\theta =\underset{-1}{\overset{1}{\int }}\exp \left(\frac{r\rho \left(\eta \right)\xi }{2a_l\left(t-\eta \right)}\right)d\xi =$

$=\frac{2a_l\left(t-\eta \right)}{r\rho \left(\eta \right)}\exp \left(-\frac{r^2+{\rho }^2\left(\eta \right)}{4a_l\left(t-\eta \right)}\right)\left(\exp \left(\frac{r\rho \left(\eta \right)}{2a_l\left(t-\eta \right)}\right)-\exp \left(-\frac{r\rho \left(\eta \right)}{2a_l\left(t-\eta \right)}\right)\right)=$

$=\frac{2a_l\left(t-\eta \right)}{r\rho \left(\eta \right)}\left(\exp \left(-\frac{{\left(r-\rho \left(\eta \right)\right)}^2}{4a_l\left(t-\eta \right)}\right)-\exp \left(-\frac{{\left(r+\rho \left(\eta \right)\right)}^2}{4a_l\left(t-\eta \right)}\right)\right),$

$\underset{0}{\overset{\pi }{\int }}\frac{\partial }{\partial \rho }\exp \left(-\frac{r^2+{\rho }^2\left(\eta \right)-2r\rho \left(\eta \right)\cos \theta }{4a_l\left(t-\eta \right)}\right)\sin \theta d\theta =$

$=\exp \left(-\frac{r^2+{\rho }^2\left(\eta \right)}{4a_l\left(t-\eta \right)}\right)\frac{1}{2a_l\left(t-\eta \right)}\underset{0}{\overset{\pi }{\int }}\left(r\cos \theta -\rho \left(\eta \right)\right)\exp \left(\frac{r\rho \left(\eta \right)\cos \theta }{2a_l\left(t-\eta \right)}\right)\sin \theta d\theta =$

$=\exp \left(-\frac{r^2+{\rho }^2\left(\eta \right)}{4a_l\left(t-\eta \right)}\right)\frac{1}{2a_l\left(t-\eta \right)}\underset{-1}{\overset{1}{\int }}r\xi \exp \left(\frac{r\rho \left(\eta \right)\xi }{2a_l\left(t-\eta \right)}\right)d\xi -$

$-\exp \left(-\frac{r^2+{\rho }^2\left(\eta \right)}{4a_l\left(t-\eta \right)}\right)\frac{\rho \left(\eta \right)}{2a_l\left(t-\eta \right)}\underset{-1}{\overset{1}{\int }}\exp \left(\frac{r\rho \left(\eta \right)\xi }{2a_l\left(t-\eta \right)}\right)d\xi =$

$\begin{array}{l}=\exp \left(-\frac{r^2+{\rho }^2\left(\eta \right)}{4a_l\left(t-\eta \right)}\right)\frac{1}{\rho \left(\eta \right)}\times \\ \times \left(\exp \left(\frac{r\rho \left(\eta \right)}{2a_l\left(t-\eta \right)}\right)+\exp \left(-\frac{r\rho \left(\eta \right)}{2a_l\left(t-\eta \right)}\right)-\underset{-1}{\overset{1}{\int }}\exp \left(\frac{r\rho \left(\eta \right)\xi }{2a_l\left(t-\eta \right)}\right)d\xi \right)-\end{array}$

$-\exp \left(-\frac{r^2+{\rho }^2\left(\eta \right)}{4a_l\left(t-\eta \right)}\right)\frac{1}{r}\left(\exp \left(\frac{r\rho \left(\eta \right)}{2a_l\left(t-\eta \right)}\right)-\exp \left(-\frac{r\rho \left(\eta \right)}{2a_l\left(t-\eta \right)}\right)\right)=$

$=\frac{1}{\rho \left(\eta \right)}\left(\exp \left(-\frac{{\left(r-\rho \left(\eta \right)\right)}^2}{4a_l\left(t-\eta \right)}\right)+\exp \left(-\frac{{\left(r+\rho \left(\eta \right)\right)}^2}{4a_l\left(t-\eta \right)}\right)\right)-$

$-\frac{2a_l\left(t-\eta \right)}{r{\rho }^2\left(\eta \right)}\left(\exp \left(-\frac{{\left(r-\rho \left(\eta \right)\right)}^2}{4a_l\left(t-\eta \right)}\right)-\exp \left(-\frac{{\left(r+\rho \left(\eta \right)\right)}^2}{4a_l\left(t-\eta \right)}\right)\right)$

$-\frac{1}{r}\left(\exp \left(-\frac{{\left(r-\rho \left(\eta \right)\right)}^2}{4a_l\left(t-\eta \right)}\right)-\exp \left(-\frac{{\left(r+\rho \left(\eta \right)\right)}^2}{4a_l\left(t-\eta \right)}\right)\right).$

Подставив полученные результаты в (19), получим

${\Theta }_l\left(r,t\right)=$$\frac{1}{2r\sqrt{\pi a_l}}\times$

$\times \underset{0}{\overset{t}{\int }}\rho \left(\eta \right)\frac{d\rho \left(\eta \right)}{d\eta }\left(\frac{q}{{С}_l{\rho }_l}-{\Theta }_l\left(\eta \right)\right)\left(\exp \left(-\frac{{\left(r-\rho \left(\eta \right)\right)}^2}{4a_l\left(t-\eta \right)}\right)-\exp \left(-\frac{{\left(r+\rho \left(\eta \right)\right)}^2}{4a_l\left(t-\eta \right)}\right)\right)\frac{d\eta }{\sqrt{t-\eta }}+$ (20)

$\begin{array}{l}+\frac{1}{4r\sqrt{\pi a_l}}\underset{0}{\overset{t}{\int }}\rho \left(\eta \right){\Theta }_l\left(\eta \right)\left(\left(r-\rho \left(\eta \right)\right)\exp \left(-\frac{{\left(r-\rho \left(\eta \right)\right)}^2}{4a_l\left(t-\eta \right)}\right)+\right.\\ +\left.\left(r+\rho \left(\eta \right)\right)\exp \left(-\frac{{\left(r+\rho \left(\eta \right)\right)}^2}{4a_l\left(t-\eta \right)}\right)\right)\frac{d\eta }{\sqrt{{\left(t-\eta \right)}^3}}-\end{array}$

$-\frac{1}{2r}\sqrt{\frac{a_l}{\pi }}\underset{0}{\overset{t}{\int }}{\Theta }_l\left(\eta \right)\left(\exp \left(-\frac{{\left(r-\rho \left(\eta \right)\right)}^2}{4a_l\left(t-\eta \right)}\right)-\exp \left(-\frac{{\left(r+\rho \left(\eta \right)\right)}^2}{4a_l\left(t-\eta \right)}\right)\right)\frac{d\eta }{\sqrt{t-\eta }}.$            

Выражение (20) для температурного поля является решением следующей краевой задачи, вытекающей из формул (1)–(9) и рассмотренной в [1]:

$\frac{\partial {\Theta }_l\left(r,t\right)}{\partial t}=a_l\left(\frac{{\partial }^2{\Theta }_l\left(r,t\right)}{\partial r^2}+\frac{2}{r}\frac{\partial {\Theta }_l\left(r,t\right)}{\partial r}\right),\rho \left(t\right)<r<\infty ,t>\mathrm{0,}$ (21)

${\Theta }_l\left(\rho \left(t\right),t\right)={\Theta }_{eq}-\frac{1}{k}\frac{d\rho \left(t\right)}{dt},t>\mathrm{0,}$ (22)

$q\frac{d\rho \left(t\right)}{dt}=-{\lambda }_l\frac{\partial {\Theta }_l\left(\rho \left(t\right),t\right)}{\partial r},t>\mathrm{0,}$ (23)

 ${\Theta }_l\left(r\mathrm{,0}\right)=\mathrm{0,}r>\mathrm{0,}$ (24)

 ${\Theta }_l\left(\infty ,t\right)=\mathrm{0,}t\ge 0.$ (25)

Краевая задача (21)–(25) является математической моделью процесса роста сферического кристалла из переохлажденного расплава [1].

В том случае, когда процесс затвердевания расплава контролируется теплопроводностью, температура на межфазной границе близка температуре кристаллизации, то есть переохлаждение на межфазной границе практически снимается и ${\Theta }_l\left(\eta \right)={\Theta }_{eq}=T_{eq}-T_{li}=const$, где $T_{eq}-T_{lн}$ есть исходное переохлаждение расплава, из формулы (19) получается известное автомодельное решение с законом движения межфазной границы $\rho \left(t\right)=2\beta \sqrt{a_{l}t}$. Для этого необходимо найти $\partial {\Theta }_l\left(\rho \left(t\right),t\right)$ исходя из выражения (20) и полагая затем $\rho \left(t\right)=2\beta \sqrt{a_{l}t}$, далее вычислить достаточно сложные интегралы и подставить результат в уравнение (24). В результате придем к следующему уравнению для параметра β, определяющего кинетику роста сферического кристалла из переохлажденного расплава (выкладки мы не приводим ввиду их громоздкости):

$2{\beta }^2\left(1-\beta \sqrt{\pi }{e}^{{\beta }^2}erfc\beta \right)=\frac{{\Theta }_{eq}{C}_l}{q{\rho }_l}$. (26)

Здесь $erfc\beta =\frac{2}{\sqrt{\pi }}\underset{\beta }{\overset{\infty }{\int }}{e}^{-x^2}dx$.

Значения левой части уравнения (26) заключены между нулем и единицей [1]. Следовательно, уравнение (26) имеет решение только при условии $0\le \left({\Theta }_{eq}{C}_l\right)/\left(q{\rho }_l\right)<1$, которое, таким образом, является и условием существования автомодельного решения.

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

$\frac{\partial {\Theta }_s\left(x,t\right)}{\partial t}=a_s\frac{{\partial }^2{\Theta }_s\left(x,t\right)}{\partial x^2},-\infty <x<\rho \left(t\right),t>\mathrm{0,}$ (27)

${\Theta }_s\left(x\mathrm{,0}\right)=\mathrm{0,}-\infty <x<\mathrm{0,}$ (28)

 ${\Theta }_s\left(\rho \left(t\right),t\right)={\Theta }_{sk},t>\mathrm{0,}$ (29)

${\Theta }_s\left(-\infty ,t\right)=\mathrm{0,}t>\mathrm{0,}$ (30)

$\frac{\partial {\Theta }_l\left(x,t\right)}{\partial t}=a_l\frac{{\partial }^2{\Theta }_l\left(x,t\right)}{\partial x^2},\rho \left(t\right)<x<\infty ,t>\mathrm{0,}$ (31)

 ${\Theta }_l\left(x\mathrm{,0}\right)=\mathrm{0,0}<x<\infty ,$ (32)

 ${\Theta }_l\left(\rho \left(t\right),t\right)={\Theta }_{lk},t>0$, (33)

 ${\Theta }_l\left(\infty ,t\right)=\mathrm{0,}t>0$, (34)

 $q\frac{d\rho \left(t\right)}{dt}={\lambda }_s\frac{\partial {\Theta }_s\left(\rho \left(t\right),t\right)}{\partial x}-{\lambda }_l\frac{\partial {\Theta }_l\left(\rho \left(t\right),t\right)}{\partial x},t>0$. (35)

Здесь ${\Theta }_s\left(x,t\right)=T_s\left(x,t\right)-T_{si}$,$T_s\left(x,t\right)-$ температура в твердой фазе в точке x в момент времени t, $T_{si}-$ начальная температура в твердой фазе;${\Theta }_l\left(x,t\right)=T_l\left(x,t\right)-T_{li}$ ,$T_l\left(x,t\right)-$ температура в жидкой фазе в точке x в момент времени t, $T_{li}-$начальная температура в жидкой фазе;${\Theta }_{sk}=T_{eq}-T_{si}$ ,${\Theta }_{lk}=T_{eq}-T_{li}$.

Запишем решения краевых задач исходя из формул (16), (17) и предполагая, что температура на межфазной границе $x=\rho \left(t\right)$ не зависит от координат y и z. Тогда, интегрируя в (16), (17) по y и z, получим

${\Theta }_s\left(x,t\right)=$

$=\frac{1}{2\sqrt{\pi a_s}}\underset{0}{\overset{t}{\int }}\left(a_s\frac{\partial {\Theta }_s\left(\rho \left(\eta \right),\eta \right)}{\partial x}+\frac{d\rho \left(\eta \right)}{d\eta }{\Theta }_s\left(\rho \left(\eta \right),\eta \right)\right)\exp \left(-\frac{{\left(x-\rho \left(\eta \right)\right)}^2}{4a_s\left(t-\eta \right)}\right)\frac{d\eta }{\sqrt{t-\eta }}-$ (36)

$-\frac{1}{4\sqrt{\pi a_s}}\underset{0}{\overset{t}{\int }}\left(x-\rho \left(\eta \right)\right)\exp \left(-\frac{{\left(x-\rho \left(\eta \right)\right)}^2}{4a_s\left(t-\eta \right)}\right){\Theta }_s\left(\rho \left(\eta \right),\eta \right)\frac{d\eta }{\sqrt{{\left(t-\eta \right)}^3}},$

${\Theta }_l\left(x,t\right)=$

$=-\frac{1}{2\sqrt{\pi a_l}}\underset{0}{\overset{t}{\int }}\left(a_l\frac{\partial {\Theta }_l\left(\rho \left(\eta \right),\eta \right)}{\partial x}+\frac{d\rho \left(\eta \right)}{d\eta }{\Theta }_l\left(\rho \left(\eta \right),\eta \right)\right)\exp \left(-\frac{{\left(x-\rho \left(\eta \right)\right)}^2}{4a_l\left(t-\eta \right)}\right)\frac{d\eta }{\sqrt{t-\eta }}+$ (37)

$+\frac{1}{4\sqrt{\pi a_l}}\underset{0}{\overset{t}{\int }}\left(x-\rho \left(\eta \right)\right)\exp \left(-\frac{{\left(x-\rho \left(\eta \right)\right)}^2}{4a_l\left(t-\eta \right)}\right){\Theta }_l\left(\rho \left(\eta \right),\eta \right)\frac{d\eta }{\sqrt{{\left(t-\eta \right)}^3}}.$             

Рассмотрим автомодельный случай, когда процесс затвердевания контролируется теплопроводностью, а не кинетическими процессами на межфазной границе. В этом случае $\rho \left(t\right)=2\beta \sqrt{t}$, ${\Theta }_s\left(\rho \left(t\right),t\right)={\Theta }_{sk},{\Theta }_l\left(\rho \left(t\right),t\right)={\Theta }_{lk}$.

Тогда, так как $\frac{d\rho \left(t\right)}{dt}=\frac{\beta }{\sqrt{t}}$, то

$\frac{\partial {\Theta }_s\left(\rho \left(t\right),t\right)}{\partial x}=\frac{{\beta }_s}{\sqrt{t}}{\Theta }_{sk}$, (38)

$\frac{\partial {\Theta }_l\left(\rho \left(t\right),t\right)}{\partial x}=\frac{{\beta }_l}{\sqrt{t}}{\Theta }_{lk}$. (39)

Подставив (38), (39) в уравнение (33), получим

 $q\beta ={\lambda }_s{\beta }_s{\Theta }_{sk}-{\lambda }_l{\beta }_l{\Theta }_{lk}.$ (40)

Уравнение (40) содержит три неизвестных величины $\beta ,{\beta }_s,{\beta }_l$. Для того чтобы найти параметр $\beta$, определяющий кинетику процесса кристаллизации, необходимо определить величины ${\beta }_s,{\beta }_l$. Их мы найдем, удовлетворяя уравнения (38), (39) с помощью полученных решений (36), (37). Однако прежде чем найти $\frac{\partial {\Theta }_s\left(\rho \left(t\right),t\right)}{\partial x},\frac{\partial {\Theta }_l\left(\rho \left(t\right),t\right)}{\partial x}$, преобразуем выражения (36), (37) с учетом того, что $d\frac{\left(x-\rho \left(\eta \right)\right)}{2a\sqrt{t-\eta }}=\left(\frac{\left(x-\rho \left(\eta \right)\right)}{4a\sqrt{{\left(t-\eta \right)}^3}}-\frac{1}{2a\sqrt{t-\eta }}\frac{d\rho \left(\eta \right)}{d\eta }\right)d\eta$.

Это позволяет объединить второе и третье слагаемые в (36), (37). В результате такого объединения формулы (36), (37) можно привести с учетом того, что ${\Theta }_s\left(\rho \left(t\right),t\right)={\Theta }_{sk},{\Theta }_l\left(\rho \left(t\right),t\right)={\Theta }_{lk}$, к следующему виду:

$\begin{array}{l}{\Theta }_s\left(x,t\right)=-\frac{1}{2}{\Theta }_{sk}erfc\frac{x}{2\sqrt{a_{s}t}}+\\ +\frac{a_s}{2\sqrt{\pi a_s}}\underset{0}{\overset{t}{\int }}\frac{\partial {\Theta }_s\left(\rho \left(\eta \right),\eta \right)}{\partial x}\exp \left(-\frac{{\left(x-\rho \left(\eta \right)\right)}^2}{4a_s\left(t-\eta \right)}\right)\frac{d\eta }{\sqrt{t-\eta }},\end{array}$ (41)

 $\begin{array}{l}{\Theta }_l\left(x,t\right)=\frac{1}{2}{\Theta }_{lk}erfc\frac{x}{2\sqrt{a_{s}t}}-\\ -\frac{a_l}{2\sqrt{\pi a_l}}\underset{0}{\overset{t}{\int }}\frac{\partial {\Theta }_l\left(\rho \left(\eta \right),\eta \right)}{\partial x}\exp \left(-\frac{{\left(x-\rho \left(\eta \right)\right)}^2}{4a_l\left(t-\eta \right)}\right)\frac{d\eta }{\sqrt{t-\eta }}.\end{array}$ (42)

Тогда

$\frac{\partial {\Theta }_s\left(\rho \left(t\right),t\right)}{\partial x}=$

$\begin{array}{l}={\Theta }_{sk}\exp \left(-\frac{{\beta }^2}{a_s}\right)\frac{1}{2\sqrt{\pi a_{s}t}}+\frac{{\beta }_s{\Theta }_{sk}}{2\sqrt{t}}-\\ -\frac{\beta {\beta }_s}{2\sqrt{\pi a_s}}{\Theta }_{sk}\underset{0}{\overset{t}{\int }}\frac{1}{\sqrt{\eta }}\frac{1}{\sqrt{t}+\sqrt{\eta }}\exp \left(-\frac{{\beta }^2\left(\sqrt{t}-\sqrt{\eta }\right)}{a_s\left(\sqrt{t}+\sqrt{\eta }\right)}\right)\frac{d\eta }{\sqrt{t-\eta }}\end{array}$ (43)

 $\frac{\partial {\Theta }_l\left(\rho \left(t\right),t\right)}{\partial x}=$

$\begin{array}{l}=-{\Theta }_{lk}\exp \left(-\frac{{\beta }^2}{a_l}\right)\frac{1}{2\sqrt{\pi a_{l}t}}+\\ +\frac{{\beta }_l{\Theta }_{lk}}{2\sqrt{t}}+\frac{\beta {\beta }_l}{2\sqrt{\pi a_l}}{\Theta }_{lk}\underset{0}{\overset{t}{\int }}\frac{1}{\sqrt{\eta }}\frac{1}{\sqrt{t}+\sqrt{\eta }}\exp \left(-\frac{{\beta }^2\left(\sqrt{t}-\sqrt{\eta }\right)}{a_l\left(\sqrt{t}+\sqrt{\eta }\right)}\right)\frac{d\eta }{\sqrt{t-\eta }}.\end{array}$ (44)

В выражения (43), (44) входят однотипные интегралы, которые также появляются при нахождении автомодельного решения рассмотренного выше однофазного случая. Интегралы такого типа подробно рассмотрены в [1] и не приводятся в известных таблицах определенных интегралов. Поэтому рассмотрим метод их вычисления. Для этого ввиду их однотипности ограничимся рассмотрением интегралов

$S=\underset{0}{\overset{t}{\int }}\frac{1}{\sqrt{\eta }}\frac{1}{\sqrt{t}+\sqrt{\eta }}\exp \left(-\frac{{\gamma }^2\left(\sqrt{t}-\sqrt{\eta }\right)}{\left(\sqrt{t}+\sqrt{\eta }\right)}\right)\frac{d\eta }{\sqrt{t-\eta }}$

и для его вычисления вводим переменную интегрирования $\xi =\eta /t$, после чего получим

$S=\frac{1}{\sqrt{t}}\underset{0}{\overset{1}{\int }}\frac{1}{\sqrt{\xi }}\frac{1}{1+\sqrt{\xi }}\exp \left(-\frac{{\gamma }^2\left(1-\sqrt{\xi }\right)}{\left(1+\sqrt{\xi }\right)}\right)\frac{d\xi }{\sqrt{1-\xi }}$.

Полученный интеграл вычисляется подстановкой $\frac{1-\sqrt{\xi }}{1+\sqrt{\xi }}=y^2$.

Тогда, учитывая, что

$\xi ={\left(\frac{1-y^2}{1+y^2}\right)}^2\mathrm{,1}-\sqrt{\xi }=\frac{2y^2}{1+y^2},d\xi =-8\frac{y\left(1-y^2\right)}{{\left(1+y^2\right)}^3}dy\mathrm{,1}-\xi =\frac{4y^2}{{\left(1+y^2\right)}^2}$,

получим после соответствующего преобразования подынтегрального выражения

$S=\frac{1}{\sqrt{t}}\underset{0}{\overset{1}{\int }}\exp \left(-{\gamma }^2{y}^2\right)dy=\frac{1}{\gamma \sqrt{t}}2\sqrt{\pi }erf\gamma$,

где $erf\gamma =\frac{2}{\sqrt{\pi }}\underset{0}{\overset{\gamma }{\int }}{e}^{-x^2}dx$.

Учтем полученный результат в формулах (43), (44). Получим с учетом (38), (39) соответственно

${\beta }_s=\exp \left(-\frac{{\beta }^2}{a_s}\right)\frac{1}{2\sqrt{\pi a_s}}+\frac{{\beta }_s}{2}-\frac{{\beta }_s}{2}erf\frac{\beta }{\sqrt{a_s}},$ (45)

${\beta }_l=-\exp \left(-\frac{{\beta }^2}{a_l}\right)\frac{1}{2\sqrt{\pi a_l}}+\frac{{\beta }_l}{2}+\frac{{\beta }_l}{2}erf\frac{\beta }{\sqrt{a_l}}.$ (46)

Из уравнений (45), (46) получим следующие выражения для искомых величин ${\beta }_s,{\beta }_l$:

${\beta }_s=\frac{\exp \left(-\frac{{\beta }^2}{a_s}\right)\frac{1}{\sqrt{\pi a_s}}}{1+erf\frac{\beta }{\sqrt{a_s}}},$ (47)

${\beta }_l=-\frac{\exp \left(-\frac{{\beta }^2}{a_l}\right)\frac{1}{\sqrt{\pi a_l}}}{1-erf\frac{\beta }{\sqrt{a_l}}}.$ (48)

Подставив выражения (47), (48) в (40), получим искомое уравнение для параметра $\beta$, определяющее кинетику процесса кристаллизации:

$q\beta =\frac{{\lambda }_s{\Theta }_{sk}\exp \left(-\frac{{\beta }^2}{a_s}\right)\frac{1}{\sqrt{\pi a_s}}}{1+erf\frac{\beta }{\sqrt{a_s}}}+\frac{{\lambda }_l{\Theta }_{lk}\exp \left(-\frac{{\beta }^2}{a_l}\right)\frac{1}{\sqrt{\pi a_l}}}{1-erf\frac{\beta }{\sqrt{a_l}}}.$ (49)

Уравнение (49) позволяет сделать некоторые выводы относительно кинетики кристаллизации расплава на подложке. В частности, из уравнения (49) следует, что при ${\Theta }_{sk}<\mathrm{0,}{\Theta }_{lk}<0$, то есть когда температура кристаллизации (плавления) $T_{eq}$ меньше и начальной температуры подложки $T_{si}$, и начальной температуры расплава $T_{li}$, кристаллизация невозможна, а происходит, наоборот, процесс плавления подложки.

При условии же, когда ${\Theta }_{sk}>\mathrm{0,}{\Theta }_{lk}>0$, то есть когда температура кристаллизации $T_{eq}$ выше и начальной температуры подложки $T_{si}$, и начальной температуры расплава $T_{li}$, процесс кристаллизации гарантированно имеет место и происходит с тем большей скоростью, чем больше указанные начальные переохлаждения подложки ${\Theta }_{sk}$ и расплава ${\Theta }_{lk}$. В этом случае теплоотвод от межфазной границы происходит и в подложку, и в переохлажденный расплав.

Рассмотрим случай, когда ${\Theta }_{sk}<\mathrm{0,}{\Theta }_{lk}>0$, то есть когда начальная температура подложки $T_{si}$ больше температуры кристаллизации $T_{eq}$, а начальная температура расплава меньше температуры кристаллизации $T_{eq}$. Тогда при достаточно малых $\beta$ уравнение (49) примет вид

$\beta =\frac{1}{q\sqrt{\pi }}\left({\Theta }_{sk}\sqrt{{\lambda }_s{C}_s{\rho }_s}+{\Theta }_{lk}\sqrt{{\lambda }_l{C}_l{\rho }_l}\right)$. (50)

Здесь $C_s,{\rho }_s,C_l,{\rho }_l$ – соответственно удельная теплоемкость и плотность твердой и жидкой фаз.

Из уравнения (50) следует, что расплав кристаллизуется на подложке при условии , то есть при условии

$\left|{\Theta }_{sk}\right|\sqrt{{\lambda }_s{C}_s{\rho }_s}<{\Theta }_{lk}\sqrt{{\lambda }_l{C}_l{\rho }_l}$.

Это условие означает, что расплав кристаллизуется на подложке в результате теплоотвода через расплав. При условии

$\left|{\Theta }_{sk}\right|\sqrt{{\lambda }_s{C}_s{\rho }_s}>{\Theta }_{lk}\sqrt{{\lambda }_l{C}_l{\rho }_l}$,

наоборот, плавится подложка, а в случае, когда

$\left|{\Theta }_{sk}\right|\sqrt{{\lambda }_s{C}_s{\rho }_s}={\Theta }_{lk}\sqrt{{\lambda }_l{C}_l{\rho }_l}$,

система расплав – подложка находится в состоянии кинетического равновесия, так как $\beta =0$.

В противоположном случае, когда ${\Theta }_{sk}>\mathrm{0,}{\Theta }_{lk}<0$, расплав кристаллизуется на подложке в результате теплоотвода через подложку при условии $\beta >0$, то есть при условии

${\Theta }_{sk}\sqrt{{\lambda }_s{C}_s{\rho }_s}>\left|{\Theta }_{lk}\right|\sqrt{{\lambda }_l{C}_l{\rho }_l}$.

При условии

${\Theta }_{sk}\sqrt{{\lambda }_s{C}_s{\rho }_s}<\left|{\Theta }_{lk}\right|\sqrt{{\lambda }_l{C}_l{\rho }_l}$,

наоборот, плавится подложка.

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

Выводы

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

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

Кроме того, метод позволяет также получить интегральные представления решений задач, моделирующих кристаллизацию слитков, то есть кристаллизацию в больших объемах [1, 2], что представляет отдельный технологический интерес.

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

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

About the authors

Valentin V. Shevelev

MIREA – Russian technological University, Institute of fine chemical technologies. M.V. Lomonosov; National research technological University "MISIS"

Author for correspondence.
Email: valeshevelev@yandex.ru

Professor

Russian Federation, 78, pr-ct Vernadsky, Moscow, 119454; 4, Leninsky pr-ct, Moscow, 119049

References