| Текущая версия |
Ваш текст |
| Строка 1: |
Строка 1: |
| | '''''Выпускная квалификационная работа''''' | | '''''Выпускная квалификационная работа''''' |
| − |
| |
| − | '''Направление:''' 01.03.03 – «Механика и математическое моделирование»
| |
| | | | |
| | '''Выполнил:''' студент группы 43604/1 [[Старобинский Егор|Е. Б. Старобинский]] | | '''Выполнил:''' студент группы 43604/1 [[Старобинский Егор|Е. Б. Старобинский]] |
| | | | |
| − | '''Научный руководитель:''' доктор физ.-мат. наук, член-корр. РАН [[Кривцов Антон | А. М. Кривцов]] | + | '''Руководитель:''' доктор физ.-мат. наук, член-корр. РАН [[Кривцов Антон | А. М. Кривцов]] |
| − | | |
| − | Материалы:
| |
| − | * [[Медиа: StarobinskiiThesis.pdf|диплом]]
| |
| − | * [[Медиа: StarobinskiiPoster.pdf|постер]]
| |
| − | * [[Медиа: StarobinskiiPresentation.pdf|презентация]]
| |
| | | | |
| | == Введение == | | == Введение == |
| | 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 name="Fermi1965">Enrico Fermi, J. Pasta, and S. Ulam. Studies of nonlinear problems. Los Alamos Report LA-1940, 978, 1965.</ref> для системы дифференциальных уравнений, описывающую одномерный кристалл с нелинейным взаимодействием между частицами. С помощью синусоидальной волны задавалось начальное возмущение, ожидалось, что с течением времени энергия волны равномерно распределится по всем частицам. Оказалось, что сначала волна действительно теряет свою форму и энергию, но по прошествии достаточного времени почти точно возвращается к исходному состоянию. В дальнейшем это явление было названо «парадоксом Ферми-Паста-Улама-Цингоу», так как исследователям удалось доказать неспособность физики 20-го века описывать взаимодействия на наноуровне. | + | В 1953 году на суперкомпьютере Maniac I группа учёных, состоящая из Энрико Ферми, Станислава Улама, Джона Паста и Мэри Цингоу, численно решила задачу Коши<ref>Enrico Fermi, J. Pasta, and S. Ulam. Studies of nonlinear problems. Los Alamos Report LA-1940, 978, 1965.</ref> для системы дифференциальных уравнений, описывающую одномерный кристалл с нелинейным взаимодействием между частицами. С помощью синусоидальной волны задавалось начальное возмущение, ожидалось, что с течением времени энергия волны равномерно распределится по всем частицам. Оказалось, что сначала волна действительно теряет свою форму и энергию, но по прошествии достаточного времени почти точно возвращается к исходному состоянию. В дальнейшем это явление было названо «парадоксом Ферми-Паста-Улама-Цингоу», так как исследователям удалось доказать неспособность физики 20-го века описывать взаимодействия на наноуровне. |
| | | | |
| | В данной работе рассматривается одна из нелинейных моделей кристалла, описанных Энрико Ферми. Ключевым отличием является введение тепловой энергии через задание дисперсии начальных скоростей. | | В данной работе рассматривается одна из нелинейных моделей кристалла, описанных Энрико Ферми. Ключевым отличием является введение тепловой энергии через задание дисперсии начальных скоростей. |
| Строка 21: |
Строка 14: |
| | == Модель одномерного кристалла == | | == Модель одномерного кристалла == |
| | | | |
| − | === Метод исследования ===
| + | Попл |
| − | | |
| − | Объектом исследования является модель одномерного теплоизолированного кристалла с заданной на нём механической волной. Начальная температура кристалла задаётся случайными скоростями частиц. Дифференциальные уравнения динамики цепочки решаются методом центральных разностей. Производится сравнение решений для бегущей и стоячей волн. Исследуется закон перехода механической энергии волны в тепловую энергию кристалла.
| |
| − | | |
| − | Для численного расчета формулируемых задач создан ряд специальных программ на языках С++, Bash и Mathematica.
| |
| − | | |
| − | === Постановка задачи ===
| |
| − | | |
| − | Рассмотрим модель одномерного кристалла: цепочку одинаковых частиц массы <math>m</math>, соединённых одинаковыми нелинейными пружинами (квадратичная нелинейность).
| |
| − | | |
| − | Каждой частице присвоим свой целочисленный индекс. Принимая <math>u_k</math> за перемещение <math>k</math>-й частицы, а <math>C_1</math> и <math>C_2</math> — за коэффициенты при линейном и квадратичном членах в разложении силы взаимодействия <math>k</math>-й частицы с двумя ближайшими соседями, получим следующее дифференциальное уравнение:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | m\ddot u_k = C_1\left(\left(u_{k + 1} - u_{k}\right) - \left(u_{k} - u_{k - 1}\right)\right) +
| |
| − | C_2\left(\left(u_{k + 1} - u_{k}\right)^2 - \left(u_{k} - u_{k - 1}\right)^2\right),
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | где <math>\dot f</math> — производная <math>f</math> по времени
| |
| − | | |
| − | Перепишем это выражение, раскрыв скобки и введя новые обозначения для коэффициентов при линейной и нелинейной части:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \ddot u_k = \omega_0^2\left(u_{k + 1} - 2u_{k} + u_{k - 1}\right) + \alpha_0^2\left(u_{k + 1}^2 - 2u_{k + 1}u_{k} + 2u_{k}u_{k - 1} - u_{k - 1}^2\right),
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \omega_0 = \sqrt{\frac{C_1}{m}}, \quad \alpha_0 = \sqrt{\frac{C_2}{m}}.
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | Здесь <math>\omega_0</math> соответствует циклической частоте колебаний.
| |
| − | | |
| − | Заметим, что
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \begin{split}
| |
| − | & u_{k + 1}^2 - 2u_{k + 1}u_{k} + 2u_{k}u_{k - 1} - u_{k - 1}^2 = \\
| |
| − | & = u_{k + 1}^2 - 2u_{k + 1}u_{k} + 2u_{k}u_{k - 1} - u_{k - 1}^2 \pm u_{k + 1}u_{k-1}= \\
| |
| − | & = \left(u_{k + 1}^2 - 2u_{k + 1}u_{k} + u_{k - 1}u_{k + 1}\right) - \left(u_{k + 1}u_{k - 1} - 2u_{k}u_{k - 1} + u_{k - 1}^2\right) = \\
| |
| − | & = u_{k + 1}\left(u_{k + 1} - 2u_{k} + u_{k - 1}\right) - u_{k - 1}\left(u_{k + 1} - 2u_{k} + u_{k - 1}\right) = \\
| |
| − | & = \left(u_{k + 1} - 2u_{k} + u_{k - 1}\right)\left( u_{k + 1} - u_{k - 1} \right).
| |
| − | \end{split}
| |
| − | \label{replace}
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | Воспользуемся этой заменой, и тогда уравнения динамики цепочки примут следующий вид:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \ddot u_k = \left(u_{k + 1} - 2u_{k} + u_{k - 1} \right) \left( \omega_0^2 + \alpha_0^2 \left( u_{k + 1} - u_{k - 1} \right) \right).
| |
| − | \label{chainEquations}
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | Аналогичное уравнение использовалось в оригинальной работе Ферми, Паста и Улама<ref name="Fermi1965"/>.
| |
| − | | |
| − | Механическую энергию системы в начальный момент времени (то есть, при <math>t = 0</math>) будем задавать с помощью синусоидальной волны. Из-за теплового движения частиц волна будет терять свою форму, а её энергия перейдёт в тепловую энергию системы. Для описания этого процесса численно решим уравнения динамики кристалла.
| |
| − | | |
| − | Чтобы замкнуть систему из <math>k</math> уравнений динамики, определим граничные условия, а также по <math>k</math> начальных условий на перемещения и скорости. Для задания начальных скоростей воспользуемся генератором случайных чисел с равномерным распределением. При этом дисперсию скоростей будем задавать в соответствии с требуемым температурным профилем.
| |
| − | | |
| − | === Начальные и граничные условия ===
| |
| − | | |
| − | Чтобы соответствовать идеальной модели, исследуемая цепочка должна иметь бесконечную длину. Для проведения численного моделирования будем рассматривать кристалл конечной длины с заданным условием периодичности (граничные условия Борна-Кармана):
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \label{boundaryConditions}
| |
| − | u\!\left(x + L\right) = u\!\left(x\right), \quad L \gg a,
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | где <math>L</math> — длина кристалла, <math>a</math> — равновесное расстояние между частицами.
| |
| − | | |
| − | Таким образом, условие теплоизоляции выполняется.
| |
| − | | |
| − | В рамках поставленной задачи будем рассматривать два вида движения в кристалле: движение механической волны и тепловое движение (тепловой шум). Соответственно, начальные условия для скоростей частиц также будут складываться из двух компонент. Меняя начальные условия, рассмотрим решение задачи для следующих постановок:
| |
| − | * на кристалле задана бегущая волна;
| |
| − | * на кристалле задана стоячая волна.
| |
| − | | |
| − | ==== Начальные условия для бегущей волны ====
| |
| − | | |
| − | Рассмотрим уравнение бегущей волны (пусть волна распространяется в одном направлении):
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \label{runningWaveEquation1}
| |
| − | \dot u\!\left(x\right) = A \sin\!\left(k_0 x - \omega_0 t + \phi \right),
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | где <math>A</math> — амплитуда колебаний, <math>k_0</math> — волновое число, <math>\phi</math> — фаза колебаний.
| |
| − | | |
| − | Величину <math>u\!\left(x\right)</math> на момент времени <math>t</math> получим интегрированием уравнения волны:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \label{runningWaveEquation2}
| |
| − | u\!\left(x\right) = \int\limits_0^t \! \dot u\!\left(x\right) \,\mathrm{d}t = \frac{1}{\omega_0} A \cos\!\left(k_0 x - \omega_0 t + \phi \right).
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | Пусть длина волны совпадает с длиной кристалла, тогда
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | k_0~\mathrel{\stackrel{\rm def}=} \frac{2 \pi}{\lambda} = \frac{2 \pi}{L},
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | где <math>\lambda</math> — длина волны.
| |
| − | | |
| − | Фазу колебаний <math>\phi</math> положим равной нулю и запишем уравнения для <math>u</math> и <math>\dot u</math> для начального момента времени (<math>t=0</math>):
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \label{initialMechanicalVelocity1}
| |
| − | \dot u\!\left(x\right)|_{t=0} = A \sin\!\left(\frac{2 \pi x}{L} \right).
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \label{initialMechanicalDisplacement1}
| |
| − | u\!\left(x\right)|_{t=0} = \frac{1}{\omega_0} A \cos\!\left(\frac{2 \pi x}{L}\right).
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | На рисунке 1 приведён график начальных перемещений и скоростей частиц кристалла при заданной бегущей волне.
| |
| − | | |
| − | [[File:StarobinskiiThesis1.png|framed|center|Рисунок 1. Начальные условия для бегущей волны]]
| |
| − | | |
| − | ==== Начальные условия для стоячей волны ====
| |
| − | | |
| − | Начальные условия, описанные в предыдущем пункте, позволяют задать бегущую волну. Получим подобные условия также для стоячей волны. Для этого рассмотрим две волны амплитуды <math>\frac{A}{2}</math>, бегущие в противоположных направлениях с разностью фаз <math>\pi</math>:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | v_1\!\left(x\right) = \frac{A}{2} \sin\!\left(k_0 x - \omega_0 t \right), \quad v_2\!\left(x\right) = \frac{A}{2} \sin\!\left(- k_0 x - \omega_0 t + \pi \right),
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | Результатом наложения этих двух волн будет стоячая волна, описываемая следующим уравнением:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \label{standingWaveEquation1}
| |
| − | \dot u\!\left(x\right) = v_1\!\left(x\right) + v_1\!\left(x\right) = A \cos\!\left(\omega_0 t \right)\sin\left(k_0 x \right).
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | Аналогично, проинтегрируем это уравнение по времени <math>t</math>:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | u\!\left(x\right) = \int\limits_0^t \! \dot u\!\left(x\right) \,\mathrm{d}t = \frac{A}{\omega_0} \sin\!\left(\omega_0 t \right)\sin\left(k_0 x \right).
| |
| − | \label{standingWaveEquation2}
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | Подставляя в уравнения для <math>u</math> и <math>\dot u</math> <math>k_0 = \frac{2 \pi}{L}</math> и <math>t = 0</math>, получим начальные условия для стоячей волны:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \label{initialMechanicalVelocity2}
| |
| − | \dot u\!\left(x\right)|_{t=0} = A \sin\!\left(\frac{2 \pi x}{L} \right).
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \label{initialMechanicalDisplacement2}
| |
| − | u\!\left(x\right)|_{t=0} = 0.
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | Соответствующий график начальных перемещений и скоростей частиц цепочки приведён на рисунке 2.
| |
| − | | |
| − | [[File:StarobinskiiThesis2.png|framed|center|Рисунок 2. Начальные условия для стоячей волны]]
| |
| − | | |
| − | Функции скоростей в начальный момент времени для бегущей и стоячей волн совпадают, различаются только начальные перемещения.
| |
| − | | |
| − | ==== Начальные условия теплового движения ====
| |
| − | | |
| − | Введём оператор усреднения функции по её аргументу:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \left<f\right> = \frac{1}{b - a}\int\limits_{a}^{b} \! f \! \left( x \right) \, \mathrm{d}x,
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | где <math>\left[a; b\right]</math> — отрезок, на котором проводится усреднение.
| |
| − | | |
| − | Оператор нахождения среднего по координате для цепочки длины <math>L</math> в таком случае записывается следующим образом:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \left<f\right> = \frac{1}{L}\int\limits_{0}^{L} \! f \! \left( x \right) \, \mathrm{d}x,
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | Чтобы получить выражение для тепловой компоненты скорости, запишем определение кинетической температуры в одномерной постановке:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \frac{1}{2}kT\!\left(x\right) = \frac{m \left< v^{2} \right>}{2},
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | где <math>k</math> — постоянная Больцмана, <math>T(x)</math> — профиль температуры, <math>v</math> — скорость.
| |
| − | | |
| − | Из этого равенства выразим <math>\left< v^{2} \right></math>:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \left< v^{2} \right> = \frac{kT\!\left(x\right)}{m}.
| |
| − | \label{vTemp2Average}
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | Пусть в начальный момент времени температура в кристалле распределена равномерно. Соответствующее распределение для скоростей будем задавать с помощью случайного числа <math>\rho</math> в диапазоне <math>\left[-1;1\right]</math>. Определим среднее значение <math>\rho^{2}</math>:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \left< \rho^{2} \right> = \frac{1}{2}\int\limits_{-1}^{1} \! \rho^{2} \, \mathrm{d}\rho = \frac{1}{2} \frac{\rho^{3}}{3} \bigg|_{-1}^{1} = \frac{1}{3}.
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | Таким образом, <math>3 \left< \rho^{2} \right> = 1</math>. Воспользуемся этим соотношением для построения требуемого профиля <math>T \! \left( x \right)</math>:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | T \! \left( x \right) = 3 \left< \rho^{2} \right> T_0,
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | где <math>T_0</math> — значение температуры при <math>t = 0</math>.
| |
| − | | |
| − | Перепишем выражение для <math>\left< v^{2} \right></math>:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \left< v^{2} \right> = \frac{ 3 k T_0}{m} \left< \rho^{2} \right>.
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | Выразим <math>v^{2}\!\left(x\right)</math> и <math>v\!\left(x\right)</math>:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | v^{2}\!\left(x\right) = \frac{ 3 k T_0}{m} \rho^{2}\!\left(x\right),
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | v\!\left(x\right) = \sqrt{\frac{3 k T_0}{m}} \rho\!\left(x\right).
| |
| − | \label{vTemp}
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | Введём величину дисперсии скоростей <math>\sigma_{v}^{2}</math> так, чтобы при <math>\sigma_{v}^{2} = 1</math> в начальный момент времени механическая энергия волны равнялась тепловой энергии цепочки:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | T_0 = \frac{m}{k} A^{2} \sigma_{v}^{2}
| |
| − | \label{initialTemp}
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | Тогда значение <math>\sigma_{v}^{2} = 10</math> можно трактовать следующим образом: при <math>t~=~0</math> энергия теплового шума в кристалле в <math>10</math> раз превосходит механическую энергию волны.
| |
| − | | |
| − | С учётом этого соотношения получил начальный профиль скоростей:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | v\!\left(x\right)|_{t = 0} = A \sqrt{3 \sigma_{v}^{2}} \rho \! \left( x \right).
| |
| − | \label{initialThermalVelocity}
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | Таким образом, были получены выражения для обеих компонент скорости: скорость механической волны задаётся синусоидальной функцией, тепловой шум определяется стохастической функцией. Итоговое условие на начальные скорости в кристалле принимает следующий вид:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \left. \dot u(x) \right |_{t = 0} = A \left( \sin\!\left( \frac{2 \pi x}{L} \right) + \sqrt{3 \sigma_{v}^{2}} \rho \! \left( x \right) \right).
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | Начальные перемещения задаются в соответствии с типом рассматриваемой волны.
| |
| − | | |
| − | === Вычисление механической энергии ===
| |
| − | | |
| − | В начальный момент скорость имеет распределение, близкое к синусоидальному. С течением времени цепочка теряет эту форму (в силу нелинейного взаимодействия между частицами), а механическая энергия переходит в тепловую.
| |
| − | | |
| − | Механическая энергия системы <math>E\!\left(x\right)</math> складывается из кинетической энергии <math>K\!\left(x\right)</math> и потенциальной энергии <math>\Pi\!\left(x\right)</math> следующий образом:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | K\!\left( x \right) = \frac{m}{2}\dot u^2\!\left(x\right),
| |
| − | \quad \Pi \!\left( x \right) = \frac{C_1}{2}\varepsilon^2\!\left(x\right),
| |
| − | \quad E\!\left(x\right) = K\!\left(x\right) + \Pi\!\left(x\right),
| |
| − | \label{energyEquations}
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | где <math>\varepsilon\!\left(x\right)</math> — возникающие деформации.
| |
| − | | |
| − | Так как нас интересует энергия не всей системы (включающей тепловой шум), а только заданной механической волны, будем искать энергию первой синусоидальной моды. Перепишем выражения для энергий:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | K\!\left( x \right) = \frac{m}{2}\left(\int\limits_L \! \dot u\!\left(x\right) \sin\!\left( \frac{2 \pi x}{L} \right) \,\mathrm{d}x\right)^2
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \Pi \!\left( x \right) = \frac{C_1}{2}\left(\int\limits_L \! \varepsilon\!\left(x\right) \cos\!\left( \frac{2 \pi x}{L} \right) \,\mathrm{d}x\right)^2
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | E\!\left(x\right) = \frac{m}{2}\left(\int\limits_L \! \dot u\!\left(x\right) \sin\!\left( \frac{2 \pi x}{L} \right) \,\mathrm{d}x\right)^2
| |
| − | + \frac{C_1}{2}\left(\int\limits_L \! \varepsilon\!\left(x\right) \cos\!\left( \frac{2 \pi x}{L} \right) \,\mathrm{d}x\right)^2
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | Первая синусоидальная мода также возникает и в тепловом шуме. Математическое ожидание от её энергии вычислим, усреднив по ансамблю реализаций систему без заданной волны.
| |
| − | | |
| − | Определим значения <math>K\!\left(x\right)</math> и <math>\Pi\!\left(x\right)</math> при <math>t = 0</math>:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | K\!\left( x \right)|_{t=0} = \frac{m}{2}\left(A \int\limits_L \! \sin^2\!\left( \frac{2 \pi x}{L} \right) \,\mathrm{d}x\right)^2,
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \label{initialPotentialEnergy}
| |
| − | \Pi\!\left( x \right)|_{t=0} = \frac{C_1}{2}\left( \int\limits_L \! \varepsilon\!\left(x\right)|_{t = 0} \cos\!\left( \frac{2 \pi x}{L} \right) \,\mathrm{d}x\right)^2
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | Интеграл в данной формуле равен нулю в случае стоячей волны (в начальный момент времени нет деформаций). Покажем, что этот интеграл также будет равен нулю в случае бегущей волны:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \begin{split}
| |
| − | & \int\limits_L \! \varepsilon\!\left(x\right)|_{t = 0} \cos\!\left( \frac{2 \pi x}{L} \right) \,\mathrm{d}x =
| |
| − | A \int\limits_L \! \sin\!\left(\frac{2 \pi x}{L}\right) \cos\!\left( \frac{2 \pi x}{L} \right) \,\mathrm{d}x = \\
| |
| − | & = \frac{A}{2} \int\limits_L \! \sin\!\left(\frac{4 \pi x}{L}\right) \,\mathrm{d}x = 0.
| |
| − | \end{split}
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | То есть, в начальный момент времени потенциальная энергия равна нулю, а механическая энергия задаваемой волны равна
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \label{initialMechanicalEnergy}
| |
| − | E\!\left( x \right)|_{t=0} = \frac{m}{2}\left(A \int\limits_L \! \sin^2\!\left( \frac{2 \pi x}{L} \right) \,\mathrm{d}x\right)^2 = \frac{m A^{2} L^{2}}{8}
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | Начальное значение механической энергии выберем в качестве масштаба энергии.
| |
| − | | |
| − | Чтобы определить закон, по которому происходит переход механической энергии в тепловую, введём безразмерный параметр <math>E^*\!\left( x \right) = \frac{E\!\left( x \right)}{E\!\left( x \right)|_{t=0}}</math> — нормированное значение механической энергии цепочки:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | E^* = \left( \frac{\int\limits_L \! \dot u\!\left(x\right) \sin\!\left( \frac{2 \pi x}{L} \right) \,\mathrm{d}x}{A \int\limits_L \! \sin^2 \left( \frac{2 \pi x}{L} \right) \,\mathrm{d}x} \right)^2 +
| |
| − | \omega_0^2 \left( \frac{\int\limits_L \! \varepsilon\!\left(x\right) \cos\!\left( \frac{2 \pi x}{L} \right) \,\mathrm{d}x}{A \int\limits_L \! \sin^2 \left( \frac{2 \pi x}{L} \right) \,\mathrm{d}x} \right)^2.
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | Решение задачи свелось к нахождению <math>E^*\!\left( x \right)</math>.
| |
| − | | |
| − | === Численное решение дифференциального уравнения ===
| |
| − | | |
| − | Для решения системы уравнений движения частиц воспользуемся методом центральных разностей. Перепишем граничные условия:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | u_{k + N} = u_{k}, \quad \dot u_{k + N} = \dot u_{k},
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | где <math>N \gg 1</math> — число независимых частиц.
| |
| − | | |
| − | Вычисление правой части уравнения для <math>\ddot u</math> аналогично решению с помощью разностной схемы. Рассчитанные таким образом ускорения частиц пересчитываем в скорости:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \dot u_k\!\left(t + \Delta t\right) = \dot u_k\!\left(t\right) + \ddot u_k\!\left(t\right)\Delta t.
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | Затем скорости пересчитываем в смещения частиц. Важно отметить, что по данному методу расчёт скоростей производится на промежуточных шагах по времени относительно расчёта смещений.
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | u_k\!\left(t + \Delta t\right) = u_k\!\left(t\right) + \dot u_k\!\left(t\right)\Delta t.
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | Значение <math>\varepsilon\!\left(x\right)</math> вычисляем следующим образом:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \varepsilon = \sum_{i = 1}^{N - 1} \left(u_{i} - u_{i - 1}\right).
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | === Параметры расчёта ===
| |
| − | | |
| − | Для достоверного численного решения дифференциального уравнения относительно <math>\left(u\!\left(x\right), \dot u\!\left(x\right)\right)</math> длина цепочки должна стремиться к бесконечной, чтобы соответствовать модели идеального кристалла. При проведении численных экспериментов удовлетворительная точность расчётов достигалась на цепочках из 10 миллионов частиц и более.
| |
| − | | |
| − | Статистические характеристики системы также можно вычислить путём усреднения бесконечного множества реализаций (каждая реализация — задача Коши для цепочки). Благодаря этому подходу, предложенному в работа<ref>A. M. Krivtsov. Energy distribution in one-dimensional crystal. Dokl. Phys. 60 (9), page 407, 2015.</ref><ref>A. M. Krivtsov. On unsteady heat conduction in a harmonic crystal. ArXiv e-prints, September 2015.</ref>, достигается эффективность проводимых численных экспериментов.
| |
| − | | |
| − | Результаты, приводимые далее, получены усреднением <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>m</math> на пружине жёсткости <math>C_1</math>:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \quad T_0~\mathrel{\stackrel{\rm def}=} \frac{2 \pi}{\omega_0} = 2 \pi \sqrt{\frac{m}{C_1}}.
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | '''Таблица 1. Параметры при моделировании стоячей волны'''
| |
| − | {| class="wikitable"
| |
| − | |-
| |
| − | ! <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.02</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.02</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.02</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.02</math>
| |
| − | |}
| |
| − | | |
| − | '''Таблица 2. Параметры при моделировании бегущей волны'''
| |
| − | {| class="wikitable"
| |
| − | |-
| |
| − | ! <math>\sigma_v^2</math> !! <math>\Delta \sigma_v^2</math> !! <math>\Delta t</math> !! Постоянные параметры
| |
| − | |-
| |
| − | | <math>[0; 10]</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; 195]</math> || <math>5</math> || <math>0.05~T_0</math> || <math>\omega_0 = 1, \quad \alpha_0 = 1, \quad A = 0.01</math>
| |
| − | |-
| |
| − | | <math>[200; 300]</math> || <math>5</math> || <math>0.005~T_0</math> || <math>\omega_0 = 1, \quad \alpha_0 = 1, \quad A = 0.01</math>
| |
| − | |}
| |
| − | | |
| − | Значения параметров проводимых экспериментов приведены в таблицах 1 и 2. Моделируемое время ограничено: <math>0 \leq t \leq 7 \cdot 10^4~T_0</math>.
| |
| − | | |
| − | == Результаты ==
| |
| − | | |
| − | === Преобразование механической энергии стоячей волны ===
| |
| − | | |
| − | На рисунке 3 приведены вычисленные значения нормированной механической энергии стоячей волны <math>E^*\!\left(\tau\right)</math> для нескольких значений <math>\sigma_v^2</math> (здесь и далее <math>\tau = \frac{t}{T_0}</math> – безразмерное время).
| |
| − | | |
| − | [[File:StarobinskiiThesis3.png|framed|center|Рисунок 3. Изменение нормированной механической энергии стоячей волны в зависимости от <math>t</math>]]
| |
| − | | |
| − | Для одномерного кристалла при отсутствии теплового шума (<math>\sigma_v^2 = 0</math>) была установлена<ref name="presTsvetkovAPM16">D. V. Tsvetkov and A. M. Krivtsov. Energy distribution in one-dimensional crystal. Book of papers 24th International Congress of Theoretical and Applied Mechanics, pages 2450–2451, 2016.</ref> следующая зависимость для механической энергии волны:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | E^* \sim e^{-\tau^2}.
| |
| − | \label{expectationTau2}
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | Подобная зависимость не выполняется при ненулевой дисперсии скоростей частиц цепочки: механическая энергия убывает быстрее при бóльших значениях <math>\sigma_v^2</math>. Предположим, что <math>E^*</math> уменьшается по модифицированному экспоненциальному закону и показатель степени при <math>\tau</math> в общем случае неизвестен:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | E^* = e^{-\alpha \tau^\beta},
| |
| − | \label{expectationTauBeta}
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | где <math>\beta</math> — функция от <math>\sigma_v^2</math>, <math>\alpha</math> — функция от <math>\sigma_v^2</math> и <math>\lambda</math>.
| |
| − | | |
| − | В рамках гипотезы <math>\beta</math> имеет определённое значение для каждого значения <math>\sigma_v^2</math>. При <math>\sigma_v^2 = 0</math> показатель <math>\beta = 2</math><ref name="presTsvetkovAPM16"/>, по мере же увеличения <math>\sigma_v^2</math>, согласно рисунку 3, величина <math>\beta</math> должна уменьшаться.
| |
| − | | |
| − | Если предположение верно, то показатель <math>\beta</math> в зависимости от величины теплового шума можно определить, перестроив график энергии в логарифмических осях.
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \begin{split}
| |
| − | E^* & = e^{-\alpha \tau^\beta} \Rightarrow \\
| |
| − | \Rightarrow \ln\!\left(E^*\right) & = -\alpha \tau^\beta \Rightarrow \\
| |
| − | \Rightarrow \ln\!\left(-\ln\!\left(E^*\right)\right) & = \ln\!\left(\alpha \tau^\beta\right) \Rightarrow \\
| |
| − | \Rightarrow \ln\!\left(\ln\!\left(\frac{1}{E^*}\right)\right) & = \ln\!\left(\alpha\right) + \ln\!\left(\tau^\beta\right) \Rightarrow \\
| |
| − | \Rightarrow \ln\!\left(\ln\!\left(\frac{1}{E^*}\right)\right) & = \ln\!\left(\alpha\right) + \beta\ln\!\left(\tau\right).
| |
| − | \end{split}
| |
| − | \label{expectationTauBetaLine}
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | Обозначим <math>\ln\!\left(\tau\right)</math> за <math>x</math>, а <math>\ln\!\left(\ln\!\left(\frac{1}{E^*}\right)\right)</math> за <math>y</math>, тогда уравнение сводится к линейному. Угол наклона прямой в осях <math>\left(\ln\!\left(\tau\right); \ln\!\left(\ln\!\left(\frac{1}{E^*}\right)\right) \right)</math>, выраженный в радианах, будет численно равен <math>\beta</math> для определённого значения дисперсии.
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | y\!\left(x\right) = \ln\!\left(\alpha\right) + \beta x
| |
| − | \label{lineEquation}
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | [[File:StarobinskiiThesis4.png|framed|center|Рисунок 4. Изменение механической энергии стоячей волны с течением времени в логарифмических осях]]
| |
| − | | |
| − | Соответствующее преобразование графика энергии приведено на рисунке 4. Получившиеся функции успешно аппроксимируются наклонными прямыми, что подтверждает высказанную ранее гипотезу.
| |
| − | | |
| − | [[File:StarobinskiiThesis5.png|framed|center|Рисунок 5. Значения <math>\beta</math> при разных значениях дисперсии, стоячая волна]]
| |
| − | | |
| − | Найденные значения <math>\beta</math> приведены на рисунке 5. Тепловая энергия в реальном материале может превосходить механическую на несколько порядков. Чтобы оценить значение <math>\beta</math> в таком случае, построим аппроксимацию этого графика и определим предельное значение показателя при условии <math>\sigma_v^2 \rightarrow \infty</math>.
| |
| − | | |
| − | [[File:StarobinskiiThesis6.png|framed|center|Рисунок 6. Значения функции <math>\beta \! \left(\sigma_{v}^{2} \right)</math> в логарифмических осях, стоячая волна]]
| |
| − | | |
| − | Перестроим график из осей <math>\left(\sigma_{v}^{2} ; \beta \right)</math> в <math>\left(\ln\!\left(\sigma_{v}^{2}\right) ; \ln\!\left(\beta - 1\right) \right)</math>. Результат приведён на рисунке 6. Начиная с <math>x = 1.5</math> значения могут быть аппроксимированы прямой:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | y\!\left(x\right) = \ln\!\left(1.08\right) - 0.25 x,
| |
| − | \label{standingWaveApproximation}
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | где <math>x = \ln\!\left(\sigma_{v}^{2}\right)</math>, <math>y = \ln\!\left(\beta\right)</math>.
| |
| − | | |
| − | Выражая функцию <math>\beta</math>, получаем следующее:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \beta \! \left(\sigma_{v}^{2} \right) = 1 + \frac{1.08}{\sqrt[4]{\sigma_{v}^{2}}} = 1 + \frac{1.08}{\sqrt{\sigma_{v}}},
| |
| − | \label{standingWaveApproximation2}
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | где <math>\sigma_{v}</math> — среднеквадратическое отклонение скоростей частиц, отнесённое к своему начальному значению (при <math>\tau~=~0</math>).
| |
| − | | |
| − | [[File:StarobinskiiThesis7.png|framed|center|Рисунок 7. Аппроксимация графика <math>\beta</math>, стоячая волна]]
| |
| − | | |
| − | [[File:StarobinskiiThesis8.png|framed|center|Рисунок 8. Аппроксимация графика <math>\beta</math>, стоячая волна]]
| |
| − | | |
| − | Аппроксимированный график <math>\beta</math> от <math>\sigma_{v}^{2}</math> приведён на рисунках 7 и 8. Найдём предельное значение <math>\beta</math>:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \lim_{\sigma_{v}^{2} \to \infty} \beta \! \left(\sigma_{v}^{2} \right) =
| |
| − | \lim_{\sigma_{v}^{2} \to \infty} 1 + \frac{1.08}{\sqrt[4]{\sigma_{v}^{2}}} = 1 + \frac{1.08}{\infty} = 1.
| |
| − | \label{standingWaveLimit}
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | Можно считать, для больших значений <math>\sigma_{v}^{2}</math> механическая энергия будет переходить в тепловую по закону <math>E^* = e^{-\alpha \tau}</math>. Сходимость <math>\beta</math> к единице также видна на рисунке 8. Уже при <math>\sigma_{v} = 12 \cdot 10^3</math> ошибка в определении <math>\beta</math> единицей составит менее процента.
| |
| − | | |
| − | [[File:StarobinskiiThesis9.png|framed|center|Рисунок 9. Аппроксимация графика <math>\alpha</math>, стоячая волна]]
| |
| − | | |
| − | Исследуем поведение <math>\alpha</math> в зависимости от величины теплового шума. Для этого воспользуемся уравнениями прямых, полученных аппроксимацией механической энергии <math>E^*</math> в логарифмических осях (см. рисунок 6). Подставив в уравнения прямых <math>x = 0</math>, определим значения <math>\ln\!\left(\alpha\right)</math> как функции от <math>\sigma_{v}^{2}</math>. Результат приведён на рисунке 9. Была получена следующая формула для аппроксимации <math>\alpha</math>:
| |
| − | | |
| − | <math>
| |
| − | \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}}.
| |
| − | \label{standingWaveAlpha}
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | Аналогично выражению для <math>\beta</math>, это уравнение имеет предельное значение при <math>\sigma_{v}^{2}~\to~\infty</math>:
| |
| − | | |
| − | <math>
| |
| − | \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).
| |
| − | \label{standingWaveAlphaLimit}
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | Отсюда следует, что
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \lim_{\sigma_{v}^{2} \to \infty} \alpha = 0.14.
| |
| − | \label{standingWaveAlphaLimit2}
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | При <math>\sigma_{v}^{2} = 12 \cdot 10^3</math> значение <math>\alpha</math> отличается от предельного примерно на <math>0.2\%</math>.
| |
| − | | |
| − | Дополнительное рассмотрение цепочек с кратным числом частиц привело к следующему результату: уменьшение длины кристалла в <math>2</math> раза приводит к увеличению предельного значения <math>\alpha</math> в <math>4</math> раза. То есть, <math>\alpha \sim k^{2}</math> (<math>k</math> — волновое число, <math>k \sim \frac{1}{\lambda}</math>). Таким образом, в реальных материалах изменение нормированной механической энергии <math>E^*</math> может быть описано экспоненциальным законом:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | E^*\!\left(\tau\right) = e^{-0.14 k_{*}^{2} \tau},
| |
| − | \label{energySolution}
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | где <math>k_{*} = \frac{k}{k_{1000}}</math>, <math>k_{1000}</math> — значение <math>k</math> при <math>N = 1000</math>.
| |
| − | | |
| − | Сравним полученную зависимость с решением уравнения продольных колебаний стержня с внутренним трением. В используемых обозначениях уравнение колебаний может быть записано таким образом:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \frac{\partial^2 u}{\partial \tau^2} = \frac{E T_{0}^{2}}{\rho a^2} \frac{\partial^2 u}{\partial x^2}
| |
| − | + 2 \gamma \frac{\partial^3 u}{\partial \tau \partial x^2},
| |
| − | \label{oscillationEquation}
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | где <math>E</math> — модуль Юнга, <math>\rho</math> — плотность, <math>2 \gamma</math> — безразмерный коэффициент при диссипативном члене, <math>x</math> — безразмерная переменная расстояния.
| |
| − | | |
| − | Решение дифференциального уравнения имеет вид:
| |
| − | | |
| − | <math>
| |
| − | \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),
| |
| − | \label{oscillationEquationSolution}
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | где <math>\omega_{c} = k_{c} \sqrt{\frac{E T_{0}^{2}}{\rho a^2} - \gamma}</math>, <math>k_{c} = \frac{a}{2} k</math>, <math>A_{c} = \frac{a}{T0} A</math>.
| |
| − | | |
| − | Данный результат был представлен в работе К. Ю. Аристовича<ref>К. Ю. Аристович. Зависимость макроскопических параметров колебаний кристаллического стержня от микроструктуры. Труды XXI международной конференции «Математическое моделирование в механике сплошных сред. Методы граничных и конечных элементов», 2:23–35, 2006.</ref> и хорошо согласуется с формулой для стоячей волны (см. раздел 2.3.2), переписанной в безразмерных величинах:
| |
| − | | |
| − | <math>
| |
| − | \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),
| |
| − | \label{standingWaveEquation3}
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | где <math>\hat{\omega_{c}} = T_{0} \omega_{0}.</math>
| |
| − | | |
| − | Сопоставим эти два решения и определим значения <math>\gamma</math>:
| |
| − | | |
| − | <math>
| |
| − | \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.
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | Покажем необратимость преобразования механической энергии в тепловую. Парадокс Ферми-Паста-Улама продемонстрировал явление возврата формы и механической энергии волны<ref name="Fermi1965"/>. Обратный переход энергии цепочки в механическую энергию был описан для одномерного нелинейного кристалла без теплового движения<ref>Д. В. Цветков. Распределение тепла в одномерном кристалле. Диссертация магистра, pages 19–23, 2015.</ref>. Такой переход также был экспериментально получен в данной работе. Рассмотрим влияние теплового шума на этот процесс.
| |
| − | | |
| − | [[File:StarobinskiiThesis10.png|framed|center|Рисунок 10. Изменение энергии стоячей волны]]
| |
| − | | |
| − | На рисунке 10 можно заметить дополнительный изгиб графика энергии <math>E^*\!\left(T\right)</math> при малых значениях <math>\sigma_v^2</math>. Так, для <math>\sigma_v^2 = 0</math> возврат механической энергии наблюдается примерно в диапазоне от <math>\tau=60</math> до <math>\tau=75</math>. Для <math>\sigma_v^2 = 1</math> этот диапазон уже: от <math>\tau=65</math> до <math>\tau=70</math>, при этом прирост энергии значительно меньше. Для высоких значений <math>\sigma_v^2</math> обратный переход энергии не удаётся продемонстрировать. Можно утверждать, что процесс преобразования механической энергии в тепловую в больших системах с высокими значениями отношения тепловой энергии к механической (соответствует реальным системам) является необратимым.
| |
| − | | |
| − | === Преобразование механической энергии бегущей волны ===
| |
| − | | |
| − | Рассмотрим график убывания механической энергии бегущей волны с течением времени (рисунок 13). Качественно он совпадает с графиком энергии стоячей волны (рисунок 3). При отсутствии теплового шума для энергии волны также выполняется соотношение <math>E^* \sim e^{-\tau^2}</math>.
| |
| − | | |
| − | [[File:StarobinskiiThesis11.png|framed|center|Рисунок 11. Изменение энергии бегущей волны в зависимости от <math>\tau</math>]]
| |
| − | | |
| − | Чтобы найти зависимость <math>\beta\!\left(\sigma_v^2\right)</math>, аналогичным образом перестроим график энергии <math>E^{*}\!\left(\tau\right)</math> в логарифмических осях (см. рисунок 12) и определим угол наклона <math>\beta</math>.
| |
| − | | |
| − | [[File:StarobinskiiThesis11.5.png|framed|center|Рисунок 12. Изменение механической энергии стоячей волны с течением времени в логарифмических осях]]
| |
| − | | |
| − | | |
| − | Вид <math>\beta</math> приведён на рисунке 13. Аппроксимация также была получена переводом осей <math>\left(\sigma_{v}^{2} ; \beta \right)</math> в <math>\left(\ln\!\left(\sigma_{v}^{2}\right) ; \ln\!\left(\beta\right)\right)</math> (см. рисунок 14).
| |
| − | | |
| − | [[File:StarobinskiiThesis12.png|framed|center|Рисунок 13. Аппроксимация графика <math>\beta</math>, бегущая волна]]
| |
| − | | |
| − | [[File:StarobinskiiThesis13.png|framed|center|Рисунок 14. Значения функции <math>\beta \! \left(\sigma_{v}^{2} \right)</math> в логарифмических осях, бегущая волна]]
| |
| − | | |
| − | Для <math>\beta</math> был получен следующий вид уравнения аппроксимации:
| |
| − | | |
| − | <math>
| |
| − | \begin{equation}
| |
| − | \beta \! \left(\sigma_{v}^{2} \right) = 1 + \frac{1.55}{\sqrt[4]{\sigma_{v}^{2}}} = 1 + \frac{1.55}{\sqrt{\sigma_{v}}}.
| |
| − | \label{runningWaveApproximation}
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | При условии, что <math>\sigma_{v}^{2} \to \infty</math>, эта функция стремится к единице, аналогично аппроксимирующей функции в случае стоячей волны.
| |
| − | | |
| − | Значения <math>\alpha</math> для различных <math>\sigma_{v}^{2}</math> также были получены (рисунок 15). Предел аппроксимирующей функции равен соответствующему значению для аппроксимирующей функции стоячей волны. Таким образом, в реальных материалах экспоненциальная формула для механической энергии выполняется вне зависимости от типа заданной механической волны.
| |
| − | | |
| − | [[File:StarobinskiiThesis14.png|framed|center|Рисунок 15. Значения функции <math>\alpha \! \left(\sigma_{v}^{2} \right)</math> в логарифмических осях, бегущая волна ]]
| |
| − | | |
| − | На рисунке 16 приведён график значений <math>\beta</math> для стоячей волны, растянутый по горизонтальной оси в 4 раза. В таком масштабе экспериментальные значения <math>\beta</math> для бегущей и стоячей волн могут быть аппроксимированы одной и той же кривой (функцией для бегущей волны). То есть, в случае бегущей волны было получено в 4 раза более медленное убывание механической энергии, чем в случае стоячей. Такая же пропорциональность наблюдается для значений <math>\alpha</math> (см. рисунок 17).
| |
| − | | |
| − | <math>
| |
| − | \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}}.
| |
| − | \label{runningWaveAlpha}
| |
| − | \end{equation}
| |
| − | </math>
| |
| − | | |
| − | [[File:StarobinskiiThesis15.png|framed|center|Рисунок 16. Сравнение значений <math>\beta</math> для стоячей волны (в масштабе 4:1) и для бегущей волны]]
| |
| − | | |
| − | [[File:StarobinskiiThesis16.png|framed|center|Рисунок 17. Сравнение значений <math>\alpha</math> для стоячей волны (в масштабе 4:1) и для бегущей волны]]
| |
| − | | |
| − | Также одним из результатов поставленных компьютерных экспериментов является обнаружение влияния теплового шума на развал скорости задаваемой механической волны. В случае нулевой дисперсии (<math>\sigma_v^2 = 0</math>) наклон графика скорости волны возрастает вплоть до критического значения, после чего образуется дополнительный колебательный процесс (рисунок 19, слева). Свойства этого процесса были изучен<ref>С. Д. Александров. Исследование свойств солитона в нелинейном одномерном кристалле. Выпускная квалификационная работа, pages 13–18, 2016.</ref>. Однако увеличение <math>\sigma_v^2</math>, как было показано ранее, приводит к росту скорости убывания энергии волны и, как следствие, качественно влияет на процесс развала (рисунок 19, справа).
| |
| − | | |
| − | [[File:StarobinskiiThesis17.png|framed|center|Рисунок 18. Скорость механической волны на момент <math>\tau~=~115.5</math> для <math>\sigma_v^2~=~0</math> (слева) и <math>\sigma_v^2~=~7</math> (справа)]]
| |
| − | | |
| − | Нагляднее это влияние можно рассмотреть на рисунке 18, где приведено развитие скоростей механической волны для разных значений дисперсии по состоянию на один и тот же момент времени. При <math>\sigma_v^2 = 5</math> побочный колебательный процесс выражен неявно, а при <math>\sigma_v^2 = 20</math> вовсе отсутствует. Объясняется это тем, что затухание волны происходит слишком быстро, и наклон графика не достигает критического значения.
| |
| − | | |
| − | При достаточно <math>\sigma_{v}^{2}</math> амплитуда волны продолжает падать, но форма остаётся практически неизменной. На рисунке 17 видно, что график скоростей в кристалле со значением <math>\sigma_v^2~=~7</math> и при <math>\tau~=~115.5</math> по форме близок к графику начального возмущения (см. рисунок 1).
| |
| − | | |
| − | [[File:StarobinskiiThesis18.png|framed|center|Рисунок 19. Скорость механической волны на момент <math>\tau = 27.0</math> для разных значений <math>\sigma_v^2</math>]]
| |
| − | | |
| − | | |
| − | [[File:StarobinskiiThesis19.png|framed|center|Рисунок 20. Эволюция скорости механической волны для <math>\sigma_v^2~=~0</math> (слева) и <math>\sigma_v^2~=~1</math> (справа)]]
| |
| | | | |
| | == Заключение == | | == Заключение == |
| Строка 765: |
Строка 41: |
| | | | |
| | Показано, что энергия стоячей волны убывает в 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: |
Строка 48: |
| | Результаты работы были получены с использованием вычислительных ресурсов суперкомпьютерного центра Санкт-Петербургского политехнического университета Петра Великого. Каждый расчёт производился на отдельном ядре 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> |
