Редактирование: Преобразование механической энергии в тепловую в одномерном кристалле
Внимание! Вы не авторизовались на сайте. Ваш 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-го века описывать взаимодействия на наноуровне. |
В данной работе рассматривается одна из нелинейных моделей кристалла, описанных Энрико Ферми. Ключевым отличием является введение тепловой энергии через задание дисперсии начальных скоростей. | В данной работе рассматривается одна из нелинейных моделей кристалла, описанных Энрико Ферми. Ключевым отличием является введение тепловой энергии через задание дисперсии начальных скоростей. | ||
Строка 31: | Строка 24: | ||
Рассмотрим модель одномерного кристалла: цепочку одинаковых частиц массы <math>m</math>, соединённых одинаковыми нелинейными пружинами (квадратичная нелинейность). | Рассмотрим модель одномерного кристалла: цепочку одинаковых частиц массы <math>m</math>, соединённых одинаковыми нелинейными пружинами (квадратичная нелинейность). | ||
− | Каждой частице присвоим свой целочисленный индекс. Принимая <math>u_k</math> за перемещение <math>k</math>-й частицы, а <math>C_1</math> и <math>C_2</math> | + | Каждой частице присвоим свой целочисленный индекс. Принимая <math>u_k</math> за перемещение <math>k</math>-й частицы, а <math>C_1</math> и <math>C_2</math>~--- за коэффициенты при линейном и квадратичном членах в разложении силы взаимодействия $k$-й частицы с двумя ближайшими соседями, получим следующее дифференциальное уравнение: |
<math> | <math> | ||
Строка 39: | Строка 32: | ||
\end{equation} | \end{equation} | ||
</math> | </math> | ||
− | + | где <math>\dot f</math>~--- производная <math>f</math> по времени | |
− | где <math>\dot f</math> | ||
Перепишем это выражение, раскрыв скобки и введя новые обозначения для коэффициентов при линейной и нелинейной части: | Перепишем это выражение, раскрыв скобки и введя новые обозначения для коэффициентов при линейной и нелинейной части: | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 49: | Строка 40: | ||
\end{equation} | \end{equation} | ||
</math> | </math> | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 55: | Строка 45: | ||
\end{equation} | \end{equation} | ||
</math> | </math> | ||
− | |||
Здесь <math>\omega_0</math> соответствует циклической частоте колебаний. | Здесь <math>\omega_0</math> соответствует циклической частоте колебаний. | ||
Заметим, что | Заметим, что | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 72: | Строка 60: | ||
\end{equation} | \end{equation} | ||
</math> | </math> | ||
− | + | Воспользуемся заменой~(\ref{replace}), и тогда уравнения динамики цепочки примут следующий вид: | |
− | Воспользуемся | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 82: | Строка 68: | ||
</math> | </math> | ||
− | Аналогичное уравнение использовалось в оригинальной работе Ферми, Паста и Улама | + | Аналогичное уравнение использовалось в оригинальной работе Ферми, Паста и Улама\cite{fermi1965crystal}. |
Механическую энергию системы в начальный момент времени (то есть, при <math>t = 0</math>) будем задавать с помощью синусоидальной волны. Из-за теплового движения частиц волна будет терять свою форму, а её энергия перейдёт в тепловую энергию системы. Для описания этого процесса численно решим уравнения динамики кристалла. | Механическую энергию системы в начальный момент времени (то есть, при <math>t = 0</math>) будем задавать с помощью синусоидальной волны. Из-за теплового движения частиц волна будет терять свою форму, а её энергия перейдёт в тепловую энергию системы. Для описания этого процесса численно решим уравнения динамики кристалла. | ||
− | Чтобы замкнуть систему из <math>k</math> уравнений | + | Чтобы замкнуть систему из <math>k</math> уравнений~(\ref{chainEquations}), определим граничные условия, а также по <math>k</math> начальных условий на перемещения и скорости. Для задания начальных скоростей воспользуемся генератором случайных чисел с равномерным распределением. При этом дисперсию скоростей будем задавать в соответствии с требуемым температурным профилем. |
=== Начальные и граничные условия === | === Начальные и граничные условия === | ||
− | Чтобы соответствовать идеальной модели, исследуемая цепочка должна иметь бесконечную длину. Для проведения численного моделирования будем рассматривать кристалл конечной длины с заданным условием периодичности (граничные условия Борна-Кармана): | + | Чтобы соответствовать идеальной модели, исследуемая цепочка должна иметь бесконечную длину. Для проведения численного моделирования будем рассматривать кристалл конечной длины с заданным условием периодичности (граничные условия Борна~--- Кармана): |
<math> | <math> | ||
Строка 98: | Строка 84: | ||
\end{equation} | \end{equation} | ||
</math> | </math> | ||
− | + | где <math>L</math>~--- длина кристалла, <math>a</math>~--- равновесное расстояние между частицами. | |
− | где <math>L</math> | ||
Таким образом, условие теплоизоляции выполняется. | Таким образом, условие теплоизоляции выполняется. | ||
Строка 108: | Строка 93: | ||
==== Начальные условия для бегущей волны ==== | ==== Начальные условия для бегущей волны ==== | ||
− | |||
Рассмотрим уравнение бегущей волны (пусть волна распространяется в одном направлении): | Рассмотрим уравнение бегущей волны (пусть волна распространяется в одном направлении): | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 117: | Строка 100: | ||
\end{equation} | \end{equation} | ||
</math> | </math> | ||
+ | где <math>A</math>~--- амплитуда колебаний, <math>k_0</math>~--- волновое число, <math>\phi</math>~--- фаза колебаний. | ||
− | + | Величину <math>u\!\left(x\right)</math> на момент времени <math>t</math> получим интегрированием уравнения~(\ref{runningWaveEquation1}): | |
− | |||
− | Величину <math>u\!\left(x\right)</math> на момент времени <math>t</math> получим интегрированием уравнения | ||
<math> | <math> | ||
Строка 130: | Строка 112: | ||
Пусть длина волны совпадает с длиной кристалла, тогда | Пусть длина волны совпадает с длиной кристалла, тогда | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 136: | Строка 117: | ||
\end{equation} | \end{equation} | ||
</math> | </math> | ||
+ | где <math>\lambda</math>~--- длина волны. | ||
− | + | Фазу колебаний <math>\phi</math> положим равной нулю и запишем уравнения~(\ref{runningWaveEquation1}) и~(\ref{runningWaveEquation2}) для начального момента времени (<math>t=0</math>): | |
− | |||
− | Фазу колебаний <math>\phi</math> положим равной нулю и запишем уравнения | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 147: | Строка 126: | ||
\end{equation} | \end{equation} | ||
</math> | </math> | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 155: | Строка 133: | ||
</math> | </math> | ||
− | На рисунке | + | На рисунке~\ref{runningWaveInitialConditions} приведён эскиз начальных перемещений и скоростей частиц кристалла при заданной бегущей волне. |
− | [ | + | <math> |
+ | \begin{figure}[ht!] | ||
+ | \center{\input{Data/InitialConditions/rw.tex}} | ||
+ | \caption{Начальные условия для бегущей волны} | ||
+ | \label{runningWaveInitialConditions} | ||
+ | \end{figure} | ||
+ | </math> | ||
==== Начальные условия для стоячей волны ==== | ==== Начальные условия для стоячей волны ==== | ||
− | + | Начальные условия~(\ref{initialMechanicalVelocity1}) и~(\ref{initialMechanicalDisplacement1}) позволяют задать бегущую волну. Получим подобные условия также для стоячей волны. Для этого рассмотрим две волны амплитуды <math>\frac{A}{2}</math>, бегущие в противоположных направлениях с разностью фаз <math>\pi</math>: | |
− | Начальные условия | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 170: | Строка 152: | ||
Результатом наложения этих двух волн будет стоячая волна, описываемая следующим уравнением: | Результатом наложения этих двух волн будет стоячая волна, описываемая следующим уравнением: | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 178: | Строка 159: | ||
</math> | </math> | ||
− | Аналогично, проинтегрируем | + | Аналогично, проинтегрируем уравнение~(\ref{standingWaveEquation1}) по времени <math>t</math>: |
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 187: | Строка 167: | ||
</math> | </math> | ||
− | Подставляя в | + | Подставляя в~(\ref{standingWaveEquation1}) и~(\ref{standingWaveEquation2}) <math>k_0 = \frac{2 \pi}{L}</math> и <math>t = 0</math>, получим начальные условия для стоячей волны: |
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 195: | Строка 174: | ||
\end{equation} | \end{equation} | ||
</math> | </math> | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 203: | Строка 181: | ||
</math> | </math> | ||
− | Соответствующий график начальных перемещений и скоростей частиц цепочки приведён на рисунке | + | Соответствующий график начальных перемещений и скоростей частиц цепочки приведён на рисунке~\ref{standingWaveInitialConditions}. |
− | [ | + | <math> |
+ | \begin{figure}[ht!] | ||
+ | \center{\input{Data/InitialConditions/sw.tex}} | ||
+ | \caption{Начальные условия для стоячей волны} | ||
+ | \label{standingWaveInitialConditions} | ||
+ | \end{figure} | ||
+ | </math> | ||
− | Функции скоростей в начальный момент времени для бегущей и стоячей волн совпадают, различаются только начальные перемещения. | + | Функции скоростей в начальный момент времени для бегущей~(\ref{initialMechanicalVelocity1}) и стоячей~(\ref{initialMechanicalVelocity2}) волн совпадают, различаются только начальные перемещения. |
==== Начальные условия теплового движения ==== | ==== Начальные условия теплового движения ==== | ||
Введём оператор усреднения функции по её аргументу: | Введём оператор усреднения функции по её аргументу: | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 218: | Строка 201: | ||
\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: | Строка 217: | ||
\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: | Строка 236: | ||
Таким образом, <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: | Строка 241: | ||
\end{equation} | \end{equation} | ||
</math> | </math> | ||
+ | где <math>T_0</math>~--- значение температуры при <math>t = 0</math>. | ||
− | + | Перепишем выражение~(\ref{vTemp2Average}): | |
− | |||
− | Перепишем выражение | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 275: | Строка 251: | ||
Выразим <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>: | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 281: | Строка 256: | ||
\end{equation} | \end{equation} | ||
</math> | </math> | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 290: | Строка 264: | ||
Введём величину дисперсии скоростей <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: | Строка 273: | ||
Тогда значение <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: | Строка 281: | ||
</math> | </math> | ||
− | Таким образом, были получены выражения для обеих компонент скорости: скорость механической волны задаётся | + | Таким образом, были получены выражения для обеих компонент скорости: скорость механической волны задаётся функцией~(\ref{initialMechanicalVelocity1}), тепловой шум определяется функцией (\ref{initialThermalVelocity}). Итоговое условие на начальные скорости в кристалле принимает следующий вид: |
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 317: | Строка 288: | ||
</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: | Строка 302: | ||
\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: | Строка 310: | ||
\end{equation} | \end{equation} | ||
</math> | </math> | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 349: | Строка 315: | ||
\end{equation} | \end{equation} | ||
</math> | </math> | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 356: | Строка 321: | ||
\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: | Строка 329: | ||
\end{equation} | \end{equation} | ||
</math> | </math> | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 374: | Строка 336: | ||
</math> | </math> | ||
− | Интеграл в | + | Интеграл в формуле (\ref{initialPotentialEnergy}) равен нулю в случае стоячей волны (в начальный момент времени нет деформаций). Покажем, что этот интеграл также будет равен нулю в случае бегущей волны: |
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 387: | Строка 348: | ||
То есть, в начальный момент времени потенциальная энергия равна нулю, а механическая энергия задаваемой волны равна | То есть, в начальный момент времени потенциальная энергия равна нулю, а механическая энергия задаваемой волны равна | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 395: | Строка 355: | ||
</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: | Строка 369: | ||
=== Численное решение дифференциального уравнения === | === Численное решение дифференциального уравнения === | ||
− | Для решения системы уравнений движения | + | Для решения системы уравнений движения~(\ref{chainEquations}) воспользуемся методом центральных разностей. Перепишем граничные условия~(\ref{boundaryConditions}): |
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 417: | Строка 375: | ||
\end{equation} | \end{equation} | ||
</math> | </math> | ||
+ | где <math>N \gg 1</math>~--- число независимых частиц. | ||
− | + | Вычисление правой части~(\ref{chainEquations}) аналогично решению с помощью разностной схемы. Рассчитанные таким образом ускорения частиц пересчитываем в скорости: | |
− | |||
− | Вычисление правой части | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 428: | Строка 384: | ||
</math> | </math> | ||
− | Затем скорости пересчитываем в смещения частиц. Важно отметить, что по данному методу расчёт скоростей производится на промежуточных шагах по времени относительно расчёта смещений. | + | Затем скорости пересчитываем в смещения частиц. Важно отметить, что по данному методу расчёт скоростей производится на промежуточных шагах по времени относительно расчёта смещений. |
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 437: | Строка 392: | ||
Значение <math>\varepsilon\!\left(x\right)</math> вычисляем следующим образом: | Значение <math>\varepsilon\!\left(x\right)</math> вычисляем следующим образом: | ||
− | |||
<math> | <math> | ||
\begin{equation} | \begin{equation} | ||
Строка 448: | Строка 402: | ||
Для достоверного численного решения дифференциального уравнения относительно <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> | ||
− | + | Параметры при моделировании стоячей волны | |
{| class="wikitable" | {| class="wikitable" | ||
|- | |- | ||
! <math>\sigma_v^2</math> !! <math>\Delta \sigma_v^2</math> !! <math>\Delta t</math> !! Постоянные параметры | ! <math>\sigma_v^2</math> !! <math>\Delta \sigma_v^2</math> !! <math>\Delta t</math> !! Постоянные параметры | ||
|- | |- | ||
− | | <math>[0; 14]</math> || <math>1</math> || <math>0.05~T_0</math> || <math>\omega_0 = 1, \quad \alpha_0 = 1, \quad A = 0. | + | | <math>[0; 14]</math> || <math>1</math> || <math>0.05~T_0</math> || <math>\omega_0 = 1, \quad \alpha_0 = 1, \quad A = 0.01</math> |
|- | |- | ||
− | | <math>[15; 60]</math> || <math>5</math> || <math>0.05~T_0</math> || <math>\omega_0 = 1, \quad \alpha_0 = 1, \quad A = 0. | + | | <math>[15; 60]</math> || <math>5</math> || <math>0.05~T_0</math> || <math>\omega_0 = 1, \quad \alpha_0 = 1, \quad A = 0.01</math> |
|- | |- | ||
− | | <math>[65; 75]</math> || <math>5</math> || <math>0.005~T_0</math> || <math>\omega_0 = 1, \quad \alpha_0 = 1, \quad A = 0. | + | | <math>[65; 75]</math> || <math>5</math> || <math>0.005~T_0</math> || <math>\omega_0 = 1, \quad \alpha_0 = 1, \quad A = 0.01</math> |
|- | |- | ||
− | | <math>[80; 100]</math> || <math>5</math> || <math>0.005~T_0</math> || <math>\omega_0 = 1, \quad \alpha_0 = 1, \quad A = 0. | + | | <math>[80; 100]</math> || <math>5</math> || <math>0.005~T_0</math> || <math>\omega_0 = 1, \quad \alpha_0 = 1, \quad A = 0.01</math> |
|} | |} | ||
− | + | Параметры при моделировании бегущей волны | |
{| class="wikitable" | {| class="wikitable" | ||
|- | |- | ||
Строка 488: | Строка 442: | ||
|} | |} | ||
− | Значения параметров проводимых экспериментов приведены в таблицах | + | Значения параметров проводимых экспериментов приведены в таблицах~\ref{standingWaveParameters} и~\ref{runningWaveParameters}. Моделируемое время ограничено: <math>0 \leq t \leq 7 \cdot 10^4~T_0</math>. |
== Результаты == | == Результаты == | ||
Строка 494: | Строка 448: | ||
=== Преобразование механической энергии стоячей волны === | === Преобразование механической энергии стоячей волны === | ||
− | На рисунке | + | На рисунке~\ref{standingWave} приведены вычисленные значения нормированной механической энергии стоячей волны $E^*\!\left(\tau\right)$ для нескольких значений $\sigma_v^2$ (здесь и далее $\tau = \frac{t}{T_0}$ – безразмерное время). |
− | [ | + | \begin{figure}[!htb] |
+ | \center{\input{Data/Energy/sw.tex}} | ||
+ | \caption{Изменение нормированной механической энергии стоячей волны в зависимости от $t$} | ||
+ | \label{standingWave} | ||
+ | \end{figure} | ||
− | Для одномерного кристалла при отсутствии теплового шума ( | + | Для одномерного кристалла при отсутствии теплового шума ($\sigma_v^2 = 0$) была установлена \cite{presTsvetkovAPM16} следующая зависимость для механической энергии волны: |
− | |||
\begin{equation} | \begin{equation} | ||
E^* \sim e^{-\tau^2}. | E^* \sim e^{-\tau^2}. | ||
\label{expectationTau2} | \label{expectationTau2} | ||
\end{equation} | \end{equation} | ||
− | |||
− | Подобная зависимость не выполняется при ненулевой дисперсии скоростей частиц цепочки: механическая энергия убывает быстрее при бóльших значениях | + | Подобная зависимость не выполняется при ненулевой дисперсии скоростей частиц цепочки: механическая энергия убывает быстрее при бóльших значениях $\sigma_v^2$. Предположим, что $E^*$ уменьшается по модифицированному экспоненциальному закону и показатель степени при $\tau$ в общем случае неизвестен: |
− | |||
\begin{equation} | \begin{equation} | ||
E^* = e^{-\alpha \tau^\beta}, | E^* = e^{-\alpha \tau^\beta}, | ||
\label{expectationTauBeta} | \label{expectationTauBeta} | ||
\end{equation} | \end{equation} | ||
− | + | где $\beta$~--- функция от $\sigma_v^2$, $\alpha$~--- функция от $\sigma_v^2$ и $\lambda$. | |
− | + | В рамках гипотезы~(\ref{expectationTauBeta}) $\beta$ имеет определённое значение для каждого значения $\sigma_v^2$. При $\sigma_v^2 = 0$ показатель $\beta = 2$~\cite{presTsvetkovAPM16}, по мере же увеличения $\sigma_v^2$, согласно рисунку~\ref{standingWave}, величина $\beta$ должна уменьшаться. | |
− | + | Если выражение~(\ref{expectationTauBeta}) верно, то показатель $\beta$ в зависимости от величины теплового шума можно определить, перестроив график энергии~\ref{standingWave} в логарифмических осях. | |
− | |||
− | |||
\begin{equation} | \begin{equation} | ||
\begin{split} | \begin{split} | ||
Строка 533: | Строка 486: | ||
\label{expectationTauBetaLine} | \label{expectationTauBetaLine} | ||
\end{equation} | \end{equation} | ||
− | |||
− | Обозначим | + | Обозначим $\ln\!\left(\tau\right)$ за $x$, а $\ln\!\left(\ln\!\left(\frac{1}{E^*}\right)\right)$ за $y$, тогда уравнение~(\ref{expectationTauBetaLine}) сводится к линейному. Угол наклона прямой в осях $\left(\ln\!\left(\tau\right); \ln\!\left(\ln\!\left(\frac{1}{E^*}\right)\right) \right)$, выраженный в радианах, будет численно равен $\beta$ для определённого значения дисперсии. |
− | |||
− | |||
\begin{equation} | \begin{equation} | ||
y\!\left(x\right) = \ln\!\left(\alpha\right) + \beta x | y\!\left(x\right) = \ln\!\left(\alpha\right) + \beta x | ||
\label{lineEquation} | \label{lineEquation} | ||
\end{equation} | \end{equation} | ||
− | |||
− | [ | + | \begin{figure}[!htb] |
+ | \center{\input{Data/Energy/swLine.tex}} | ||
+ | \caption{Изменение механической энергии стоячей волны с течением времени в логарифмических осях} | ||
+ | \label{standingWaveLine} | ||
+ | \end{figure} | ||
− | Соответствующее преобразование графика энергии приведено на рисунке | + | Соответствующее~(\ref{expectationTauBetaLine}) преобразование графика энергии приведено на рисунке~\ref{standingWaveLine}. Получившиеся функции действительно аппроксимbруются наклонными прямыми, что подтверждает высказанную гипотезу~(\ref{expectationTauBeta}). |
− | [ | + | \begin{figure}[!htb] |
+ | \center{\input{Data/Beta/swNoLine.tex}} | ||
+ | \caption{Значения $\beta$ при разных значениях дисперсии, стоячая волна} | ||
+ | \label{standingWaveBeta} | ||
+ | \end{figure} | ||
− | Найденные значения | + | Найденные значения $\beta$ приведены на рисунке~\ref{standingWaveBeta}. Тепловая энергия в реальном материале может превосходить механическую на несколько порядков. Чтобы оценить значение $\beta$ в таком случае, построим аппроксимацию этого графика и определим предельное значение показателя при условии $\sigma_v^2 \rightarrow \infty$. |
− | + | Перестроим график из осей $\left(\sigma_{v}^{2} ; \beta \right)$ в $\left(\ln\!\left(\sigma_{v}^{2}\right) ; \ln\!\left(\beta - 1\right) \right)$. Результат приведён на рисунке~\ref{standingWaveBetaLog}. Начиная с $x = 1.5$ значения могут быть аппроксимированы прямой: | |
− | |||
− | |||
− | |||
\begin{equation} | \begin{equation} | ||
y\!\left(x\right) = \ln\!\left(1.08\right) - 0.25 x, | y\!\left(x\right) = \ln\!\left(1.08\right) - 0.25 x, | ||
\label{standingWaveApproximation} | \label{standingWaveApproximation} | ||
\end{equation} | \end{equation} | ||
− | + | где $x = \ln\!\left(\sigma_{v}^{2}\right)$, $y = \ln\!\left(\beta\right)$. | |
− | + | Выражая из~(\ref{standingWaveApproximation}) функцию $\beta$, получаем следующее: | |
− | |||
− | |||
− | |||
\begin{equation} | \begin{equation} | ||
\beta \! \left(\sigma_{v}^{2} \right) = 1 + \frac{1.08}{\sqrt[4]{\sigma_{v}^{2}}} = 1 + \frac{1.08}{\sqrt{\sigma_{v}}}, | \beta \! \left(\sigma_{v}^{2} \right) = 1 + \frac{1.08}{\sqrt[4]{\sigma_{v}^{2}}} = 1 + \frac{1.08}{\sqrt{\sigma_{v}}}, | ||
\label{standingWaveApproximation2} | \label{standingWaveApproximation2} | ||
\end{equation} | \end{equation} | ||
− | + | где $\sigma_{v}$~--- среднеквадратическое отклонение скоростей частиц, отнесённое к своему начальному значению (при $\tau~=~0$). | |
− | |||
− | где | ||
− | [ | + | \begin{figure}[!htb] |
+ | \center{\input{Data/Beta/swLog.tex}} | ||
+ | \caption{Значения функции $\beta \! \left(\sigma_{v}^{2} \right)$ в логарифмических осях, стоячая волна} | ||
+ | \label{standingWaveBetaLog} | ||
+ | \end{figure} | ||
− | [ | + | \begin{figure}[!htb] |
+ | \begin{subfigure}{\linewidth} | ||
+ | \center{\input{Data/Beta/sw.tex}} | ||
+ | \caption{} | ||
+ | \end{subfigure} | ||
+ | \par\medskip | ||
+ | \begin{subfigure}{\linewidth} | ||
+ | \center{\input{Data/Beta/swLongLine.tex}} | ||
+ | \caption{} | ||
+ | \end{subfigure} | ||
+ | \caption{Аппроксимация графика $\beta$, стоячая волна} | ||
+ | \label{standingWaveBetaApproximation} | ||
+ | \end{figure} | ||
− | Аппроксимированный график | + | Аппроксимированный график $\beta$ от $\sigma_{v}^{2}$ приведён на рисунке~\ref{standingWaveBetaApproximation} (сверху). Найдём предельное значение $\beta$: |
− | |||
\begin{equation} | \begin{equation} | ||
\lim_{\sigma_{v}^{2} \to \infty} \beta \! \left(\sigma_{v}^{2} \right) = | \lim_{\sigma_{v}^{2} \to \infty} \beta \! \left(\sigma_{v}^{2} \right) = | ||
Строка 588: | Строка 552: | ||
\label{standingWaveLimit} | \label{standingWaveLimit} | ||
\end{equation} | \end{equation} | ||
− | |||
− | Можно считать, для больших значений | + | Можно считать, для больших значений $\sigma_{v}^{2}$ механическая энергия будет переходить в тепловую по закону $E^* = e^{-\alpha \tau}$. Сходимость $\beta$ к единице также видна на рисунке~\ref{standingWaveBetaApproximation} (снизу). Уже при $\sigma_{v} = 12 \cdot 10^3$ ошибка в определении $\beta$ единицей составит менее процента. |
− | [ | + | \begin{figure}[!htb] |
+ | \center{\input{Data/Alpha/swLog.tex}} | ||
+ | \caption{Аппроксимация графика $\alpha$, стоячая волна} | ||
+ | \label{standingWaveAlphaApproximation} | ||
+ | \end{figure} | ||
− | Исследуем поведение | + | Исследуем поведение $\alpha$ в зависимости от величины теплового шума. Для этого воспользуемся уравнениями прямых, полученных аппроксимацией механической энернигии $E^*$ в логарифмических осях (см. рисунок~\ref{standingWaveBetaLog}). Подставив в уравнение (\ref{lineEquation}) $x = 0$, определим значения $\ln\!\left(\alpha\right)$ как функции от $\sigma_{v}^{2}$. Результат приведён на рисунке~\ref{standingWaveAlphaApproximation}. Была получена следующая формула для аппроксимации $\alpha$: |
− | |||
\begin{equation} | \begin{equation} | ||
\ln\!\left(\alpha\right) = \ln\!\left(0.14\right) - \frac{50}{5 + \sqrt{\sigma_{v}^{2}}} = \ln\!\left(0.14\right) - \frac{50}{5 + \sigma_{v}}. | \ln\!\left(\alpha\right) = \ln\!\left(0.14\right) - \frac{50}{5 + \sqrt{\sigma_{v}^{2}}} = \ln\!\left(0.14\right) - \frac{50}{5 + \sigma_{v}}. | ||
\label{standingWaveAlpha} | \label{standingWaveAlpha} | ||
\end{equation} | \end{equation} | ||
− | |||
− | Аналогично выражению для | + | Аналогично выражению для $\beta$, уравнение \ref{standingWaveAlpha} имеет предельное значение при $\sigma_{v}^{2}~\to~\infty$: |
− | |||
\begin{equation} | \begin{equation} | ||
\lim_{\sigma_{v}^{2} \to \infty} \ln\!\left(\alpha\right) = \lim_{\sigma_{v}^{2} \to \infty} \left( \ln\!\left(0.14\right) - \frac{50}{5 + \sqrt{\sigma_{v}^{2}}} \right) = \ln\!\left(0.14\right). | \lim_{\sigma_{v}^{2} \to \infty} \ln\!\left(\alpha\right) = \lim_{\sigma_{v}^{2} \to \infty} \left( \ln\!\left(0.14\right) - \frac{50}{5 + \sqrt{\sigma_{v}^{2}}} \right) = \ln\!\left(0.14\right). | ||
\label{standingWaveAlphaLimit} | \label{standingWaveAlphaLimit} | ||
\end{equation} | \end{equation} | ||
− | |||
− | |||
− | |||
− | + | Отсюда следует, что | |
\begin{equation} | \begin{equation} | ||
\lim_{\sigma_{v}^{2} \to \infty} \alpha = 0.14. | \lim_{\sigma_{v}^{2} \to \infty} \alpha = 0.14. | ||
\label{standingWaveAlphaLimit2} | \label{standingWaveAlphaLimit2} | ||
\end{equation} | \end{equation} | ||
− | |||
− | При | + | При $\sigma_{v}^{2} = 12 \cdot 10^3$ значение $\alpha$ отличается от предельного примерно на $0.2\%$. |
− | Дополнительное рассмотрение цепочек с кратным числом частиц привело к следующему результату: уменьшение длины кристалла в | + | Дополнительное рассмотрение цепочек с кратным числом частиц привело к следующему результату: уменьшение длины кристалла в $2$ раза приводит к увеличению предельного значения $\alpha$ в $4$ раза. То есть, $\alpha \sim k^{2}$ ($k$~--- волновое число, $k \sim \frac{1}{\lambda}$). Таким образом, в реальных материалах изменение нормированной механической энергии $E^*$ может быть описано экспоненциальным законом: |
− | |||
− | |||
\begin{equation} | \begin{equation} | ||
E^*\!\left(\tau\right) = e^{-0.14 k_{*}^{2} \tau}, | E^*\!\left(\tau\right) = e^{-0.14 k_{*}^{2} \tau}, | ||
\label{energySolution} | \label{energySolution} | ||
\end{equation} | \end{equation} | ||
− | + | где $k_{*} = \frac{k}{k_{1000}}$, $k_{1000}$~--- значение $k$ при $N = 1000$. | |
− | |||
− | где | ||
Сравним полученную зависимость с решением уравнения продольных колебаний стержня с внутренним трением. В используемых обозначениях уравнение колебаний может быть записано таким образом: | Сравним полученную зависимость с решением уравнения продольных колебаний стержня с внутренним трением. В используемых обозначениях уравнение колебаний может быть записано таким образом: | ||
− | |||
\begin{equation} | \begin{equation} | ||
\frac{\partial^2 u}{\partial \tau^2} = \frac{E T_{0}^{2}}{\rho a^2} \frac{\partial^2 u}{\partial x^2} | \frac{\partial^2 u}{\partial \tau^2} = \frac{E T_{0}^{2}}{\rho a^2} \frac{\partial^2 u}{\partial x^2} | ||
Строка 642: | Строка 597: | ||
\label{oscillationEquation} | \label{oscillationEquation} | ||
\end{equation} | \end{equation} | ||
− | + | где $E$~--- модуль Юнга, $\rho$~--- плотность, $2 \gamma$~--- безразмерный коэффициент при диссипативном члене, $x$~--- безразмерная переменная расстояния. | |
− | + | Решение дифференциального уравнения~(\ref{oscillationEquation}) имеет вид: | |
− | |||
− | Решение дифференциального уравнения имеет вид: | ||
− | |||
− | |||
\begin{equation} | \begin{equation} | ||
u = \frac{A_{c}}{\omega_{c}} e^{-\gamma k_{c}^{2} \tau} \sin\!\left( \omega_{c} \tau \right) \sin\!\left( k_{c} x \right), | u = \frac{A_{c}}{\omega_{c}} e^{-\gamma k_{c}^{2} \tau} \sin\!\left( \omega_{c} \tau \right) \sin\!\left( k_{c} x \right), | ||
\label{oscillationEquationSolution} | \label{oscillationEquationSolution} | ||
\end{equation} | \end{equation} | ||
− | + | где $\omega_{c} = k_{c} \sqrt{\frac{E T_{0}^{2}}{\rho a^2} - \gamma}$, $k_{c} = \frac{a}{2} k$, $A_{c} = \frac{a}{T0} A$. | |
− | |||
− | где | ||
− | |||
− | |||
− | + | Данный результат был представлен в работе~\cite{aristovich2006oscillations} и хорошо согласуется с формулой~(\ref{standingWaveEquation2}), переписанной в безразмерных величинах: | |
\begin{equation} | \begin{equation} | ||
u = \frac{A_{c}}{\hat{\omega_{c}}} e^{- 0.14 k_{*}^{2} \tau} \sin\!\left( \hat{\omega_{c}} \tau \right) \sin\!\left( k_{*} x \right), | u = \frac{A_{c}}{\hat{\omega_{c}}} e^{- 0.14 k_{*}^{2} \tau} \sin\!\left( \hat{\omega_{c}} \tau \right) \sin\!\left( k_{*} x \right), | ||
\label{standingWaveEquation3} | \label{standingWaveEquation3} | ||
\end{equation} | \end{equation} | ||
− | + | где $\hat{\omega_{c}} = T_{0} \omega_{0}.$ | |
− | |||
− | где | ||
− | |||
− | |||
− | + | Сопоставим решение~(\ref{oscillationEquationSolution}) с~(\ref{standingWaveEquation3}) и определим значений $\gamma$: | |
\begin{equation} | \begin{equation} | ||
\gamma = 0.14 \left(\frac{k_{*}}{k_{c}} \right)^{2} = 0.14 \left( \frac{k}{k_{1000}} \frac{1}{k a} \right)^{2} = 0.14 \left(\frac{1000}{N} \frac{N}{2 \pi} \right)^{2} \approx 3550. | \gamma = 0.14 \left(\frac{k_{*}}{k_{c}} \right)^{2} = 0.14 \left( \frac{k}{k_{1000}} \frac{1}{k a} \right)^{2} = 0.14 \left(\frac{1000}{N} \frac{N}{2 \pi} \right)^{2} \approx 3550. | ||
\end{equation} | \end{equation} | ||
− | |||
− | Покажем необратимость преобразования механической энергии в тепловую. Парадокс Ферми-Паста-Улама продемонстрировал явление возврата формы и механической энергии волны | + | Покажем необратимость преобразования механической энергии в тепловую. Парадокс Ферми-Паста-Улама продемонстрировал явление возврата формы и механической энергии волны \cite{fermi1965crystal}. Обратный переход энергии цепочки в механическую энергию был описан для одномерного нелинейного кристалла без теплового движения \cite{tsvetkov2015heatIn1D}. Такой переход также был экспериментально получен в данной работе. Рассмотрим влияние теплового шума на этот процесс. |
− | + | На рисунке~\ref{returningEnergy} можно заметить дополнительный изгиб графика энергии $E^*\!\left(T\right)$ при малых значениях $\sigma_v^2$. Так, для $\sigma_v^2 = 0$ возврат механической энергии наблюдается примерно в диапазоне от $\tau=60$ до $\tau=75$. Для $\sigma_v^2 = 1$ этот диапазон уже: от $\tau=65$ до $\tau=70$, при этом прирост энергии значительно меньше. Для высоких значений $\sigma_v^2$ обратный переход энергии не удаётся продемонстрировать. Можно утверждать, что процесс преобразования механической энергии в тепловую в больших системах с высокими значениями отнощения тепловой энергии к механической (соответствует реальным системам) является необратимым. | |
− | + | \begin{figure}[!htb] | |
+ | \center{\input{Data/Energy/swLong.tex}} | ||
+ | \caption{Изменение энергии стоячей волны} | ||
+ | \label{returningEnergy} | ||
+ | \end{figure} | ||
=== Преобразование механической энергии бегущей волны === | === Преобразование механической энергии бегущей волны === | ||
− | Рассмотрим график убывания механической энергии бегущей волны с течением времени (рисунок | + | Рассмотрим график убывания механической энергии бегущей волны с течением времени (рисунок~\ref{runningWave}). Качественно он совпадает с графиком энергии стоячей волны (рисунок~\ref{standingWave}). При отсутствии теплового шума для энергии волны также выполняется равенство~(\ref{expectationTau2}). |
− | |||
− | |||
− | Чтобы найти зависимость | + | Чтобы найти зависимость $\beta\!\left(\sigma_v^2\right)$, аналогичным образом перестроим график энергии $E^{*}\!\left(\tau\right)$ в логарифмических осях (см. рисунок~\ref{runningWaveLine}) и определим угол наклона $\beta$. |
− | + | Вид $\beta$ приведён на рисунке~\ref{runningWaveBetaApproximation}. Аппроксимация также была получена переводом осей $\left(\sigma_{v}^{2} ; \beta \right)$ в $\left(\ln\!\left(\sigma_{v}^{2}\right) ; \ln\!\left(\beta\right)\right)$ (см. рисунок~\ref{runningWaveBetaLog}). | |
+ | \begin{figure}[!htb] | ||
+ | \center{\input{Data/Energy/rw.tex}} | ||
+ | \caption{Изменение энергии бегущей волны в зависимости от $\tau$} | ||
+ | \label{runningWave} | ||
+ | \end{figure} | ||
− | + | \begin{figure}[!htb] | |
+ | \center{\input{Data/Energy/rwLine.tex}} | ||
+ | \caption{Изменение механической энергии стоячей волны с течением времени в логарифмических осях} | ||
+ | \label{runningWaveLine} | ||
+ | \end{figure} | ||
− | + | Для $\beta$ был получен следующий вид уравнения аппроксимации: | |
− | |||
− | |||
− | |||
− | |||
− | |||
\begin{equation} | \begin{equation} | ||
\beta \! \left(\sigma_{v}^{2} \right) = 1 + \frac{1.55}{\sqrt[4]{\sigma_{v}^{2}}} = 1 + \frac{1.55}{\sqrt{\sigma_{v}}}. | \beta \! \left(\sigma_{v}^{2} \right) = 1 + \frac{1.55}{\sqrt[4]{\sigma_{v}^{2}}} = 1 + \frac{1.55}{\sqrt{\sigma_{v}}}. | ||
\label{runningWaveApproximation} | \label{runningWaveApproximation} | ||
\end{equation} | \end{equation} | ||
− | |||
− | |||
− | |||
− | + | При условии, что $\sigma_{v}^{2} \to \infty$, значение~(\ref{runningWaveApproximation}) стремится к единице, аналогично~(\ref{standingWaveLimit}) в случае стоячей волны. | |
− | + | Значения $\alpha$ для различных $\sigma_{v}^{2}$ также были получены (рисунок~\ref{runningWaveAlphaLog}). Предел аппроксимирующей функции~(\ref{standingWaveAlpha}) равен соответствующему значению~(\ref{standingWaveAlphaLimit2}) для стоячей волны. Таким образом, в реальных материалах формула~(\ref{energySolution}) для механической энергии выполняется вне зависимости от типа заданной механической волны. | |
− | На рисунке | + | На рисунке~\ref{runningAndStandinsWavesBeta} приведён график значений $\beta$ для стоячей волны, расстянутый по горизонтальной оси в 4 раза. В таком масштабе экспериментальные значения $\beta$ для бегущей и стоячей волн могут быть аппроксимированы одной и той же кривой (заданной уравнением~(\ref{runningWaveApproximation})). То есть, в случае бегущей волны было получено в 4 раза более медленное убывание механической энергии, чем в случае стоячей. Такая же пропорциональность наблюдается для значений $\alpha$ (см. рисунок~\ref{runningAndStandinsWavesAlpha}). |
− | |||
\begin{equation} | \begin{equation} | ||
\ln\!\left(\alpha\right) = \ln\!\left(0.14\right) - \frac{97.5}{9.5 + \sqrt{\sigma_{v}^{2}}} = \ln\!\left(0.14\right) - \frac{97.5}{9.5 + \sigma_{v}}. | \ln\!\left(\alpha\right) = \ln\!\left(0.14\right) - \frac{97.5}{9.5 + \sqrt{\sigma_{v}^{2}}} = \ln\!\left(0.14\right) - \frac{97.5}{9.5 + \sigma_{v}}. | ||
\label{runningWaveAlpha} | \label{runningWaveAlpha} | ||
\end{equation} | \end{equation} | ||
− | |||
− | [[ | + | \begin{figure}[!htb] |
+ | \center{\input{Data/Beta/rw.tex}} | ||
+ | \caption{Аппроксимация графика $\beta$, бегущая волна} | ||
+ | \label{runningWaveBetaApproximation} | ||
+ | \end{figure} | ||
+ | |||
+ | \begin{figure}[!htb] | ||
+ | \center{\input{Data/Beta/rwLog.tex}} | ||
+ | \caption{Значения функции $\beta \! \left(\sigma_{v}^{2} \right)$ в логарифмических осях, бегущая волна} | ||
+ | \label{runningWaveBetaLog} | ||
+ | \end{figure} | ||
+ | |||
+ | \begin{figure}[!htb] | ||
+ | \center{\input{Data/Alpha/rwLog.tex}} | ||
+ | \caption{Значения функции $\alpha \! \left(\sigma_{v}^{2} \right)$ в логарифмических осях, бегущая волна} | ||
+ | \label{runningWaveAlphaLog} | ||
+ | \end{figure} | ||
+ | |||
+ | \begin{figure}[!htb] | ||
+ | \center{\input{Data/Beta/rw_sw.tex}} | ||
+ | \caption{Сравнение значений $\beta$ для стоячей волны (в масштабе 4:1) и для бегущей волны} | ||
+ | \label{runningAndStandinsWavesBeta} | ||
+ | \end{figure} | ||
− | |||
− | Также одним из результатов поставленных компьютерных экспериментов является обнаружение влияния теплового шума на развал скорости задаваемой механической волны. В случае нулевой дисперсии ( | + | Также одним из результатов поставленных компьютерных экспериментов является обнаружение влияния теплового шума на развал скорости задаваемой механической волны. В случае нулевой дисперсии ($\sigma_v^2 = 0$) наклон графика скорости волны возрастает вплоть до критического значения, после чего образуется дополнительный колебательный процесс (рисунок~\ref{waveEvolution}, слева). Свойства этого процесса были изучены \cite{alexandrov2016solitons}. Однако увеличение $\sigma_v^2$, как было показано ранее, приводит к росту скорости убывания энергии волны и, как следствие, качественно влияет на процесс развала (рисунок~\ref{waveEvolution}, справа). |
− | + | Нагляднее это влияние можно рассмотреть на рисунке~\ref{waveAt27T0}, где приведено развитие скоростей механической волны для разных значений дисперсии по состоянию на один и тот же момент времени. При $\sigma_v^2 = 5$ побочный колебательный процесс выражен неявно, а при $\sigma_v^2 = 20$ вовсе отсутствует. Объясняется это тем, что затухание волны происходит слишком быстро, и наклон графика не достигает критического значения. | |
− | + | При достаточно больших $\sigma_{v}^{2}$ амплитуда волны продолжает падать, но форма остаётся практически неизменной. На рисунке~\ref{waveAt115T0} видно, что график скоростей в кристалле со значением $\sigma_v^2~=~7$ и при $\tau~=~115.5$ по форме близок к графику начального возмущения (см. рисунок~\ref{runningWaveInitialConditions}). | |
− | + | \begin{figure}[!htb] | |
+ | \center{\input{Data/Alpha/rw_swLog.tex}} | ||
+ | \caption{Сравнение значений $\alpha$ для стоячей волны (в масштабе 4:1) и для бегущей волны} | ||
+ | \label{runningAndStandinsWavesAlpha} | ||
+ | \end{figure} | ||
− | [ | + | \begin{figure}[!htb] |
+ | \begin{subfigure}{\linewidth} | ||
+ | \scalebox{0.6}{\input{Data/WaveBreaking/f0t770.tex}}\hfill | ||
+ | \scalebox{0.6}{\input{Data/WaveBreaking/f7t770.tex}} | ||
+ | \end{subfigure} | ||
+ | \par\medskip | ||
+ | \caption{Скорость механической волны на момент $\tau~=~115.5$ для $\sigma_v^2~=~0$ (слева) и $\sigma_v^2~=~7$ (справа)} | ||
+ | \label{waveAt115T0} | ||
+ | \end{figure} | ||
+ | \begin{figure}[!htb] | ||
+ | \begin{subfigure}{\linewidth} | ||
+ | \scalebox{0.6}{\input{Data/WaveBreaking/f0t180.tex}}\hfill | ||
+ | \scalebox{0.6}{\input{Data/WaveBreaking/f1t180.tex}} | ||
+ | \caption{$\sigma_v^2 = 0$ (слева), $\sigma_v^2 = 1$ (справа)} | ||
+ | \end{subfigure} | ||
+ | \par\medskip | ||
+ | \begin{subfigure}{\linewidth} | ||
+ | \scalebox{0.6}{\input{Data/WaveBreaking/f5t180.tex}}\hfill | ||
+ | \scalebox{0.6}{\input{Data/WaveBreaking/f20t180.tex}} | ||
+ | \caption{$\sigma_v^2 = 5$ (слева), $\sigma_v^2 = 20$ (справа)} | ||
+ | \end{subfigure} | ||
+ | \par\medskip | ||
+ | \caption{Скорость механической волны на момент $\tau = 27.0$ для разных значений $\sigma_v^2$} | ||
+ | \label{waveAt27T0} | ||
+ | \end{figure} | ||
− | [ | + | \begin{figure}[!htb] |
+ | \begin{subfigure}{\linewidth} | ||
+ | \scalebox{0.5}{\input{Data/WaveBreaking/f0t1.tex}}\hfill | ||
+ | \scalebox{0.5}{\input{Data/WaveBreaking/f1t1.tex}} | ||
+ | \caption{$\tau = 0$} | ||
+ | \end{subfigure} | ||
+ | \par\medskip | ||
+ | \begin{subfigure}{\linewidth} | ||
+ | \scalebox{0.5}{\input{Data/WaveBreaking/f0t100.tex}}\hfill | ||
+ | \scalebox{0.5}{\input{Data/WaveBreaking/f1t100.tex}} | ||
+ | \caption{$\tau = 15.0$} | ||
+ | \end{subfigure} | ||
+ | \par\medskip | ||
+ | \begin{subfigure}{\linewidth} | ||
+ | \scalebox{0.5}{\input{Data/WaveBreaking/f0t150.tex}}\hfill | ||
+ | \scalebox{0.5}{\input{Data/WaveBreaking/f1t150.tex}} | ||
+ | \caption{$\tau = 22.5$} | ||
+ | \end{subfigure} | ||
+ | \par\medskip | ||
+ | \begin{subfigure}{\linewidth} | ||
+ | \scalebox{0.5}{\input{Data/WaveBreaking/f0t250.tex}}\hfill | ||
+ | \scalebox{0.5}{\input{Data/WaveBreaking/f1t250.tex}} | ||
+ | \caption{$\tau = 37.5$} | ||
+ | \end{subfigure} | ||
+ | \par\medskip | ||
+ | \begin{subfigure}{\linewidth} | ||
+ | \scalebox{0.5}{\input{Data/WaveBreaking/f0t350.tex}}\hfill | ||
+ | \scalebox{0.5}{\input{Data/WaveBreaking/f1t350.tex}} | ||
+ | \caption{$\tau = 52.5$} | ||
+ | \end{subfigure} | ||
+ | \par\medskip | ||
+ | \caption{Эволюция скорости механической волны для $\sigma_v^2~=~0$ (слева) и $\sigma_v^2~=~1$ (справа)} | ||
+ | \label{waveEvolution} | ||
+ | \end{figure} | ||
== Заключение == | == Заключение == | ||
Строка 765: | Строка 790: | ||
Показано, что энергия стоячей волны убывает в 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: | Строка 797: | ||
Результаты работы были получены с использованием вычислительных ресурсов суперкомпьютерного центра Санкт-Петербургского политехнического университета Петра Великого. Каждый расчёт производился на отдельном ядре 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> |