Редактирование: Преобразование механической энергии в тепловую в одномерном кристалле
Внимание! Вы не авторизовались на сайте. Ваш IP-адрес будет публично видимым, если вы будете вносить любые правки. Если вы войдёте или создадите учётную запись, правки вместо этого будут связаны с вашим именем пользователя, а также у вас появятся другие преимущества.
Правка может быть отменена. Пожалуйста, просмотрите сравнение версий, чтобы убедиться, что это именно те изменения, которые вас интересуют, и нажмите «Записать страницу», чтобы изменения вступили в силу.
Текущая версия | Ваш текст | ||
Строка 1: | Строка 1: | ||
'''''Выпускная квалификационная работа''''' | '''''Выпускная квалификационная работа''''' | ||
− | |||
− | |||
'''Выполнил:''' студент группы 43604/1 [[Старобинский Егор|Е. Б. Старобинский]] | '''Выполнил:''' студент группы 43604/1 [[Старобинский Егор|Е. Б. Старобинский]] | ||
− | ''' | + | '''Руководитель:''' доктор физ.-мат. наук, член-корр. РАН [[Кривцов Антон | А. М. Кривцов]] |
− | |||
− | |||
− | |||
− | |||
− | |||
== Введение == | == Введение == | ||
Cовременные нанотехнологии позволяют получать практически идеальные (бездефектные) материалы<ref>Yenny et al. Hernandez. High-yield production of graphene by liquid-phase exfoliation of graphite. Nature Nanotechnology, 3:563–568, 2008.</ref><ref>Xiaolin et al. Li. Highly conducting graphene sheets and langmuir-blodgett films. Nature Nanotechnology, 3:538–542, 2008.</ref>, отличающиеся по своим свойствам от материалов с дефектами<ref>L. Shi et al. Evaluating broader impacts of nanoscale thermal transport research. Nanoscale and Microscale Thermophysical Engineering, 19:127–165, 2015.</ref>. В частности, тепловые процессы в таких наноструктурах протекают по иным, более сложным законам, чем для тел макроуровня. Знание этих законов и особенностей поведения наносистем имеет большое практическое значение при разработке новых устройств<ref>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.</ref><ref>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.</ref> и расширении области применения наноматериалов, в том числе на промышленном уровне. Однако описание наносистем с помощью существующих моделей сопряжено с существенными противоречиями экспериментальным данным, что связано как с отсутствием цельной теории, позволяющей описывать процессы наноуровня, так и с трудоёмкостью проведения натурных экспериментов. Компромиссом является численное моделирование наноструктур, при котором поведение материала является прямым следствием уравнений динамики. Минус такого подхода – необходимость рассчитывать крупные фрагменты материала, чтобы за исследумое время возмущение не успевало дойти до границ. Следовательно, требуются большие вычислительные мощности. | Cовременные нанотехнологии позволяют получать практически идеальные (бездефектные) материалы<ref>Yenny et al. Hernandez. High-yield production of graphene by liquid-phase exfoliation of graphite. Nature Nanotechnology, 3:563–568, 2008.</ref><ref>Xiaolin et al. Li. Highly conducting graphene sheets and langmuir-blodgett films. Nature Nanotechnology, 3:538–542, 2008.</ref>, отличающиеся по своим свойствам от материалов с дефектами<ref>L. Shi et al. Evaluating broader impacts of nanoscale thermal transport research. Nanoscale and Microscale Thermophysical Engineering, 19:127–165, 2015.</ref>. В частности, тепловые процессы в таких наноструктурах протекают по иным, более сложным законам, чем для тел макроуровня. Знание этих законов и особенностей поведения наносистем имеет большое практическое значение при разработке новых устройств<ref>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.</ref><ref>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.</ref> и расширении области применения наноматериалов, в том числе на промышленном уровне. Однако описание наносистем с помощью существующих моделей сопряжено с существенными противоречиями экспериментальным данным, что связано как с отсутствием цельной теории, позволяющей описывать процессы наноуровня, так и с трудоёмкостью проведения натурных экспериментов. Компромиссом является численное моделирование наноструктур, при котором поведение материала является прямым следствием уравнений динамики. Минус такого подхода – необходимость рассчитывать крупные фрагменты материала, чтобы за исследумое время возмущение не успевало дойти до границ. Следовательно, требуются большие вычислительные мощности. | ||
− | В 1953 году на суперкомпьютере Maniac I группа учёных, состоящая из Энрико Ферми, Станислава Улама, Джона Паста и Мэри Цингоу, численно решила задачу Коши<ref | + | В 1953 году на суперкомпьютере Maniac I группа учёных, состоящая из Энрико Ферми, Станислава Улама, Джона Паста и Мэри Цингоу, численно решила задачу Коши<ref>Enrico Fermi, J. Pasta, and S. Ulam. Studies of nonlinear problems. Los Alamos Report LA-1940, 978, 1965.</ref> для системы дифференциальных уравнений, описывающую одномерный кристалл с нелинейным взаимодействием между частицами. С помощью синусоидальной волны задавалось начальное возмущение, ожидалось, что с течением времени энергия волны равномерно распределится по всем частицам. Оказалось, что сначала волна действительно теряет свою форму и энергию, но по прошествии достаточного времени почти точно возвращается к исходному состоянию. В дальнейшем это явление было названо «парадоксом Ферми-Паста-Улама-Цингоу», так как исследователям удалось доказать неспособность физики 20-го века описывать взаимодействия на наноуровне. |
В данной работе рассматривается одна из нелинейных моделей кристалла, описанных Энрико Ферми. Ключевым отличием является введение тепловой энергии через задание дисперсии начальных скоростей. | В данной работе рассматривается одна из нелинейных моделей кристалла, описанных Энрико Ферми. Ключевым отличием является введение тепловой энергии через задание дисперсии начальных скоростей. | ||
Строка 29: | Строка 22: | ||
=== Постановка задачи === | === Постановка задачи === | ||
− | Рассмотрим модель одномерного кристалла: цепочку одинаковых частиц массы | + | Рассмотрим модель одномерного кристалла: цепочку одинаковых частиц массы $m$, соединённых одинаковыми нелинейными пружинами (квадратичная нелийность). |
− | Каждой частице присвоим свой целочисленный индекс. Принимая | + | Каждой частице присвоим свой целочисленный индекс. Принимая $u_k$ за перемещение $k$-й частицы, а $C_1$ и $C_2$~--- за коэффициенты при линейном и квадратичном членах в разложении силы взаимодействия $k$-й частицы с двумя ближайшими соседями, получим следующее дифференциальное уравнение: |
<math> | <math> | ||
Строка 39: | Строка 32: | ||
\end{equation} | \end{equation} | ||
</math> | </math> | ||
− | + | где <math>\dot f</math>~--- производная $f$ по времени | |
− | где <math>\dot f</math> | ||
Перепишем это выражение, раскрыв скобки и введя новые обозначения для коэффициентов при линейной и нелинейной части: | Перепишем это выражение, раскрыв скобки и введя новые обозначения для коэффициентов при линейной и нелинейной части: | ||
− | |||
\begin{equation} | \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), | \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} | \end{equation} | ||
− | |||
− | |||
\begin{equation} | \begin{equation} | ||
\omega_0 = \sqrt{\frac{C_1}{m}}, \quad \alpha_0 = \sqrt{\frac{C_2}{m}}. | \omega_0 = \sqrt{\frac{C_1}{m}}, \quad \alpha_0 = \sqrt{\frac{C_2}{m}}. | ||
\end{equation} | \end{equation} | ||
− | |||
− | Здесь | + | Здесь $\omega_0$ соответствует циклической частоте колебаний. |
Заметим, что | Заметим, что | ||
− | |||
\begin{equation} | \begin{equation} | ||
\begin{split} | \begin{split} | ||
Строка 71: | Строка 58: | ||
\label{replace} | \label{replace} | ||
\end{equation} | \end{equation} | ||
− | |||
− | |||
− | |||
− | + | Воспользуемся заменой~(\ref{replace}), и тогда уравнения динамики цепочки примут следующий вид: | |
\begin{equation} | \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). | \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} | \label{chainEquations} | ||
\end{equation} | \end{equation} | ||
− | |||
− | Аналогичное уравнение использовалось в оригинальной работе Ферми, Паста и Улама | + | Аналогичное уравнение использовалось в оригинальной работе Ферми, Паста и Улама\cite{fermi1965crystal}. |
− | Механическую энергию системы в начальный момент времени (то есть, при | + | Механическую энергию системы в начальный момент времени (то есть, при $t = 0$) будем задавать с помощью синусоидальной волны. Из-за теплового движения частиц волна будет терять свою форму, а её энергия перейдёт в тепловую энергию системы. Для описания этого процесса численно решим уравнения динамики кристалла. |
− | Чтобы замкнуть систему из | + | Чтобы замкнуть систему из $k$ уравнений~(\ref{chainEquations}), определим граничные условия, а также по $k$ начальных условий на перемещения и скорости. Для задания начальных скоростей воспользуемся генератором случайных чисел с равномерным распределением. При этом дисперсию скоростей будем задавать в соответствии с требуемым температурным профилем. |
=== Начальные и граничные условия === | === Начальные и граничные условия === | ||
− | Чтобы соответствовать идеальной модели, исследуемая цепочка должна иметь бесконечную длину. Для проведения численного моделирования будем рассматривать кристалл конечной длины с заданным условием периодичности (граничные условия Борна-Кармана): | + | Чтобы соответствовать идеальной модели, исследуемая цепочка должна иметь бесконечную длину. Для проведения численного моделирования будем рассматривать кристалл конечной длины с заданным условием периодичности (граничные условия Борна~--- Кармана): |
− | |||
\begin{equation} | \begin{equation} | ||
\label{boundaryConditions} | \label{boundaryConditions} | ||
u\!\left(x + L\right) = u\!\left(x\right), \quad L \gg a, | u\!\left(x + L\right) = u\!\left(x\right), \quad L \gg a, | ||
\end{equation} | \end{equation} | ||
− | + | где $L$~--- длина кристалла, $a$~--- равновесное расстояние между частицами. | |
− | |||
− | где | ||
Таким образом, условие теплоизоляции выполняется. | Таким образом, условие теплоизоляции выполняется. | ||
В рамках поставленной задачи будем рассматривать два вида движения в кристалле: движение механической волны и тепловое движение (тепловой шум). Соответственно, начальные условия для скоростей частиц также будут складываться из двух компонент. Меняя начальные условия, рассмотрим решение задачи для следующих постановок: | В рамках поставленной задачи будем рассматривать два вида движения в кристалле: движение механической волны и тепловое движение (тепловой шум). Соответственно, начальные условия для скоростей частиц также будут складываться из двух компонент. Меняя начальные условия, рассмотрим решение задачи для следующих постановок: | ||
− | + | \begin{itemize} | |
− | + | \item на кристалле задана бегущая волна; | |
+ | \item на кристалле задана стоячая волна. | ||
+ | \end{itemize} | ||
==== Начальные условия для бегущей волны ==== | ==== Начальные условия для бегущей волны ==== | ||
− | |||
Рассмотрим уравнение бегущей волны (пусть волна распространяется в одном направлении): | Рассмотрим уравнение бегущей волны (пусть волна распространяется в одном направлении): | ||
− | |||
− | |||
\begin{equation} | \begin{equation} | ||
\label{runningWaveEquation1} | \label{runningWaveEquation1} | ||
\dot u\!\left(x\right) = A \sin\!\left(k_0 x - \omega_0 t + \phi \right), | \dot u\!\left(x\right) = A \sin\!\left(k_0 x - \omega_0 t + \phi \right), | ||
\end{equation} | \end{equation} | ||
− | + | где $A$~--- амплитуда колебаний, $k_0$~--- волновое число, $\phi$~--- фаза колебаний. | |
− | |||
− | где | ||
− | Величину | + | Величину $u\!\left(x\right)$ на момент времени $t$ получим интегрированием уравнения~(\ref{runningWaveEquation1}): |
− | |||
\begin{equation} | \begin{equation} | ||
\label{runningWaveEquation2} | \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). | 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} | \end{equation} | ||
− | |||
Пусть длина волны совпадает с длиной кристалла, тогда | Пусть длина волны совпадает с длиной кристалла, тогда | ||
− | |||
− | |||
\begin{equation} | \begin{equation} | ||
k_0~\mathrel{\stackrel{\rm def}=} \frac{2 \pi}{\lambda} = \frac{2 \pi}{L}, | k_0~\mathrel{\stackrel{\rm def}=} \frac{2 \pi}{\lambda} = \frac{2 \pi}{L}, | ||
\end{equation} | \end{equation} | ||
− | + | где $\lambda$~--- длина волны. | |
− | |||
− | где | ||
− | Фазу колебаний | + | Фазу колебаний $\phi$ положим равной нулю и запишем уравнения~(\ref{runningWaveEquation1}) и~(\ref{runningWaveEquation2}) для начального момента времени ($t=0$): |
− | |||
− | |||
\begin{equation} | \begin{equation} | ||
\label{initialMechanicalVelocity1} | \label{initialMechanicalVelocity1} | ||
\dot u\!\left(x\right)|_{t=0} = A \sin\!\left(\frac{2 \pi x}{L} \right). | \dot u\!\left(x\right)|_{t=0} = A \sin\!\left(\frac{2 \pi x}{L} \right). | ||
\end{equation} | \end{equation} | ||
− | |||
− | |||
− | |||
\begin{equation} | \begin{equation} | ||
\label{initialMechanicalDisplacement1} | \label{initialMechanicalDisplacement1} | ||
u\!\left(x\right)|_{t=0} = \frac{1}{\omega_0} A \cos\!\left(\frac{2 \pi x}{L}\right). | u\!\left(x\right)|_{t=0} = \frac{1}{\omega_0} A \cos\!\left(\frac{2 \pi x}{L}\right). | ||
\end{equation} | \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} | \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), | 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} | \end{equation} | ||
− | |||
Результатом наложения этих двух волн будет стоячая волна, описываемая следующим уравнением: | Результатом наложения этих двух волн будет стоячая волна, описываемая следующим уравнением: | ||
− | |||
− | |||
\begin{equation} | \begin{equation} | ||
\label{standingWaveEquation1} | \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). | \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} | \end{equation} | ||
− | |||
− | |||
− | |||
− | + | Аналогично, проинтегрируем уравнение~(\ref{standingWaveEquation1}) по времени $t$: | |
\begin{equation} | \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). | 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} | \label{standingWaveEquation2} | ||
\end{equation} | \end{equation} | ||
− | |||
− | Подставляя в | + | Подставляя в~(\ref{standingWaveEquation1}) и~(\ref{standingWaveEquation2}) $k_0 = \frac{2 \pi}{L}$ и $t = 0$, получим начальные условия для стоячей волны: |
− | |||
− | |||
\begin{equation} | \begin{equation} | ||
\label{initialMechanicalVelocity2} | \label{initialMechanicalVelocity2} | ||
\dot u\!\left(x\right)|_{t=0} = A \sin\!\left(\frac{2 \pi x}{L} \right). | \dot u\!\left(x\right)|_{t=0} = A \sin\!\left(\frac{2 \pi x}{L} \right). | ||
\end{equation} | \end{equation} | ||
− | |||
− | |||
− | |||
\begin{equation} | \begin{equation} | ||
\label{initialMechanicalDisplacement2} | \label{initialMechanicalDisplacement2} | ||
u\!\left(x\right)|_{t=0} = 0. | u\!\left(x\right)|_{t=0} = 0. | ||
\end{equation} | \end{equation} | ||
− | |||
− | Соответствующий график начальных перемещений и скоростей частиц цепочки приведён на рисунке | + | Соответствующий график начальных перемещений и скоростей частиц цепочки приведён на рисунке~\ref{standingWaveInitialConditions}. |
+ | |||
+ | \begin{figure}[ht!] | ||
+ | \center{\input{Data/InitialConditions/sw.tex}} | ||
+ | \caption{Начальные условия для стоячей волны} | ||
+ | \label{standingWaveInitialConditions} | ||
+ | \end{figure} | ||
− | + | Функции скоростей в начальный момент времени для бегущей~(\ref{initialMechanicalVelocity1}) и стоячей~(\ref{initialMechanicalVelocity2}) волн совпадают, различаются только начальные перемещения. | |
− | |||
==== Начальные условия теплового движения ==== | ==== Начальные условия теплового движения ==== | ||
Введём оператор усреднения функции по её аргументу: | Введём оператор усреднения функции по её аргументу: | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 218: | Строка 175: | ||
\end{equation} | \end{equation} | ||
</math> | </math> | ||
− | + | где <math>\left[a; b\right]</math>~--- отрезок, на котором проводится усреднение. | |
− | где <math>\left[a; b\right]</math> | ||
Оператор нахождения среднего по координате для цепочки длины <math>L</math> в таком случае записывается следующим образом: | Оператор нахождения среднего по координате для цепочки длины <math>L</math> в таком случае записывается следующим образом: | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 236: | Строка 191: | ||
\end{equation} | \end{equation} | ||
</math> | </math> | ||
− | + | где <math>k</math>~--- постоянная Больцмана, <math>T(x)</math>~--- профиль температуры, <math>v</math>~--- скорость. | |
− | где <math>k</math> | ||
Из этого равенства выразим <math>\left< v^{2} \right></math>: | Из этого равенства выразим <math>\left< v^{2} \right></math>: | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 257: | Строка 210: | ||
Таким образом, <math>3 \left< \rho^{2} \right> = 1</math>. Воспользуемся этим соотношением для построения требуемого профиля <math>T \! \left( x \right)</math>: | Таким образом, <math>3 \left< \rho^{2} \right> = 1</math>. Воспользуемся этим соотношением для построения требуемого профиля <math>T \! \left( x \right)</math>: | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 263: | Строка 215: | ||
\end{equation} | \end{equation} | ||
</math> | </math> | ||
+ | где <math>T_0</math>~--- значение температуры при <math>t = 0</math>. | ||
− | + | Перепишем выражение~(\ref{vTemp2Average}): | |
− | |||
− | Перепишем выражение | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 275: | Строка 225: | ||
Выразим <math>v^{2}\!\left(x\right)</math> и <math>v\!\left(x\right)</math>: | Выразим <math>v^{2}\!\left(x\right)</math> и <math>v\!\left(x\right)</math>: | ||
− | |||
− | |||
\begin{equation} | \begin{equation} | ||
v^{2}\!\left(x\right) = \frac{ 3 k T_0}{m} \rho^{2}\!\left(x\right), | v^{2}\!\left(x\right) = \frac{ 3 k T_0}{m} \rho^{2}\!\left(x\right), | ||
\end{equation} | \end{equation} | ||
</math> | </math> | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 287: | Строка 234: | ||
\label{vTemp} | \label{vTemp} | ||
\end{equation} | \end{equation} | ||
− | < | + | <math> |
Введём величину дисперсии скоростей <math>\sigma_{v}^{2}</math> так, чтобы при <math>\sigma_{v}^{2} = 1</math> в начальный момент времени механическая энергия волны равнялась тепловой энергии цепочки: | Введём величину дисперсии скоростей <math>\sigma_{v}^{2}</math> так, чтобы при <math>\sigma_{v}^{2} = 1</math> в начальный момент времени механическая энергия волны равнялась тепловой энергии цепочки: | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 300: | Строка 246: | ||
Тогда значение <math>\sigma_{v}^{2} = 10</math> можно трактовать следующим образом: при <math>t~=~0</math> энергия теплового шума в кристалле в <math>10</math> раз превосходит механическую энергию волны. | Тогда значение <math>\sigma_{v}^{2} = 10</math> можно трактовать следующим образом: при <math>t~=~0</math> энергия теплового шума в кристалле в <math>10</math> раз превосходит механическую энергию волны. | ||
− | + | Перепишем~(\ref{vTemp}) с учётом соотношения~(\ref{initialTemp}): | |
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 309: | Строка 254: | ||
</math> | </math> | ||
− | Таким образом, были получены выражения для обеих компонент скорости: скорость механической волны задаётся | + | Таким образом, были получены выражения для обеих компонент скорости: скорость механической волны задаётся функцией~(\ref{initialMechanicalVelocity1}), тепловой шум определяется функцией (\ref{initialThermalVelocity}). Итоговое условие на начальные скорости в кристалле принимает следующий вид: |
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 317: | Строка 261: | ||
</math> | </math> | ||
− | Начальные перемещения задаются | + | Начальные перемещения при рассмотрении бегущей волны задаются формулой~(\ref{initialMechanicalDisplacement1}), при рассмотрении стоячей~--- формулой~(\ref{initialMechanicalDisplacement2}) |
=== Вычисление механической энергии === | === Вычисление механической энергии === | ||
− | + | В начальный момент скорость имеет распределение, близкое к синусоидальному. С течением времени цепочка теряет эту форму (в силу нелийненого взаимодействия между частицами), а механическая энергия переходит в тепловую. | |
− | В начальный момент скорость имеет распределение, близкое к синусоидальному. С течением времени цепочка теряет эту форму (в силу | ||
Механическая энергия системы <math>E\!\left(x\right)</math> складывается из кинетической энергии <math>K\!\left(x\right)</math> и потенциальной энергии <math>\Pi\!\left(x\right)</math> следующий образом: | Механическая энергия системы <math>E\!\left(x\right)</math> складывается из кинетической энергии <math>K\!\left(x\right)</math> и потенциальной энергии <math>\Pi\!\left(x\right)</math> следующий образом: | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 333: | Строка 275: | ||
\end{equation} | \end{equation} | ||
</math> | </math> | ||
− | + | где <math>\varepsilon\!\left(x\right)</math>~--- возникающие деформации. | |
− | где <math>\varepsilon\!\left(x\right)</math> | ||
Так как нас интересует энергия не всей системы (включающей тепловой шум), а только заданной механической волны, будем искать энергию первой синусоидальной моды. Перепишем выражения для энергий: | Так как нас интересует энергия не всей системы (включающей тепловой шум), а только заданной механической волны, будем искать энергию первой синусоидальной моды. Перепишем выражения для энергий: | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 343: | Строка 283: | ||
\end{equation} | \end{equation} | ||
</math> | </math> | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 349: | Строка 288: | ||
\end{equation} | \end{equation} | ||
</math> | </math> | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 356: | Строка 294: | ||
\end{equation} | \end{equation} | ||
</math> | </math> | ||
− | |||
Первая синусоидальная мода также возникает и в тепловом шуме. Математическое ожидание от её энергии вычислим, усреднив по ансамблю реализаций систему без заданной волны. | Первая синусоидальная мода также возникает и в тепловом шуме. Математическое ожидание от её энергии вычислим, усреднив по ансамблю реализаций систему без заданной волны. | ||
Определим значения <math>K\!\left(x\right)</math> и <math>\Pi\!\left(x\right)</math> при <math>t = 0</math>: | Определим значения <math>K\!\left(x\right)</math> и <math>\Pi\!\left(x\right)</math> при <math>t = 0</math>: | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 366: | Строка 302: | ||
\end{equation} | \end{equation} | ||
</math> | </math> | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 374: | Строка 309: | ||
</math> | </math> | ||
− | Интеграл в | + | Интеграл в формуле (\ref{initialPotentialEnergy}) равен нулю в случае стоячей волны (в начальный момент времени нет деформаций). Покажем, что этот интеграл также будет равен нулю в случае бегущей волны: |
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 387: | Строка 321: | ||
То есть, в начальный момент времени потенциальная энергия равна нулю, а механическая энергия задаваемой волны равна | То есть, в начальный момент времени потенциальная энергия равна нулю, а механическая энергия задаваемой волны равна | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 395: | Строка 328: | ||
</math> | </math> | ||
− | + | Значение (\ref{initialMechanicalEnergy}) выберем в качестве масштаба энергии. | |
− | |||
− | |||
+ | Чтобы определить закон, по которому происходит переход механической энергии в тепловую, введём безразмерный параметр <math>E^*\!\left( x \right) = \frac{E\!\left( x \right)}{E\!\left( x \right)|_{t=0}}</math>~--- нормированное значение механической энергии цепочки: | ||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 410: | Строка 342: | ||
=== Численное решение дифференциального уравнения === | === Численное решение дифференциального уравнения === | ||
− | Для решения системы уравнений движения | + | Для решения системы уравнений движения~(\ref{chainEquations}) воспользуемся методом центральных разностей. Перепишем граничные условия~(\ref{boundaryConditions}): |
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 417: | Строка 348: | ||
\end{equation} | \end{equation} | ||
</math> | </math> | ||
+ | где <math>N \gg 1</math>~--- число независимых частиц. | ||
− | + | Вычисление правой части~(\ref{chainEquations}) аналогично решению с помощью разностной схемы. Рассчитанные таким образом ускорения частиц пересчитываем в скорости: | |
− | |||
− | Вычисление правой части | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 428: | Строка 357: | ||
</math> | </math> | ||
− | Затем скорости пересчитываем в смещения частиц. Важно отметить, что по данному методу расчёт скоростей производится на промежуточных шагах по времени относительно расчёта смещений. | + | Затем скорости пересчитываем в смещения частиц. Важно отметить, что по данному методу расчёт скоростей производится на промежуточных шагах по времени относительно расчёта смещений. |
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 437: | Строка 365: | ||
Значение <math>\varepsilon\!\left(x\right)</math> вычисляем следующим образом: | Значение <math>\varepsilon\!\left(x\right)</math> вычисляем следующим образом: | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 448: | Строка 375: | ||
Для достоверного численного решения дифференциального уравнения относительно <math>\left(u\!\left(x\right), \dot u\!\left(x\right)\right)</math> длина цепочки должна стремиться к бесконечной, чтобы соответствовать модели идеального кристалла. При проведении численных экспериментов удовлетворительная точность расчётов достигалась на цепочках из 10 миллионов частиц и более. | Для достоверного численного решения дифференциального уравнения относительно <math>\left(u\!\left(x\right), \dot u\!\left(x\right)\right)</math> длина цепочки должна стремиться к бесконечной, чтобы соответствовать модели идеального кристалла. При проведении численных экспериментов удовлетворительная точность расчётов достигалась на цепочках из 10 миллионов частиц и более. | ||
− | Статистические характеристики системы также можно вычислить путём усреднения бесконечного множества реализаций (каждая реализация | + | Статистические характеристики системы также можно вычислить путём усреднения бесконечного множества реализаций (каждая реализация~--- задача Коши для цепочки). Благодаря этому подходу, предложенному в работах \cite{art1}, \cite{art2}, достигается эффективность проводимых численных экспериментов. |
− | Результаты, приводимые далее, получены усреднением <math>5 \cdot 10^4</math> реализаций цепочки c <math>N = 10^3</math> (за исключением <math>\rho</math> все параметры в рамках серии сохранялись). Такая конфигурация будет аналогична модели кристалла с 50 миллионами частиц, но допускает параллельный компьютерный расчёт реализаций. Проведение параллельного моделирования коротких цепочек на суперкомпьютерных системах привело к | + | Результаты, приводимые далее, получены усреднением <math>5 \cdot 10^4</math> реализаций цепочки c <math>N = 10^3</math> (за исключением <math>\rho</math> все параметры в рамках серии сохранялись). Такая конфигурация будет аналогична модели кристалла с 50 миллионами частиц, но допускает параллельный компьютерный расчёт реализаций. Проведение параллельного моделирования коротких цепочек на суперкомпьютерных системах привело к сушественному сокращению общего времени вычислений. |
По мере повышения величины дисперсии <math>\sigma_v^2</math> шаг по времени <math>\Delta t</math> следует уменьшать. Иначе влияние нелинейной компоненты силы приводит к разрушению системы, эта неустойчивость связана с нарушением условия Куранта. Однако, уменьшение <math>\Delta t</math> увеличивает время, требуемое для проведения численного моделирования. Так как численно <math>\sigma_v^2</math> соответствует отношению тепловой энергии в кристалле к энергии механической волны, её значение в реальных системах на несколько порядков превышает исследуемое в рамках данной работы. | По мере повышения величины дисперсии <math>\sigma_v^2</math> шаг по времени <math>\Delta t</math> следует уменьшать. Иначе влияние нелинейной компоненты силы приводит к разрушению системы, эта неустойчивость связана с нарушением условия Куранта. Однако, уменьшение <math>\Delta t</math> увеличивает время, требуемое для проведения численного моделирования. Так как численно <math>\sigma_v^2</math> соответствует отношению тепловой энергии в кристалле к энергии механической волны, её значение в реальных системах на несколько порядков превышает исследуемое в рамках данной работы. | ||
В качестве масштаба времени возьмём период колебаний материальной точки массы <math>m</math> на пружине жёсткости <math>C_1</math>: | В качестве масштаба времени возьмём период колебаний материальной точки массы <math>m</math> на пружине жёсткости <math>C_1</math>: | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
\quad T_0~\mathrel{\stackrel{\rm def}=} \frac{2 \pi}{\omega_0} = 2 \pi \sqrt{\frac{m}{C_1}}. | \quad T_0~\mathrel{\stackrel{\rm def}=} \frac{2 \pi}{\omega_0} = 2 \pi \sqrt{\frac{m}{C_1}}. | ||
\end{equation} | \end{equation} | ||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
</math> | </math> | ||
− | + | \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}. Моделируемое время ограничено: <math>0 \leq t \leq 7 \cdot 10^4~T_0</math>. | |
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | <math> | ||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
== Заключение == | == Заключение == | ||
Строка 765: | Строка 439: | ||
Показано, что энергия стоячей волны убывает в 4 раза быстрее энергии бегущей волны при одних параметрах эксперимента. Продемонстрированы предельные значения <math>\beta</math>: | Показано, что энергия стоячей волны убывает в 4 раза быстрее энергии бегущей волны при одних параметрах эксперимента. Продемонстрированы предельные значения <math>\beta</math>: | ||
− | |||
* <math>\beta = 2</math> при <math>\sigma_v^2 = 0</math>; | * <math>\beta = 2</math> при <math>\sigma_v^2 = 0</math>; | ||
* <math>\beta \rightarrow 1</math> при <math>\sigma_v^2 \rightarrow \infty</math>. | * <math>\beta \rightarrow 1</math> при <math>\sigma_v^2 \rightarrow \infty</math>. | ||
Строка 773: | Строка 446: | ||
Результаты работы были получены с использованием вычислительных ресурсов суперкомпьютерного центра Санкт-Петербургского политехнического университета Петра Великого. Каждый расчёт производился на отдельном ядре 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> – заданный масштаб времени). | Результаты работы были получены с использованием вычислительных ресурсов суперкомпьютерного центра Санкт-Петербургского политехнического университета Петра Великого. Каждый расчёт производился на отдельном ядре 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> – заданный масштаб времени). | ||
− | Автор благодарен | + | Автор благодарен Д. В. Цветкову за полезные обсуждения. |
== Список литературы == | == Список литературы == | ||
+ | <references> |