Преобразование механической энергии в тепловую в одномерном кристалле — различия между версиями

Версия 12:05, 1 августа 2017

Выпускная квалификационная работа

Выполнил: студент группы 43604/1 Е. Б. Старобинский

Руководитель: доктор физ.-мат. наук, член-корр. РАН А. М. Кривцов

Введение

Cовременные нанотехнологии позволяют получать практически идеальные (бездефектные) материалы^[1][2], отличающиеся по своим свойствам от материалов с дефектами^[3]. В частности, тепловые процессы в таких наноструктурах протекают по иным, более сложным законам, чем для тел макроуровня. Знание этих законов и особенностей поведения наносистем имеет большое практическое значение при разработке новых устройств^[4][5] и расширении области применения наноматериалов, в том числе на промышленном уровне. Однако описание наносистем с помощью существующих моделей сопряжено с существенными противоречиями экспериментальным данным, что связано как с отсутствием цельной теории, позволяющей описывать процессы наноуровня, так и с трудоёмкостью проведения натурных экспериментов. Компромиссом является численное моделирование наноструктур, при котором поведение материала является прямым следствием уравнений динамики. Минус такого подхода – необходимость рассчитывать крупные фрагменты материала, чтобы за исследумое время возмущение не успевало дойти до границ. Следовательно, требуются большие вычислительные мощности.

В 1953 году на суперкомпьютере Maniac I группа учёных, состоящая из Энрико Ферми, Станислава Улама, Джона Паста и Мэри Цингоу, численно решила задачу Коши^[6] для системы дифференциальных уравнений, описывающую одномерный кристалл с нелинейным взаимодействием между частицами. С помощью синусоидальной волны задавалось начальное возмущение, ожидалось, что с течением времени энергия волны равномерно распределится по всем частицам. Оказалось, что сначала волна действительно теряет свою форму и энергию, но по прошествии достаточного времени почти точно возвращается к исходному состоянию. В дальнейшем это явление было названо «парадоксом Ферми-Паста-Улама-Цингоу», так как исследователям удалось доказать неспособность физики 20-го века описывать взаимодействия на наноуровне.

В данной работе рассматривается одна из нелинейных моделей кристалла, описанных Энрико Ферми. Ключевым отличием является введение тепловой энергии через задание дисперсии начальных скоростей.

Модель одномерного кристалла

Метод исследования

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

Для численного расчета формулируемых задач создан ряд специальных программ на языках С++, Bash и Mathematica.

Постановка задачи

Рассмотрим модель одномерного кристалла: цепочку одинаковых частиц массы $m$, соединённых одинаковыми нелинейными пружинами (квадратичная нелийность).

Каждой частице присвоим свой целочисленный индекс. Принимая $u_k$ за перемещение $k$-й частицы, а $C_1$ и $C_2$~--- за коэффициенты при линейном и квадратичном членах в разложении силы взаимодействия $k$-й частицы с двумя ближайшими соседями, получим следующее дифференциальное уравнение:

[math] \begin{equation} m\ddot u_k = C_1\left(\left(u_{k + 1} - u_{k}\right) - \left(u_{k} - u_{k - 1}\right)\right) + C_2\left(\left(u_{k + 1} - u_{k}\right)^2 - \left(u_{k} - u_{k - 1}\right)^2\right), \end{equation} [/math] где [math]\dot f[/math]~--- производная $f$ по времени

Перепишем это выражение, раскрыв скобки и введя новые обозначения для коэффициентов при линейной и нелинейной части:

\begin{equation} \ddot u_k = \omega_0^2\left(u_{k + 1} - 2u_{k} + u_{k - 1}\right) + \alpha_0^2\left(u_{k + 1}^2 - 2u_{k + 1}u_{k} + 2u_{k}u_{k - 1} - u_{k - 1}^2\right), \end{equation}

\begin{equation} \omega_0 = \sqrt{\frac{C_1}{m}}, \quad \alpha_0 = \sqrt{\frac{C_2}{m}}. \end{equation}

Здесь $\omega_0$ соответствует циклической частоте колебаний.

Заметим, что

\begin{equation} \begin{split} & u_{k + 1}^2 - 2u_{k + 1}u_{k} + 2u_{k}u_{k - 1} - u_{k - 1}^2 = \\ & = u_{k + 1}^2 - 2u_{k + 1}u_{k} + 2u_{k}u_{k - 1} - u_{k - 1}^2 \pm u_{k + 1}u_{k-1}= \\ & = \left(u_{k + 1}^2 - 2u_{k + 1}u_{k} + u_{k - 1}u_{k + 1}\right) - \left(u_{k + 1}u_{k - 1} - 2u_{k}u_{k - 1} + u_{k - 1}^2\right) = \\ & = u_{k + 1}\left(u_{k + 1} - 2u_{k} + u_{k - 1}\right) - u_{k - 1}\left(u_{k + 1} - 2u_{k} + u_{k - 1}\right) = \\ & = \left(u_{k + 1} - 2u_{k} + u_{k - 1}\right)\left( u_{k + 1} - u_{k - 1} \right). \end{split} \label{replace} \end{equation}

Воспользуемся заменой~(\ref{replace}), и тогда уравнения динамики цепочки примут следующий вид: \begin{equation} \ddot u_k = \left(u_{k + 1} - 2u_{k} + u_{k - 1} \right) \left( \omega_0^2 + \alpha_0^2 \left( u_{k + 1} - u_{k - 1} \right) \right). \label{chainEquations} \end{equation}

Аналогичное уравнение использовалось в оригинальной работе Ферми, Паста и Улама\cite{fermi1965crystal}.

Механическую энергию системы в начальный момент времени (то есть, при $t = 0$) будем задавать с помощью синусоидальной волны. Из-за теплового движения частиц волна будет терять свою форму, а её энергия перейдёт в тепловую энергию системы. Для описания этого процесса численно решим уравнения динамики кристалла.

Чтобы замкнуть систему из $k$ уравнений~(\ref{chainEquations}), определим граничные условия, а также по $k$ начальных условий на перемещения и скорости. Для задания начальных скоростей воспользуемся генератором случайных чисел с равномерным распределением. При этом дисперсию скоростей будем задавать в соответствии с требуемым температурным профилем.

Начальные и граничные условия

Чтобы соответствовать идеальной модели, исследуемая цепочка должна иметь бесконечную длину. Для проведения численного моделирования будем рассматривать кристалл конечной длины с заданным условием периодичности (граничные условия Борна~--- Кармана):

\begin{equation} \label{boundaryConditions} u\!\left(x + L\right) = u\!\left(x\right), \quad L \gg a, \end{equation} где $L$~--- длина кристалла, $a$~--- равновесное расстояние между частицами.

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

В рамках поставленной задачи будем рассматривать два вида движения в кристалле: движение механической волны и тепловое движение (тепловой шум). Соответственно, начальные условия для скоростей частиц также будут складываться из двух компонент. Меняя начальные условия, рассмотрим решение задачи для следующих постановок: \begin{itemize} \item на кристалле задана бегущая волна; \item на кристалле задана стоячая волна. \end{itemize}

Начальные условия для бегущей волны

Рассмотрим уравнение бегущей волны (пусть волна распространяется в одном направлении): \begin{equation} \label{runningWaveEquation1} \dot u\!\left(x\right) = A \sin\!\left(k_0 x - \omega_0 t + \phi \right), \end{equation} где $A$~--- амплитуда колебаний, $k_0$~--- волновое число, $\phi$~--- фаза колебаний.

Величину $u\!\left(x\right)$ на момент времени $t$ получим интегрированием уравнения~(\ref{runningWaveEquation1}):

\begin{equation} \label{runningWaveEquation2} u\!\left(x\right) = \int\limits_0^t \! \dot u\!\left(x\right) \,\mathrm{d}t = \frac{1}{\omega_0} A \cos\!\left(k_0 x - \omega_0 t + \phi \right). \end{equation}

Пусть длина волны совпадает с длиной кристалла, тогда \begin{equation} k_0~\mathrel{\stackrel{\rm def}=} \frac{2 \pi}{\lambda} = \frac{2 \pi}{L}, \end{equation} где $\lambda$~--- длина волны.

Фазу колебаний $\phi$ положим равной нулю и запишем уравнения~(\ref{runningWaveEquation1}) и~(\ref{runningWaveEquation2}) для начального момента времени ($t=0$): \begin{equation} \label{initialMechanicalVelocity1} \dot u\!\left(x\right)|_{t=0} = A \sin\!\left(\frac{2 \pi x}{L} \right). \end{equation} \begin{equation} \label{initialMechanicalDisplacement1} u\!\left(x\right)|_{t=0} = \frac{1}{\omega_0} A \cos\!\left(\frac{2 \pi x}{L}\right). \end{equation}

На рисунке~\ref{runningWaveInitialConditions} приведён эскиз начальных перемещений и скоростей частиц кристалла при заданной бегущей волне.

\begin{figure}[ht!] \center{\input{Data/InitialConditions/rw.tex}} \caption{Начальные условия для бегущей волны} \label{runningWaveInitialConditions} \end{figure}

Начальные условия для стоячей волны

Начальные условия~(\ref{initialMechanicalVelocity1}) и~(\ref{initialMechanicalDisplacement1}) позволяют задать бегущую волну. Получим подобные условия также для стоячей волны. Для этого рассмотрим две волны амплитуды $\frac{A}{2}$, бегущие в противоположных направлениях с разностью фаз $\pi$: \begin{equation} v_1\!\left(x\right) = \frac{A}{2} \sin\!\left(k_0 x - \omega_0 t \right), \quad v_2\!\left(x\right) = \frac{A}{2} \sin\!\left(- k_0 x - \omega_0 t + \pi \right), \end{equation}

Результатом наложения этих двух волн будет стоячая волна, описываемая следующим уравнением: \begin{equation} \label{standingWaveEquation1} \dot u\!\left(x\right) = v_1\!\left(x\right) + v_1\!\left(x\right) = A \cos\!\left(\omega_0 t \right)\sin\left(k_0 x \right). \end{equation}

Аналогично, проинтегрируем уравнение~(\ref{standingWaveEquation1}) по времени $t$: \begin{equation} u\!\left(x\right) = \int\limits_0^t \! \dot u\!\left(x\right) \,\mathrm{d}t = \frac{A}{\omega_0} \sin\!\left(\omega_0 t \right)\sin\left(k_0 x \right). \label{standingWaveEquation2} \end{equation}

Подставляя в~(\ref{standingWaveEquation1}) и~(\ref{standingWaveEquation2}) $k_0 = \frac{2 \pi}{L}$ и $t = 0$, получим начальные условия для стоячей волны: \begin{equation} \label{initialMechanicalVelocity2} \dot u\!\left(x\right)|_{t=0} = A \sin\!\left(\frac{2 \pi x}{L} \right). \end{equation} \begin{equation} \label{initialMechanicalDisplacement2} u\!\left(x\right)|_{t=0} = 0. \end{equation}

Соответствующий график начальных перемещений и скоростей частиц цепочки приведён на рисунке~\ref{standingWaveInitialConditions}.

\begin{figure}[ht!] \center{\input{Data/InitialConditions/sw.tex}} \caption{Начальные условия для стоячей волны} \label{standingWaveInitialConditions} \end{figure}

Функции скоростей в начальный момент времени для бегущей~(\ref{initialMechanicalVelocity1}) и стоячей~(\ref{initialMechanicalVelocity2}) волн совпадают, различаются только начальные перемещения.

Начальные условия теплового движения

Введём оператор усреднения функции по её аргументу: \begin{equation} \left<f\right> = \frac{1}{b - a}\int\limits_{a}^{b} \! f \! \left( x \right) \, \mathrm{d}x, \end{equation} где $\left[a; b\right]$~--- отрезок, на котором проводится усреднение.

Оператор нахождения среднего по координате для цепочки длины $L$ в таком случае записывается следующим образом: \begin{equation} \left<f\right> = \frac{1}{L}\int\limits_{0}^{L} \! f \! \left( x \right) \, \mathrm{d}x, \end{equation}

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

\begin{equation} \frac{1}{2}kT\!\left(x\right) = \frac{m \left< v^{2} \right>}{2}, \end{equation} где $k$~--- постоянная Больцмана, $T(x)$~--- профиль температуры, $v$~--- скорость.

Из этого равенства выразим $\left< v^{2} \right>$: \begin{equation} \left< v^{2} \right> = \frac{kT\!\left(x\right)}{m}. \label{vTemp2Average} \end{equation}

Пусть в начальный момент времени температура в кристалле распределена равномерно. Соответствующее распределение для скоростей будем задавать с помощью случайного числа $\rho$ в диапазоне $\left[-1;1\right]$. Определим среднее значение $\rho^{2}$: \begin{equation} \left< \rho^{2} \right> = \frac{1}{2}\int\limits_{-1}^{1} \! \rho^{2} \, \mathrm{d}\rho = \frac{1}{2} \frac{\rho^{3}}{3} \bigg|_{-1}^{1} = \frac{1}{3}. \end{equation}

Таким образом, $3 \left< \rho^{2} \right> = 1$. Воспользуемся этим соотношением для построения требуемого профиля $T \! \left( x \right)$: \begin{equation} T \! \left( x \right) = 3 \left< \rho^{2} \right> T_0, \end{equation} где $T_0$~--- значение температуры при $t = 0$.

Перепишем выражение~(\ref{vTemp2Average}): \begin{equation} \left< v^{2} \right> = \frac{ 3 k T_0}{m} \left< \rho^{2} \right>. \end{equation}

Выразим $v^{2}\!\left(x\right)$ и $v\!\left(x\right)$ \begin{equation} v^{2}\!\left(x\right) = \frac{ 3 k T_0}{m} \rho^{2}\!\left(x\right). \end{equation}

\begin{equation} v\!\left(x\right) = \sqrt{\frac{3 k T_0}{m}} \rho\!\left(x\right). \label{vTemp} \end{equation}

Введём величину дисперсии скоростей $\sigma_{v}^{2}$ так, чтобы при $\sigma_{v}^{2} = 1$ в начальный момент времени механическая энергия волны равнялась тепловой энергии цепочки: \begin{equation} T_0 = \frac{m}{k} A^{2} \sigma_{v}^{2} \label{initialTemp} \end{equation}

Тогда значение $\sigma_{v}^{2} = 10$ можно трактовать следующим образом: при $t~=~0$ энергия теплового шума в кристалле в $10$ раз превосходит механическую энергию волны.

Перепишем~(\ref{vTemp}) с учётом соотношения~(\ref{initialTemp}):

\begin{equation} v\!\left(x\right)|_{t = 0} = A \sqrt{3 \sigma_{v}^{2}} \rho \! \left( x \right). \label{initialThermalVelocity} \end{equation}

Таким образом, были получены выражения для обеих компонент скорости: скорость механической волны задаётся функцией~(\ref{initialMechanicalVelocity1}), тепловой шум определяется функцией (\ref{initialThermalVelocity}). Итоговое условие на начальные скорости в присталле принимает следующий вид: \begin{equation} \left. \dot u(x) \right |_{t = 0} = A \left( \sin\!\left( \frac{2 \pi x}{L} \right) + \sqrt{3 \sigma_{v}^{2}} \rho \! \left( x \right) \right). \end{equation}

Начальные перемещения при рассмотрении бегущей волны задаются формулой~(\ref{initialMechanicalDisplacement1}), при рассмотрении стоячей~--- формулой~(\ref{initialMechanicalDisplacement2})

Вычисление механической энергии

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

Механическая энергия системы $E\!\left(x\right)$ складывается из кинетической энергии $K\!\left(x\right)$ и потенциальной энергии $\Pi\!\left(x\right)$ следующий образом: \begin{equation} K\!\left( x \right) = \frac{m}{2}\dot u^2\!\left(x\right), \quad \Pi \!\left( x \right) = \frac{C_1}{2}\varepsilon^2\!\left(x\right), \quad E\!\left(x\right) = K\!\left(x\right) + \Pi\!\left(x\right), \label{energyEquations} \end{equation} где $\varepsilon\!\left(x\right)$~--- возникающие деформации.

Так как нас интересует энергия не всей системы (включающей тепловой шум), а только заданной механической волны, будем искать энергию первой синусоидальной моды. Перепишем выражения для энергий: \begin{equation} K\!\left( x \right) = \frac{m}{2}\left(\int\limits_L \! \dot u\!\left(x\right) \sin\!\left( \frac{2 \pi x}{L} \right) \,\mathrm{d}x\right)^2 \end{equation} \begin{equation} \Pi \!\left( x \right) = \frac{C_1}{2}\left(\int\limits_L \! \varepsilon\!\left(x\right) \cos\!\left( \frac{2 \pi x}{L} \right) \,\mathrm{d}x\right)^2 \end{equation} \begin{equation} E\!\left(x\right) = \frac{m}{2}\left(\int\limits_L \! \dot u\!\left(x\right) \sin\!\left( \frac{2 \pi x}{L} \right) \,\mathrm{d}x\right)^2 + \frac{C_1}{2}\left(\int\limits_L \! \varepsilon\!\left(x\right) \cos\!\left( \frac{2 \pi x}{L} \right) \,\mathrm{d}x\right)^2 \end{equation}

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

Определим значения $K\!\left(x\right)$ и $\Pi\!\left(x\right)$ при $t = 0$:

\begin{equation} K\!\left( x \right)|_{t=0} = \frac{m}{2}\left(A \int\limits_L \! \sin^2\!\left( \frac{2 \pi x}{L} \right) \,\mathrm{d}x\right)^2, \end{equation}

\begin{equation} \label{initialPotentialEnergy} \Pi\!\left( x \right)|_{t=0} = \frac{C_1}{2}\left( \int\limits_L \! \varepsilon\!\left(x\right)|_{t = 0} \cos\!\left( \frac{2 \pi x}{L} \right) \,\mathrm{d}x\right)^2 \end{equation}

Интеграл в формуле (\ref{initialPotentialEnergy}) равен нулю в случае стоячей волны (в начальный момент времени нет деформаций). Покажем, что этот интеграл также будет равен нулю в случае бегущей волны:

\begin{equation} \begin{split} & \int\limits_L \! \varepsilon\!\left(x\right)|_{t = 0} \cos\!\left( \frac{2 \pi x}{L} \right) \,\mathrm{d}x = A \int\limits_L \! \sin\!\left(\frac{2 \pi x}{L}\right) \cos\!\left( \frac{2 \pi x}{L} \right) \,\mathrm{d}x = \\ & = \frac{A}{2} \int\limits_L \! \sin\!\left(\frac{4 \pi x}{L}\right) \,\mathrm{d}x = 0. \end{split} \end{equation}

То есть, в начальный момент времени потенциальная энергия равна нулю, а механическая энергия задаваемой волны равна \begin{equation} \label{initialMechanicalEnergy} E\!\left( x \right)|_{t=0} = \frac{m}{2}\left(A \int\limits_L \! \sin^2\!\left( \frac{2 \pi x}{L} \right) \,\mathrm{d}x\right)^2 = \frac{m A^{2} L^{2}}{8} \end{equation}

Значение (\ref{initialMechanicalEnergy}) выберем в качестве масштаба энергии.

Чтобы определить закон, по которому происходит переход механической энергии в тепловую, введём безразмерный параметр $E^*\!\left( x \right) = \frac{E\!\left( x \right)}{E\!\left( x \right)|_{t=0}}$~--- нормированное значение механической энергии цепочки: \begin{equation} E^* = \left( \frac{\int\limits_L \! \dot u\!\left(x\right) \sin\!\left( \frac{2 \pi x}{L} \right) \,\mathrm{d}x}{A \int\limits_L \! \sin^2 \left( \frac{2 \pi x}{L} \right) \,\mathrm{d}x} \right)^2 + \omega_0^2 \left( \frac{\int\limits_L \! \varepsilon\!\left(x\right) \cos\!\left( \frac{2 \pi x}{L} \right) \,\mathrm{d}x}{A \int\limits_L \! \sin^2 \left( \frac{2 \pi x}{L} \right) \,\mathrm{d}x} \right)^2. \end{equation}

Решение задачи свелось к нахождению $E^*\!\left( x \right)$.

Численное решение дифференциального уравнения

Для решения системы уравнений движения~(\ref{chainEquations}) воспользуемся методом центральных разностей. Перепишем граничные условия~(\ref{boundaryConditions}): \begin{equation} u_{k + N} = u_{k}, \quad \dot u_{k + N} = \dot u_{k}, \end{equation} где $N \gg 1$~--- число независимых частиц.

Вычисление правой части~(\ref{chainEquations}) аналогично решению с помощью разностной схемы. Рассчитанные таким образом ускорения частиц пересчитываем в скорости: \begin{equation} \dot u_k\!\left(t + \Delta t\right) = \dot u_k\!\left(t\right) + \ddot u_k\!\left(t\right)\Delta t. \end{equation}

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

\begin{equation} u_k\!\left(t + \Delta t\right) = u_k\!\left(t\right) + \dot u_k\!\left(t\right)\Delta t. \end{equation}

Значение $\varepsilon\!\left(x\right)$ вычисляем следующим образом: \begin{equation} \varepsilon = \sum_{i = 1}^{N - 1} \left(u_{i} - u_{i - 1}\right). \end{equation}

Параметры расчёта

Для достоверного численного решения дифференциального уравнения относительно $\left(u\!\left(x\right), \dot u\!\left(x\right)\right)$ длина цепочки должна стремиться к бесконечной, чтобы соответствовать модели идеального кристалла. При проведении численных экспериментов удовлетворительная точность расчётов достигалась на цепочках из 10 миллионов частиц и более.

Статистические характеристики системы также можно вычислить путём усреднения бесконечного множества реализаций (каждая реализация~--- задача Коши для цепочки). Благодаря этому подходу, предложенному в работах \cite{art1}, \cite{art2}, достигается эффективность проводимых численных экспериментов.

Результаты, приводимые далее, получены усреднением $5 \cdot 10^4$ реализаций цепочки c $N = 10^3$ (за исключением $\rho$ все параметры в рамках серии сохранялись). Такая конфигурация будет аналогична модели кристалла с 50 миллионами частиц, но допускает параллельный компьютерный расчёт реализаций. Проведение параллельного моделирования коротких цепочек на суперкомпьютерных системах привело к сушественному сокращению общего времени вычислений.

По мере повышения величины дисперсии $\sigma_v^2$ шаг по времени $\Delta t$ следует уменьшать. Иначе влияние нелинейной компоненты силы приводит к разрушению системы, эта неустойчивость связана с нарушением условия Куранта. Однако, уменьшение $\Delta t$ увеличивает время, требуемое для проведения численного моделирования. Так как численно $\sigma_v^2$ соответствует отношению тепловой энергии в кристалле к энергии механической волны, её значение в реальных системах на несколько порядков превышает исследуемое в рамках данной работы.

В качестве масштаба времени возьмём период колебаний материальной точки массы $m$ на пружине жёсткости $C_1$: \begin{equation} \quad T_0~\mathrel{\stackrel{\rm def}=} \frac{2 \pi}{\omega_0} = 2 \pi \sqrt{\frac{m}{C_1}}. \end{equation}

\begin{center} \captionof{table}{Параметры при моделировании стоячей волны} \begin{tabular}{l | l | l | l} $\sigma_v^2$ & $\Delta \sigma_v^2$ & $\Delta t$ & Постоянные параметры\\ \hline $[0; 60]$ & $1$ & $0.05~T_0$ & \multirow{4}{*}{$\omega_0 = 1, \quad \alpha_0 = 1, \quad A = 0.01$} \\ $[15; 60]$ & $5$ & $0.05~T_0$ & \\ $[65; 75]$ & $5$ & $0.005~T_0$ & \\ $[80; 100]$ & $5$ & $0.0005~T_0$ & \end{tabular} \label{standingWaveParameters} \end{center}

\begin{center} \captionof{table}{Параметры при моделировании бегущей волны} \begin{tabular}{l | l | l | l} $\sigma_v^2$ & $\Delta \sigma_v^2$ & $\Delta t$ & Постоянные параметры\\ \hline $[0; 10]$ & $1$ & $0.05~T_0$ & \multirow{3}{*}{$\omega_0 = 1, \quad \alpha_0 = 1, \quad A = 0.01$} \\ $[15; 200]$ & $5$ & $0.05~T_0$ & \\ $[200; 300]$ & $5$ & $0.005~T_0$ & \end{tabular} \label{runningWaveParameters} \end{center}

Значения параметров проводимых экспериментов приведены в таблицах~\ref{standingWaveParameters} и~\ref{runningWaveParameters}. Моделируемое время ограничено: $0 \leq t \leq 7 \cdot 10^4~T_0$.

Заключение

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

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

В качестве меры механической энергии системы [math]E^*[/math] взято отношение энергии волны к её начальному значению. Для проведения исследования введён безразмерный параметр [math]\sigma_v^2[/math], характеризующий величину дисперсии флуктуаций скоростей в цепочке и равный отношению тепловой энергии к механической в начале расчёта.

Заданная механическая энергия с течением времени переходит в тепловую. Продемонстрирована необратимость этого процесса, а также зависимость скорости перехода от [math]\sigma_v^2[/math]: при увеличении дисперсии механическая энергия убывает быстрее. Функция убывания описана с помощью модифицированного экспоненциального закона: [math]E^* = e^{-\alpha \tau ^{\beta}}[/math], где [math]\tau[/math] – безразмерное время процесса, [math]\beta[/math] и [math]\alpha[/math]– функции от [math]\sigma_v^2[/math].

Зависимость [math]\alpha \! \left( \sigma_v^2 \right)[/math] и [math]\beta \! \left( \sigma_v^2 \right)[/math] получены экспериментально для серии [math]\sigma_v^2[/math]. Предложены аппроксимации этих функций для стойчей и бегущей волн.

[math] \beta \! \left(\sigma_{v}^{2} \right) = 1 + \frac{1.08}{\sqrt{\sigma_{v}}}, \quad \alpha \! \left(\sigma_{v}^{2} \right) = 0.14 \exp\!\left(- \frac{50}{5 + \sigma_{v}}\right), [/math]

[math] \beta \! \left(\sigma_{v}^{2} \right) = 1 + \frac{1.55}{\sqrt{\sigma_{v}}}, \quad \alpha \! \left(\sigma_{v}^{2} \right) = 0.14 \exp\!\left(- \frac{97.5}{9.5 + \sigma_{v}}\right), [/math]

где [math]\sigma_{v}[/math] – среднеквадратическое отклонение скоростей частиц кристалла.

Показано, что энергия стоячей волны убывает в 4 раза быстрее энергии бегущей волны при одних параметрах эксперимента. Продемонстрированы предельные значения [math]\beta[/math]:

• [math]\beta = 2[/math] при [math]\sigma_v^2 = 0[/math];
• [math]\beta \rightarrow 1[/math] при [math]\sigma_v^2 \rightarrow \infty[/math].

Показано влияние скорости убывания [math]E^*[/math] на развал бегущей волны и, как следствие, образование дополнительного колебательного процесса. При больших значениях [math]\sigma_v^2[/math] волна сохраняет свою первоначальную форму, развал не происходит из-за слишком быстрого перехода механической энергии волны в тепловую.

Результаты работы были получены с использованием вычислительных ресурсов суперкомпьютерного центра Санкт-Петербургского политехнического университета Петра Великого. Каждый расчёт производился на отдельном ядре CPU с использованием технологии MPI, одновременно задействовалось [math]1000[/math] ядер. Таким образом, для получения [math]5 \cdot 10^4[/math] реализаций вычисления повторялись [math]50[/math] раз, в среднем общее время моделирования для одного значения [math]\sigma^{2}_{v}[/math] составило [math]20[/math] минут при шаге [math]\Delta t = 0.05 T_0[/math], [math]6[/math] часов при [math]\Delta t = 0.005 T_0[/math] и [math]12[/math] часов при [math]\Delta t = 0.0005 T_0[/math] ([math]T_0[/math] – заданный масштаб времени).

Автор благодарен Д. В. Цветкову за полезные обсуждения.

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

1. ↑ Yenny et al. Hernandez. High-yield production of graphene by liquid-phase exfoliation of graphite. Nature Nanotechnology, 3:563–568, 2008.
2. ↑ Xiaolin et al. Li. Highly conducting graphene sheets and langmuir-blodgett films. Nature Nanotechnology, 3:538–542, 2008.
3. ↑ L. Shi et al. Evaluating broader impacts of nanoscale thermal transport research. Nanoscale and Microscale Thermophysical Engineering, 19:127–165, 2015.
4. ↑ Kostya S Novoselov, Andre K Geim, Sergei V Morozov, D Jiang, Y Zhang, Sergey V Dubonos, Irina V Grigorieva, and Alexandr A Firsov. Electric field effect in atomically thin carbon films. Science, 306(5696):666–669, 2004.
5. ↑ Z. Xu, G. Tai, Y. Zhou, F. Gao, and K. H. Wong. Self-charged graphene battery harvests electricity from thermal energy of the environment. ArXiv e-prints, 2012.
6. ↑ Enrico Fermi, J. Pasta, and S. Ulam. Studies of nonlinear problems. Los Alamos Report LA-1940, 978, 1965.