Редактирование: Распределение тепла в одномерном кристалле

Выполнил: [[Цветков_Денис|Цветков Д. В.]]

Научный руководитель: [[Антон_Кривцов|Кривцов А. М.]]

== Введение ==

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

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

Пример: сосуд, разделенный перегородкой (рис.1а). В левой части сосуда содержится небольшое количество молекул газа, правая часть сосуда пуста. Убираем перегородку (рис.1б). Постепенно молекулы газа занимают все пространство сосуда (рис.1в). Но, так как молекул очень мало, через некоторое время они снова одновременно окажутся в левой части сосуда (рис.1г).

<gallery widths=253px heights=133px perrow=2 caption="Сосуд с небольшим количеством молекул газа">
Файл: Tcvet_masters_B1.png|Рис.1а)
Файл: Tcvet_masters_B2.png|Рис.1б)
Файл: Tcvet_masters_B3.png|Рис.1в)
Файл: Tcvet_masters_B4.png|Рис.1г)
</gallery>

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

<gallery widths=253px heights=133px perrow=2 caption="Сосуд с большим количеством молекул газа">
Файл: Tcvet_masters_B5.png|Рис.2а)
Файл: Tcvet_masters_B6.png|Рис.2б)
</gallery>


Важно исследовать граничное время, в которое возврат системы становится практически невозможен.

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

== Сравнительный анализ решений задач теплопроводности для бесконечного слоя, находящегося под воздействием ультракороткого лазерного импульса==
=== Метод исследования===
Рассматривается задача одномерного теплоизолированного кристалла, нагретого ультракоротким лазерным импульсом. Производится сравнение следующих решений:


* стохастическое решение --- задача динамики кристалла, температура кристалла задается с помощью случайных скоростей
* численное решение дифференциального уравнения, полученного в <ref name="krivtsov_1dcrystal_2014"/>:
:<math>
\ddot T + \frac 1t \dot T = c^2 T''
</math>
* аналитическое решение, полученное в <ref name="babenkov_tcvetkov"/>:
:<math>
T(t, x) = \frac1\pi \int\limits_{-t}^t \frac{e^{-\gamma \left|x - ct \right| }I_0}{k\sqrt{t^2 - \tau^2}} \,d\tau,
</math>
где <math>I_0</math> --- интенсивность входящего пучка света, <math>\gamma</math> --- коэффициент поглощения света.

Для численного расчета и сравнения задач написано специальное программное обеспечение на языке MATLAB.

=== Постановка стохастической задачи ===

Рассмотрим одномерный кристалл: цепочку одинаковых частиц массы <math>m</math>, соединенных одинаковыми линейными пружинами с жесткостью <math>C</math>:
[[Файл: Tcvet masters Particles.png|Одномерный кристалл]]

Уравнения динамики кристалла имеют вид:
:<math>
\ddot u_n = \omega_0^2 (u_{n-1} - 2u_n + u_{n+1}),\quad \omega_0 = \sqrt{\frac{C}{m}},
</math>
где <math>u_n</math> --- перемещение <math>n</math>-й частицы; <math>n</math> --- индекс, принимающий произвольные целые значения, <math>C</math> --- жесткость связи между частицами, <math>m</math> --- масса частицы.

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

=== Начальные условия===
Требуемое распределение температуры по координате в кристалле на начальный момент времени описывается формулой:
:<math>
T_0(x) = \frac{I_0}k e^{-\gamma x}
</math>
где <math>I_0</math> --- интенсивность входящего пучка света, <math>\gamma</math> --- коэффициент поглощения света, <math>k</math> --- постоянная Больцмана.
Согласно закону Бугера <ref name="Putilov1963"/>, данные начальные условия можно интерпретировать как результат воздействия на кристалл ультракороткого лазерного импульса.



==== Задание начальных скоростей ====
Температура в кристалле задается с помощью скоростей.
Воспользуемся определением кинетической температуры для одномерного случая:
:<math>
kT(x) = m\overline{v_x^2},
</math>
где <math>k</math> --- постоянная Больцмана, <math>\overline x = f(x)</math> --- оператор усреднения по координате. Выразим <math>\overline{v_x^2}</math>:
:<math>
\overline{v_x^2} = \frac{kT(x)}m
</math>
Для того, чтобы получить требуемое распределение, будем задавать скорости с помощью случайного числа <math>\rho</math> в диапазоне <math>[-1; 1]</math>.

Посчитаем, чему равно <math>\overline{\rho^2}</math>:
:<math>
\overline{\rho^2} = \frac 1 2 \int\limits_a^b \rho^2 \,d\rho =  \frac 1 2 \left.  \frac {\rho^3}3 \right|_{-1}^1 = \frac 1 3
</math>
Значит <math>3\overline{\rho^2} = 1</math>.

Заменим <math>T(x)</math> на <math>T_0(x)</math>, генерируемое с помощью <math>\overline{\rho^2}</math>:
:<math>
\overline{v_x^2} = \frac{3k\overline{\rho^2}T_0(x)}m
</math>
Выразим <math>v_x</math>:
:<math>
v_x^2 = 3\rho^2 T_0(x)
</math>
:<math>
v_x = \rho\sqrt{3T_0(x)}
</math>

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

[[Файл: Tcvet_masters_Theme1_IC.png|600px|Начальные условия]]

==== Задание начальных перемещений ====
Если задать нулевые начальные перемещения, то потенциальная энергия кристалла будет равна нулю, и с течением времени примерно половина кинетической энергии преобразуется в потенциальную <ref name="krivtsov_1dcrystal_2014"/>. Чтобы этого избежать, начальные перемещения задаются с помощью деформаций. Деформации задаются таким образом, чтобы кинетическая и потенциальная энергии кристалла в начальный момент времени были равны.

Сначала задаются <math>N-1</math> случайных чисел, сумма которых равна нулю:
:<math>
\sum_{i=0}^{N-1}{\rho_i} = 0,
</math>
где <math>N</math> --- количество частиц в цепочке.
После рассчитывается параметр <math>\varPi_{\rho}</math>:
:<math>
\varPi_{\rho} = C\sum_{i=0}^{N-1}{\rho^2_i}
</math>
и кинетическая энергия <math>K</math>:
:<math>
K = m\sum_{i=0}^{N}{v^2_i}
</math>
Далее, с помощью отношения <math>K</math> к <math>\varPi_{\rho}</math> находятся деформации:
:<math>
\varepsilon_i = \rho_i\sqrt{\frac K\varPi_{\rho}}
</math>

=== Граничные условия ===
Граничные условия --- зеркальные. Это значит, что значения перемещения и скорости частицы на границе равны значениям перемещения и скорости соседней частицы:
:<math>
    u_1 = u_2,  \quad u_N = u_{N-1},   \\
    v_1 = v_2,  \quad v_N = v_{N-1},
</math>
где <math>N</math> --- число частиц в цепочке.


[[Файл: Tcvet masters Theme1 BC.png|600px|Граничные условия]]

При данных условиях энергия не убывает и не прибывает --- на границах стержня выполняется условие теплоизоляции.

=== Результаты ===

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

[[Файл: Tcvet masters Theme1 result.png|600px|Результаты]]

== Задача о переходе энергии механических колебаний в тепло для нелинейного одномерного кристалла ==
=== Постановка задачи ===
Рассмотрим одномерный кристалл: цепочку одинаковых частиц массы <math>m</math>, соединенных одинаковыми нелинейными пружинами с жесткостью <math>C</math>.
Уравнения динамики кристалла имеют вид:

:<math>
\ddot u_k = \omega_0^2 (u_{k-1} - 2u_k + u_{k+1})(1 + u_{k+1} - u_{k-1}),\quad \omega_0 = \sqrt{\frac{C}{m}},
</math>
где <math>u_k</math> --- перемещение <math>k</math>-й частицы; <math>k</math> --- индекс, принимающий произвольные целые значения.
Будем считать, что выполнены условия периодичности: <math>u_{k+N} = u_k</math>, где <math>N \gg 1</math> --- число независимых частиц.

Рассматриваются два случая:
* Детерминированная задача --- при <math>t = 0</math> перемещения частиц равны нулю, а скорости распределены по синусоидальному закону.
* Стохастическая задача --- к скоростям в детерминированной задаче добавляются случайные флуктуации.

=== Детерминированная задача ===
==== Начальные условия ====
Начальные скорости в кристалле задаются следующим образом (рис. \ref{Theme3_IC}):
:<math>
\left. v(x) \right |_{t = 0} = A \sin \( \frac{2\pi x}{L} \),
</math>
где <math>A</math> --- амплитуда, <math>L</math> --- длина кристалла, <math>x\in [0, L]</math>.

[[Файл: Tcvet_masters_Theme3_IC.png|600px|Начальная скорость]]

Начальные перемещения равны нулю по всей длине кристалла.

==== Метод определения формы кристалла в различные моменты времени ====
При малых <math>t</math> в распределении скоростей по длине цепочки явно прослеживается форма синуса, однако, при возрастании <math>t</math> цепочка теряет свою форму и движение частиц неотличимо от теплового. Это случается из за нелинейного взаимодействия между частицами. Однако, с течением времени цепочка вновь принимает форму синуса, искривленного тепловым движением.

<gallery caption="Трансформации синуса" perrow=3 widths="250px" heights="250px">
Файл:Tcvet_masters_Theme3_sin_trans_T0.png|а) <math>T = 0</math>
Файл:Tcvet_masters_Theme3_sin_trans_T1100.png|б) <math>T = 15</math>
Файл:Tcvet_masters_Theme3_sin_trans_T4500.png|в) <math>T = 45</math>
Файл:Tcvet_masters_Theme3_sin_trans_T9000.png|г) <math>T = 90</math>
Файл:Tcvet_masters_Theme3_sin_trans_T18000.png|д) <math>T = 180</math>
Файл:Tcvet_masters_Theme3_sin_trans_T36035.png|е) <math>T = 360</math>
</gallery>



Чтобы более точно определить, насколько форма цепочки близка к форме синуса, посчитаем параметр <math>E^*</math> --- отношение полной механической энергии <math>E</math> к начальной полной механической энергии <math>E_0</math> по следующей формуле:

:<math>
  E^* = \left(\frac{\int\limits_{L} v(x)\sin\(\frac{2\pi x}{L}\)\,dx}
  {A\int\limits_{L} \sin^2\(\frac{2\pi x}{L}\)\,dx}\right)^2
  + \left(\frac{\int\limits_{L} D(x)\cos\(\frac{2\pi x}{L}\)\,dx}
  {A\int\limits_{L} \sin^2\(\frac{2\pi x}{L}\)\,dx}\right)^2.
</math>

<math>E^* \in [0, 1]</math>. Чем больше параметр <math>E^*</math>, тем больше форма кристалла похожа на синус.



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

<gallery caption="T --- период флуктуации синуса" perrow=3 widths="250px" heights="250px">
Файл:Tcvet_masters_Theme3_sin_T0.png|а) <math>T = 0</math>
Файл:Tcvet_masters_Theme3_sin_T75.png|б) <math>T = 0.125</math>
Файл:Tcvet_masters_Theme3_sin_T150.png|в) <math>T = 0.25</math>
Файл:Tcvet_masters_Theme3_sin_T225.png|г) <math>T = 0.375</math>
Файл:Tcvet_masters_Theme3_sin_T300.png|д) <math>T = 0.5</math>
Файл:Tcvet_masters_Theme3_sin_T600.png|е) <math>T = 1</math>
</gallery>



Использованы следующие параметры системы: амплитуда <math>A = 0.01</math>, количество частиц <math>N = 100</math>;

По графику ниже видно, что через некоторое время система возвращается к состоянию, близкому к начальному. Это означает, что распределение скоростей вновь приобретает форму синуса. Также видно, что система возвращается к начальному состоянию не на <math>100\%</math>, а с некоторым отклонением.

[[Файл: Tcvet_masters_Theme3_Temperature_short.png|600px|Результаты]]


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

[[Файл: Tcvet_masters_Theme3_Temperature_long.png|600px|Результаты, большое <math>T</math>]]




==== Зависимость времени возвращения кристалла к начальному состоянию от количества частиц и амплитуды синуса ====

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

Амплитуда <math>A = 0.01</math>.

Для того, чтобы глубже понять характер зависимости между временем, требуемым для возвращения кристалла к начальному состоянию, количеством частиц и амплитудой синуса, рассчитаны значения <math>T(N)</math> для разных амплитуд. Графики зависимостей <math>T(N)</math> представлены на рисунке ниже:

[[Файл: Tcvet_masters_Theme3_time_reverse(N,A)_6_lines.png|600px|<math>T(N)</math> для различных значений амплитуды]]

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



==== Высота второго пика ====
На рисунке ниже показана зависимость высоты второго пика функции <math>E^*(T)</math> от количества частиц в цепочке. Амплитуда <math>A = 0.01</math>.
Под вторым пиком подразумевается второй по времени локальный максимум функции <math>E^*(T)<math>, который означает, что после преобразования механической энергии кристалла в тепловую, произошло обратное преобразование, и тепловая энергия кристалла в максимальной степени преобразовалась в механическую. В этот момент распределение скоростей в кристалле вновь приобретает форму синуса.

По графику можно предположить, что со временем данное значение асимптотически стремится к нулю, но рассчитать его для большего <math>N</math> затруднительно, т.к. время, требуемое на расчет такой системы, растет нелинейно с увеличением количества частиц в цепочке.

[[Файл: Tcvet_masters_Theme3_second_peak_height.png|600px|Высота второго пика от <math>N</math>]]

=== Стохастическая задача ===
==== Начальные условия ====
==== Время перехода механической энергии в тепловую ====
== Заключение ==