A1  ^0.20 Запишите выражение для $v$ в терминах $T$ и физических констант. Отсюда получите закон пропорциональности $D(T)$.

Как известно из термодинамики, $v=\sqrt{3kT/m}$, где $m$ — масса частицы, а $k$ — постоянная Больцмана. Поскольку характерное время перескока частицы между соседними междоузлиями обратно пропорционально вероятность перехода $\times$ скорость частицы, а коэффициент диффузии обратно пропорционален времени перескока, то $D(T)\propto v\exp(-E_g/kT)\propto\sqrt T\exp(-E_g/kT)$.

A2  ^0.40 Запишите, как связаны между собой среднеквадратичное перемещение частицы $\Delta r_{sq}$ с коэффициентом диффузии $D$ и временем симуляции $t$ в трёхмерном случае.

Примечание: движение частицы вдоль каждой из трёх осей в случайном блуждании можно считать независимым. Коэффициент диффузии $D$ в одномерном случае связан с размером ячейки $a$ и характерным временем перескока $\tau$ как $D=a^2 /2 \tau$.

Примечание: если вы не справились с этим пунктом, в дальнейшем считайте $\Delta r_{sq}=\sqrt{Dt}$.

Найдём сначала средний квадрат смещения частицы вдоль одной из осей (пусть это будет ось $x$). В этом случае движение аналогично рассмотренному в предварительной лекции, а средний квадрат перемещения равен:\[\overline{x^2}=na^2,\]где $n$ — число перескоков. Число перескоков может быть оценено как:\[n=t/\tau.\]Учитывая, что $D=a^2/2\tau$, получаем в итоге:\[\overline{x^2}=2Dt.\]Наконец, поскольку движение вдоль различных осей происходит независимо, $\overline{y^2}=\overline{z^2}=2Dt$. Воспользовавшись теоремой Пифагора, получаем:\[\Delta r_{sq}=\sqrt{\overline{x^2}+\overline{y^2}+\overline{z^2}}=\sqrt{6Dt}\]

A2  ^2.00 Снимите точки для нахождения $D(T)$ и линеаризуйте полученную зависимость. Найдите энергию $E_g$ в условных единицах с точностью хотя бы $2\text%$. Явно укажите, как вы достигаете этой точности.

Проводите измерения в диапазоне $T\in[1;10]$, $t\in[1;10]$ и $N\in[10^2;10^4]$.

Для нахождения коэффициента диффузии используем результат предыдущего пункта. Поскольку $\Delta r_{sq}=6Dt$, то $D=\Delta r_{sq}^2/6t$. Подставляя в результат A1, получим $\Delta r_{sq}^2/4t\sqrt T=\exp(-E_g/kT)$. Самый удобный вариант – найти угловой коэффициент зависимости $Y=\ln\left[\Delta r_{sq}^2 t^{-1}T^{-1/2}\right]$ от $X=1/T$.

Как гарантировать точность? Во-первых, очень желательно иметь $\Delta r_{sq}\gg1$. Во-вторых, $1/T$ меняется на $0.9$, а $E_g$ пропорционально угловому коэффициенту зависимости. Чтобы гарантировать точность $2\%$, нужно, чтобы сумма погрешностей при граничных температурах $\sigma_y(0.1)+\sigma_y(1.0)\le0.02\cdot0.9 E_g$. Т.к. $y=\ln[\Delta r_{sq}^2T^{-1/2}]$, то $\varepsilon_r(0.1)+\varepsilon_r(1.0)\le0.009 E_g$. Можно потребовать от всех точек $\varepsilon_r\le0.0045 E_g$.

Снимем точки и занесём в таблицу:

$t$	10	10	10	10	10	10	10	10	10	10
$T$	1.000	1.111	1.250	1.429	1.667	2.000	2.500	3.333	5.000	10.00
$N$	10000	10000	10000	10000	10000	10000	10000	10000	10000	10000
$\Delta r_{sq}$	13.750	17.009	20.882	25.999	32.154	40.362	50.733	65.932	86.227	122.632
$\sigma_r$	0.053	0.066	0.081	0.102	0.127	0.156	0.198	0.257	0.336	0.477
$X$	1.000	0.900	0.800	0.700	0.6	0.5	0.4	0.3	0.2	0.1
$Y$	2.939	3.312	3.664	4.035	4.383	4.747	5.092	5.473	5.807	6.164

Нетрудно видеть, что требуемая точность нами достигнута. Воспользовавшись InteractiveBuffer.nb, построим график $y(x)$ и найдём его угловой коэффициент. Он будет равен\[-E_g=-3.58\pm0.02\]

B1  ^0.40 Чему будет равен коэффициент $\alpha$ в такой структуре?

Примечание: не обязательно решать этот пункт, чтобы перейти к численному моделированию.

Из соображений симметрии можно предположить, что эффективный коэффициент диффузии будет таким же, как и вдоль любой из осей поликристаллической структуры с кубическими зёрнами. Используя данные в условии утверждения, можем заключить, что плоскости зерновых граней, перпендикулярные потоку, почти не влияют на эффективный коэффициент диффузии, и тогда эффективный коэффициент диффузии:\[D_{eff}\approx2\delta D_{gb}/3\implies \alpha=2/3\]

B2  ^2.00 Моделируя диффузию в поликристалле при $t\in[2;5]$, $D_{gb}/D\in[2;5]$, $\delta\in[0.005;0.05]$, получите экспериментальное значение $\alpha$. Сравните с результатом B1.

Примечание: В результате работы программы в файл Bavg.csv записывается зависимость логарифма средней концентрации $\ln\bar c$ частиц примеси от квадрата глубины проникновения $z^2$. Эту зависимость можно анализировать с помощью программы B.nb.

Насколько понятно из условия, чем больше $D_{gb}/D$ и чем меньше $\delta$, тем лучше результат. Учитывая условие возникновения диффузии кинетического типа A, лучше всего использовать большое время симуляции. Таким образом, будет смотреть на результаты моделирования рядом с точкой $(5,5,0.005)$. Угловой коэффициент графика, получаемого в B.nb, равен $k=-\frac{1}{4D_{eff}t}$, поэтому в условных единицах:\[D_{eff}=\frac1{4|k|t}=1+\alpha\delta(D_{gb}/D-1).\]Значит, линеаризуем зависимость $Y=\left[\frac1{4|k|t}-1\right]$ от $X=\delta(D_{gb}/D-1)$. Чтобы величина $X$ менялась заметно, лучше всего плавно увеличивать $\delta$.

При анализе с помощью B.nb немного уменьшим ползунок Width (до $3/4$ от максимума), чтобы уменьшить влияние возможных краевых эффектов. Сохраним конфигурацию Center и Width для всех измерений, чтобы исключить влияние лёгкой нелинейности графика. Занесём точки в таблицу ниже и построим график $y(x)$ с помощью InteractiveBuffer.nb:

$t$	5	5	5	5	5	5	5	5	5	5
$D_{gb}/D$	5	5	5	5	5	5	5	5	5	5
$\delta$	0.005	0.010	0.015	0.020	0.025	0.030	0.035	0.040	0.045	0.050
$|k|$	0.0491948 0.0493256 0.0492711	0.0485796 0.0486521 0.0486798	0.0478542 0.0479781 0.0479744	0.0472941 0.0474704 0.0472866	0.0466705 0.0465584 0.0466160	0.0461591 0.0462083 0.0458631	0.0454212 0.0458370 0.0448760	0.0445461 0.0446153 0.0448859	0.0444450 0.0443558 0.0439954	0.0438239 0.0431223 0.0432071
$Y$	0.01637 0.01367 0.01479	0.02924 0.02770 0.02712	0.04484 0.04212 0.04222	0.05721 0.05329 0.05738	0.07134 0.07392 0.07259	0.08321 0.08206 0.09020	0.10081 0.09082 0.11418	0.12243 0.12069 0.11394	0.12499 0.12725 0.13648	0.14093 0.15949 0.15722
$X$	0.02	0.04	0.06	0.08	0.10	0.12	0.14	0.16	0.18	0.20

Результаты моделирования дают результат $\alpha=0.75\pm0.09$, что в пределах погрешности совпадает с теоретическим ответом

C1  ^0.40 Что можно сказать к коэффициентах, не прибегая к численному моделированию? Найдите все значения коэффициентов, которые можно определить априори.

Единственное, что работает априори — это метод размерностей. Все величины в правой части, кроме $a$, удобно обезразмерены, а тогда $\psi = -6/5$

C2  ^4.00 Путём численного моделирования найдите остальные коэффициенты. Проводите измерения в диапазоне $t\in[0.005;0.200]$, $D_{gb}/D\in[5;50]$, $\delta\in[0.02;0.20]$.

В результате работы программы в файл Cavg.csv записывается зависимость логарифма средней концентрации $\ln\bar c$ частиц примеси от степени глубины проникновения $z^{6/5}$, а в файл Cmap.csv – зависимость логарифма концентрации $\ln c$ от степени глубины проникновения $z^{6/5}$ и координаты $y$, отсчитываемой вдоль поверхности твёрдого тела. Эти зависимости можно анализировать с помощью программы C.nb.

Учитывая условие, данное в задаче, масштаб диффузии должен быть $\sqrt{2Dt}\ll1$ и $\sqrt{2D_{gb}t}\gg\delta$. Таким образом, если взять, к примеру, начальные параметры $(0.100;25;0.10)$, то каждую из переменных можно менять по отдельности во всём диапазоне, оставаясь в кинетическом режиме B. Для каждой симуляции будем строить объёмный профиль концентрации, чтобы убедиться, что он похож на рисунок в условии. Например, профиль концентрации в $(0.100;25;0.10)$ имеет вид:

Угловой коэффициент $k=\cfrac{\mathrm d\ln\bar c}{\mathrm dz^{6/5}}$ будем искать, выделяя линейный участок на графике, построенном в C.nb. Чтобы получить показатель степени $\xi$, будем менять величину $\delta$ и построим график $\ln|k|(\ln\delta)$ в InteractiveBuffer.nb. Аналогично найдём $\eta$ и $\zeta$. Снимем точки, занесём в таблицу и построим графики:

$\delta$	0.100	0.080	0.060	0.120	0.150
$k$	-1.08567	-1.11796	-1.2865	-0.942494	-0.903676
$\ln\delta$	-2.303	-2.526	-2.813	-2.120	-1.897
$\ln |k|$	0.092	0.112	0.252	-0.059	-0.101
$\delta^{-0.39}[D_{gb}/D]^{-0.68}t^{-0.44}$	0.758	0.826	0.925	0.706	0.647

$D_{gb}/D$	25	35	50	15	10
$k$	-1.08567	-0.865926	-0.675888	-1.55123	-2.00586
$\ln[D_{gb}/D]$	3.219	3.555	3.912	2.708	2.303
$\ln |k|$	0.092	-0.144	-0.392	0.439	0.696
$\delta^{-0.39}[D_{gb}/D]^{-0.68}t^{-0.44}$	0.758	0.603	0.473	1.072	1.413

$t$	0.100	0.080	0.120	0.150	0.200
$k$	-1.08567	-1.18623	-0.99746	-0.908966	-0.799996
$\ln t$	-2.303	-2.526	-2.120	-1.897	-1.609
$\ln |k|$	0.092	0.171	-0.003	-0.095	-0.223
$\delta^{-0.39}[D_{gb}/D]^{-0.68}t^{-0.44}$	0.758	0.836	0.699	0.634	0.558

Искомые показатели степени — это просто соответствующие угловые коэффициенты построенных графиков. Итого:

C3  ^0.60 Найдите как можно точнее величину $C$.

Зная показатели степени, несложно найти $C$ как угловой коэффициент графика, в котором $X=\delta^{-0.39} [G_{gb}/D]^{-0.68} t^{-0.44}$ и $Y=|k|$. Пересчитаем точки в таблицах предыдущего пункта и построим график в InteractiveBuffer.nb:

Угловой коэффициент: