CC BY
58
УДК 539.3: 539.4
Распространение тепла в одномерном гармоническом кристалле на упругом основании
A.M. Кривцов1,2, М.Б. Бабенков1,2, Д.В. Цветков1
1 Санкт-Петербургский политехнический университет Петра Великого, Санкт-Петербург, 195251, Россия 2 Институт проблем машиноведения РАН, Санкт-Петербург, 199178, Россия
В работе получена замкнутая система дифференциально-разностных уравнений, описывающая тепловые процессы в одномерном гармоническом кристалле на упругом основании. Показано, что процесс эволюции теплового возмущения в таком кристалле описывается нестационарным дискретным уравнением, частным случаем которого является гиперболическое уравнение баллистической теплопроводности. Данное уравнение остается справедливым и для случая отрицательных значений жесткости связей между частицами кристалла во всем диапазоне его устойчивости. Фронт теплового возмущения распространяется с максимальной групповой скоростью механических волн. В работе показано, что распространение короткого теплового возмущения в кристалле на упругом основании определяется уравнением баллистической теплопроводности того же вида, что и в кристалле без упругого основания. Единственным параметром данного уравнения является величина максимальной по модулю групповой скорости, т.е. максимальная скорость, с которой в кристалле на упругом основании может распространяться энергия. Данная величина пропорциональна модулю полуразности верхней и нижней частот отсечек. Скорость распространения тепловой волны в кристалле на упругом основании с положительной жесткостью всегда меньше, чем скорость тепловой волны в кристалле без упругого основания. Установлено, что полученное уравнение справедливо как для положительных значений жесткости, так и для отрицательных значений, для которых выполняется условие устойчивости цепочки. Для иллюстрации получено точное решение динамической задачи распространения тепла для параболического начального профиля температуры, моделирующее нагрев одномерного кристалла на подложке коротким лазерным импульсом. За счет наличия дисперсии механических волн в цепочке на подложке их групповая скорость зависит от волнового числа и соотношения жесткостей связей в цепочке и упругого основания. Тепловой фронт распространяется с максимально возможной в системе групповой скоростью, которая зависит только от данного соотношения.
Ключевые слова: одномерный кристалл, теплопроводность, упругое основание, отрицательный коэффициент жесткости, групповая скорость, ковариации
DOI 10.24411/1683-805X-2019-12006
Heat propagation in a one-dimensional harmonic crystal on an elastic foundation
A.M. Krivtsov12, M.B. Babenkov12, and D.V. Tsvetkov1
1 Peter the Great St. Petersburg Polytechnic University, St. Petersburg, 195251, Russia 2 Institute for Problems in Mechanical Engineering RAS, St. Petersburg, 199178, Russia
A closed system of differential equations has been derived to describe thermal processes in a one-dimensional harmonic crystal on an elastic foundation. It is shown that the evolution of thermal perturbation in such a crystal is described by a discrete unsteady-state equation, a special case of which is the hyperbolic equation of ballistic heat conduction. This equation remains valid with negative stiffness of bonds between particles of the crystal in its entire stability range. The thermal perturbation front propagates with the maximum group velocity of mechanical waves. The propagation of a short-term thermal perturbation in the crystal on the elastic foundation is determined by the equation of ballistic thermal conductivity of the same type as in the crystal without an elastic foundation. The only parameter of this equation is the maximum group velocity (in absolute value), i.e., the maximum rate of energy propagation in the crystal on the elastic foundation. This quantity is proportional to the absolute value of the half-difference of the upper and lower cutoff frequencies. The rate of heat wave propagation in the crystal on the elastic foundation with positive stiffness is always lower than that in the crystal without an elastic foundation. The obtained equation is found to be valid both for positive stiffness values and for negative ones, for which the chain stability condition is satisfied. As an example, a dynamic problem of heat distribution is solved exactly for a parabolic initial temperature profile to model heating of a one-dimensional crystal on a foundation by a short laser pulse. Due to the dispersion of mechanical waves in the chain on the foundation, their group velocity depends on the wave number and the ratio of bond stiffnesses in the chain and the elastic foundation. The thermal front propagates with the maximum possible group velocity in the system, which depends only on this ratio.
Keywords: one-dimensional crystal, thermal conductivity, elastic foundation, negative stiffness coefficient, group velocity, covariance
© Кривцов A.M., Бабенков М.Б., Цветков Д.В., 2019
1. Введение
Современные технологии позволяют получать образцы одномерного кристалла углерода (карбина), состоящего из более чем 6 тысяч атомов углерода, собранных внутри двойной углеродной нанотрубки [1]. Одномерные кристаллы, состоящие из атомов различных веществ (так называемые «ионные кристаллы»), могут быть выращены внутри углеродных нанотрубок, обеспечивающих их устойчивость [2]. Развитие технологий, позволяющих получать одномерные структуры, ставит вопрос о физических свойствах синтезированных объектов. Согласно теоретически предсказанным свойствам [3], данные материалы являются отличными проводниками тепла и могут быть полезны, например, для контактного отвода тепла в MEMS/NEMS и в молекулярной электронике, в случаях когда естественное конвекционное и лучистое охлаждение затруднено или невозможно.
Особенностью синтеза устойчивых одномерных структур, обусловленной их высокой химической активностью, является необходимость использования молекулярного «термоса», роль которого выполняет двойная углеродная нанотурбка [2]. Моделью системы «углеродная нанотрубка - одномерный кристалл» в линейном приближении может служить одномерная цепочка масс с гармоническим потенциалом взаимодействия частиц между собой и с закрепленным основанием. При линеаризации сложного потенциала взаимодействия между частицами вещества, находящимися в поле потенциала подложки, жесткость межатомной связи может оказаться отрицательной [4] и в ряде случаев несимметричной [5]. Полученное в работе уравнение переноса тепла применимо как для положительных, так и для отрицательных жесткостей связи при выполнении условий устойчивости такой системы.
Дискретность атомарного строения вещества существенно влияет на волновые процессы в нем [6, 7], внося дисперсию в закон распространения волн. Наличие свободных электронов теплопроводности в решетке проводников создает фотоакустический эффект, влияющий на процесс распространения механических волн [8]. Механические волны, распространяющиеся в системе с потенциалом подложки, имеют фильтрующий закон дисперсии, характеризующийся нижней и верхней частотами среза (отсечки). Сходные эффекты возникают в наноструктурах, в частности в массиве параллельных друг другу наноразмерных осцилляторов. В работах [9, 10] было установлено, что закон дисперсии волн для такого массива является разрывной функцией частоты. К особенностям термомеханических процессов в дискретных системах относятся негативное тепловое расширение [11, 12] и возможность структурных переходов [13, 14]. Подобное разнообразие тепловых
явлений определяет высокую степень сложности описания теплопереноса в дискретных средах.
Результаты последних экспериментальных исследований [15-17] показали, что закон теплопроводности Фурье нарушается в низкоразмерных наноструктурах. Указанные аномалии наиболее ярко проявляются в простейших решеточных моделях, в частности одномерных гармонических кристаллах [ 18-20]. Аналитические и компьютерные исследования [21-24] показывают существенные отклонения от классической теплопроводности Фурье в идеальном монокристалле. Стоит отметить, что эти отклонения, в принципе, могут быть уменьшены или даже полностью устранены использованием специальных законов взаимодействия [25-28] или достаточно сложных структур [29, 30].
По сравнению с принятыми моделями волнового переноса тепла на нано- и микромасштабном уровне в сплошных средах [31-35], протекание тепловых процессов в дискретных средах существенно отличается. В частности, при резком тепловом воздействии на дискретную систему в ней инициируются высокочастотные колебания [36], связанные с выравниванием значений потенциальной и кинетической энергии, согласно теореме о вириале [37]. Точное решение, описывающее протекание такого процесса в гармоническом кристалле, получено в работе [38, 39]. В работах [40, 41] получено фундаментальное решение через полиномы Чебышева и функции Бесселя и исследовано распределение температуры внутри одномерного гармонического кристалла с учетом корреляций, связывающих положение частиц. В работе [42] получено решение для высокочастотных колебаний энергии в гармоническом кристалле на упругом основании и построены его асимптотики для случаев жесткой и мягкой подложки. В работе [43] исследованы высокочастотные тепловые процессы в векторных решетках. В [44] получено общее решение для задачи распространения тепла в Ш и 2D случаях для скалярных решеток, исследованы медленные и быстрые колебания, найдена скорость фронта теплового импульса.
В данной работе исследуется медленный процесс распространения тепла в кристалле, реализующийся после затухания быстрого переходного процесса [42, 45, 46]. Кристалл рассматривается как стохастическая система, случайность в которую вводится через начальные условия, а детерминированные динамические уравнения получаются для ковариаций скоростей и перемещений частиц. Подход, излагаемый в данной статье, основан на работах [38, 45, 47, 48] и позволяет получить достаточно простое динамическое уравнение переноса тепла в частных производных, которое содержит информацию о дисперсии и характере затухания тепловых волн в одномерном кристалле на упругом основании.
2. Уравнения динамики цепочки масс
Если в начальный момент времени частицы одномерной цепочки на упругом основании покоятся, а их скорости задаются случайно, то в решетке начнется процесс перехода кинетической энергии в потенциальную энергию деформирования пружин. Возникнут высокочастотные колебания кинетической и потенциальной энергий, характерные для дискретных сред [36]. Этот процесс будет протекать, пока энергия не распределится, согласно теореме о вириале [37], равномерно между кинетическими и деформационными степенями свободы. Переходный процесс выравнивания потенциальной и кинетической энергии носит высокочастотный характер и затухает за несколько десятков периодов, причем чем жестче подложка, тем медленнее происходит затухание [42]. По окончании переходного процесса в кристалле температура дискретной системы может быть связана с кинетической энергией [49].
Рассмотрим цепочку, состоящую из одинаковых масс т, связанных пружинами жесткости С0. Цепочка находится на упругом основании жесткости С1. Тогда уравнение динамики частиц цепочки имеет вид
(1)
Un = ®0 (un-1 - 0un + Un+1) -
®0 =
n
^ О =лР,
т V т
где ип — перемещения п-й частицы; п — индекс, принимающий произвольные целые значения; С0 — жесткость связей между массами; С1 — жесткость связей между массами и основанием.
Цепочка может быть устойчива и при отрицательных значениях жесткости, если выполняется условие [50]: С0 > -С[/4. В таком случае удобна альтернативная форма записи уравнения (1), имеющая вещественные коэффициенты на всем интервале допустимых значений жесткости С0 и С1:
ип = 1°2 (ип-1 - 2ип + ип+1) - 1®12 (ип-1 + 2ип + ип+1 )>
(2)
ю0 =
Ci + 4 Со
где величины ю1 < т2 являются нижней и верхней частотами отсечек и связаны с парциальной частотой т0
2 2 2 соотношением 4<В2 = ю2 - о^.
Пусть выполняются условия периодичности ик+N + ик, где N >> 1 — число независимых частиц. Определив разностный оператор второго порядка дп следующим образом:
def
Д/ = /п-1 - 2/п + /п+1, (3)
можно записать уравнение (1) в виде
Un = (ЮоА0 -ff)10)un.
(4)
Будем рассматривать начальные условия, соответствующие мгновенному тепловому возмущению в результате воздействия на кристалл ультракороткого лазерного
импульса:
ип \1 =0 = 0 ип\( =0 = СТ(х )Рп > (5)
где рп — независимые случайные величины с нулевым математическим ожиданием и единичной дисперсией; ст(х) — детерминированная девиация начальных скоростей. Интересным с практической точки зрения является случай, когда ст(х) является медленно изменяющейся функцией аргумента х.
3. Вывод уравнений для ковариационной температуры
Для получения уравнений распространения тепла воспользуемся подходом на основании корреляционного анализа [45, 46, 48, 49]. Введя линейный оператор L, перепишем уравнение движения (1) в виде
Un = LnUn > Ln =f Ю°А°
(6)
Очевидно, такому же уравнению удовлетворяет и
def ..
ъп = ип — скорость частицы.
Будем считать, что начальные значения перемещений и скоростей — случайные величины, имеющие нулевое математическое ожидание. Рассмотрим кова-риации перемещений и скоростей частиц цепочки:
^ def , . def , .
Ipq = \U pUq ), Пpq = <Vp Vq ),
(7)
где треугольными скобками обозначено математическое ожидание. Непосредственное дифференцирование соотношений (7) дает
£pq = (Кр + ^ ^pq + 2Пpq ' (8)
Пpq = р + Lq )Пpq + 2LpLq^pq •
Исключение п рц из системы приводит к уравнению 4-го порядка:
^рч - 2(LP + 4 рч +
't^>pq q ' t^pq
+ (Lp - Lq )0 Ipq = 0.
"р ^рч (9)
Отметим, что этому же уравнению удовлетворяют и ковариации скоростей прч • Для операторов Кр и Ьч выполняются соотношения
+ Кч = ю2(Др +Д2) - 2о>2,
Lp - Lq =и°(Ар -А0).
(10)
Кр - =о0(Д р -дч )
Для того чтобы перейти к длинноволновому приближению, введем пространственную координату х и ковариационный индекс и:
(11)
def p + q def
x = —о— a, n = q - p,
где а — шаг решетки кристалла. Будем считать ковариации (7) функциями х и и:
1рч =£п(х), прч =Пп(х), (12)
причем зависимость от пространственной координаты будем считать достаточно гладкой, чтобы функции допускали разложение в ряд Тейлора. Имеем
Ар£РЧ = £»-11 * + 2 ]■-2£»(*) + £»+1 ( *-1'
Г а 1 Г а ] (13)
= £»+11 *+-2 I- 2£» (*) + £»-11 * --2 1 •
Раскладывая ковариации в ряд по а до слагаемых 2-го порядка малости, получаем
АР £ рч = А?»£» + -2(£»_1 -£»+1)' + а 2
+ ~(£»-1 +£»+1) ,
а (14)
а2 £ м = А»£»+а2(£»+1 -£»-1)'+ а 2
+ "8(£»+1 + £»-1) ,
где штрихом обозначена производная по х. Складывая и вычитая эти соотношения, получаем
АР + А^ = 2А» + ^ (А2 + 2)Э 2,
4 (15)
АР -А2 = -а (А»X» )Э*, где использовано
(А2 + 2)£» = £»+1 + £»-1, (А»X»)£» = £»+1 - £»-1. (16) Первое соотношение является очевидным тождеством, второе для простоты можно считать обозначением. В результате, операторы (10) принимают вид
Lp + Lq = 2(и§А2 -ю?) +1 с?(А? + 2)Э2,
4 (17)
Lp - ^ = -юОс(а»х» )д*,
где с ^ ю0а — скорость звука. Подставив соотношение (17) в (9), получим континуальное уравнение для ко-вариаций:
Э4£» -4(ю?А2 -ю?)д?£» -Iс2(А? + 2)Э?£: +
+ ю00с2(А»Х»)? £'я=0, (18)
где точками обозначена частная производная по времени t. Будем также использовать операторную запись уравнения (18): Л£» = 0, где дифференциально-разностный оператор Л определен следующим образом:
Лй= Э4 -4(Ю0А? -Ю?)Э? -?с2(А? + 2)Э?д? +
+ ю0с2(А»Х»)?д?. (19)
Полученное уравнение представляется трудным для анализа, т.к. содержит четвертую производную по времени, вторую производную по координате и разностный оператор А?. Согласно проведенным исследованиям [36], в решеточных структурах наблюдаются два типа динамических процессов: быстро затухающие колебания без переноса энергии вдоль координаты х, возникающие сразу после приложения теплового воздействия, и распространение тепловых волн. Предпола-
гается, что волновой фронт формируется за время затухания быстрых колебательных процессов, необходимое для выравнивания значений потенциальной и кинетической энергии в кристалле [42]. Удерживая слагаемые одного порядка малости в операторе (19), удается разделить (18) на уравнения для быстрых и медленных движений.
3.1. Быстрые движения
Для быстрых движений выполняется соотношение порядков дифференциальных операторов д1 ~ ю, сд* ~ ~ О, где О<<ю ~ ю0 ~ ю1. Сохраняя в операторе (19) только старшие слагаемые порядка ю4, получаем
Л = д? -4(ю?А? -ю12)д?, (20)
что дает для £» уравнение
£» -4(ю?А?-Ю?)£» = А + Лги (21)
Отметим, что полученное уравнение не содержит производных по х, т.е. для быстрых движений пространственный перенос отсутствует. Если константы интегрирования А1 и Л? равны нулю, то уравнение (21) совпадает с полученным ранее уравнением колебаний энергий для одномерного гармонического кристалла на упругом основании [42].
3.2. Медленные движения
Для медленных движений д{ ~ сд* ~ О, где О<<ю~
~ ю0 ~ ю1. Сохраняя в операторе (19) только старшие
2 2
слагаемые порядка ю О , получаем
Л = -4(ю?А» -ю?)д? + ю?с?(А»X»)?д?2, (22)
что дает для £» уравнение
4(ю?2А» -Ю1?)£» = ю??с2(А»X»)?£». (23)
При отсутствии упругого основания (ю1 = 0) уравнение (23)принимает вид
£» -1 с2Х?£» = В1 + В?». (24)
Если константы интегрирования В1 и В? равны нулю, то уравнение (21) совпадает с полученным ранее уравнением динамики ковариаций для одномерного гармонического кристалла на упругом основании [48]. Используя обозначение X = ш?/ю0 + 2, запишем в явном виде полученное уравнение медленных движений (23) для ковариации скоростей п». Как было замечено выше, оно имеет тот же вид, что и уравнение для ковариа-ции перемещений £»:
(П»+1 -ХП» +П»-1)" = 1 с2(п»+?-2П» + П»-2)^' (25)
Данное уравнение полностью описывает макроскопический процесс распространения тепла в одномерном гармоническом кристалле на упругом основании. Уравнение содержит два параметра — скорость звука с и безразмерный параметр X, определяемый отношением жесткостей: Х = С1/ С0 + 2, где С1 — жесткость упру-
гого основания; С0 — жесткость связи между частицами кристалла. Выполняя замену кв(-1)п 6п: = тг|и [50], придем к уравнению относительно ковариационной температуры 6п:
(9й+1 +Л6„ + 6Я_1)" = -4С2(6И+2 - 26и + 6И_2)'. (26)
При допустимых значениях жесткостей _ С^/4 < С0 < и 0< С1 параметр X принимает значения из интервалов < X < -2 и 2<Л<+^. Начальные условия для ковариационной температуры имеют вид То(х), п = 0,
Дx t)|(=0 =
0, 0 < n < N,
(27)
6 п (х, г )\,=0 = 0.
Решение задачи (26), (27), полученное с помощью преобразования Фурье [51, 52] по ковариационной координате п и позиционной координате х, в (п, к, ¿^представлении выглядит следующим образом:
6 п (к, г) = %(к)! х
п
xj cos
0
(
aC0kt sin 8
A
•y/am(1 + 4C0/C1cos2 (8/2))
cos (8n)d8, (28)
где 8 — переменная интегрирования. Условие вещественности корня в формуле (28) дает ограничение на множество допустимых значений жесткостей С0 > > -^/4 устойчивой цепочки [50].
Рассматривая только некоррелированные колебания, положим п = 0. Тогда интеграл (28) может быть представлен через функции Бесселя (доказательство этого утверждения приводится в приложении А):
60 (к, г) = Т0(к) JоC*кt, (29)
где с* = а/2 \ ю2 -ю1 \ — максимальная групповая скорость механических волн в цепочке на подложке.
Если жесткость подложки много больше жесткости цепочки и положительна (С1 >> С0), то с* ~ Су]С0/С1, напротив, если С1 << С0, то с* = с. В случае равенства жесткостей С1 = С0, скорость с* = с/2(>/5 -1) = 0.62с. При С1 = 0 решение переходит в полученное ранее [45]. При отрицательных значениях жесткости С0 скорость распространения тепловых волн с* остается вещественной, пока выполняется условие С0 >-С^/4.
Из представления (29) и свойств функции Бесселя следует, что распределение кинетической температуры Т(к, г) = 60(к, г) в цепочке на подложке удовлетворяет следующей задаче:
' & (30)
T + = -¿к2Т, T\t=0 = T0 (x), T\t=0 = 0.
Обратное преобразование Фурье от (30) по к дает уравнение [45, 53]:
Т + 1т = с*2Г, Т =0 = Т (х), Т\г=0 = 0, (31)
где £ — физическое время, отсчитываемое с момента
воздействия на цепочку коротким тепловым импульсом. Скорость теплового фронта в цепочке с* равна максимальной по модулю групповой скорости механических волн в рассматриваемой системе (см. приложение А). Заключение о скорости теплового фронта можно также сделать, анализируя разложение функции Грина для цепочки (формула (26) в [52]), построенное методом стационарной фазы.
Отметим, что уравнение (31) аналогично уравнению Дарбу [54, 55]
д2Т(х, г) а-1 дТ(х, г)
- = АxT(x, r),
(32)
дг г дг
описывающему сферическое среднее решение волнового уравнения в пространстве размерности а +1 при осреднении по сфере радиуса г.
Решение начальной задачи (31) можно представить [54-56] как в виде свертки:
Т(х, г) =11 Т0(х-с*т)н( -Т )dт =
1 t
= 1J
Tq( x - c*T)
Vt2—2
dT,
так и в виде суммы прямой и обратной волн:
t (x, t)=-Lj
2п 0
Tq| x-ctsin8 | + Tq| x + ctsin
T0
2п.
J T01 x +1 sin — | d8.
(33)
d—=
(34)
Аналогичное выражение получено в работе [44]. Выполняя замену T = t sin (8/ 2) в последнем выражении, можно привести его к виду (34). Для иллюстрации решения задачи о переносе тепла (30) предположим, что бесконечный одномерный кристалл в начальный момент времени был нагрет воздействием короткого лазерного импульса, интенсивность которого убывает с удалением от центра пятна. Будем приближенно считать начальное распределение температуры параболическим:
Т|(=о = To(l2 -х2)H(х2 -12), (35)
где T0 — температура в центре лазерного пятна; l — радиус лазерного пучка; H — функция Хевисайда. Тогда, подставляя (35) в (33), получим функцию, описывающую эволюцию температурного профиля в одномерном кристалле с течением времени (рис. 1, сплошные линии). Вид решения T(x, t) приведен в приложении Б. Анализ решения задачи (30) при некоторых других начальных условиях выполнен в работе [56]. Задачи с подводом тепла в объем исследовались в работе [57].
Система дискретных уравнений (1) динамики цепочки со стохастическими начальными условиями (5), где
ст(х) = T0(l2 - х2)H(х2 -12),
ТШ
0.8-
0.4-
0.01
-0.4 -0.2
0.0
0.2
0.4
Рис. 1. Распределение температуры в бесконечном одномерном кристалле: толстыми линиями показано решение начальной задачи Коши, полученное с помощью формулы (33) и начальных условий (35), пунктирными линиями показано решение, построенное в рамках модели классической теплопроводности Фурье. Точками показано численное решение. Кривые 1-3 соответствуют распределению температуры в моменты времени ^ < ¿2 <
решалась численно методом центральных разностей с шагом интегрирования 5 -10-4 т0, где т0 = 2п/ш0. Численное решение показано на рис. 1 точками. Начальные скорости задавались генератором случайных чисел с равномерным распределением. Математическое ожидание вычислялось как среднее по кристаллу. Для получения численного решения требуемой точности использовалась цепочка, содержащая 50 реализаций по 10000 частиц каждая.
Поскольку единственной константой в уравнении распространения тепла является групповая скорость, то интерес представляет то, как она зависит от параметров системы (1). Для этого в следующем разделе выполним анализ дисперсионных соотношений, в том числе в отрицательной области значений жесткости С0.
4. Дисперсионные соотношения
Дисперсионное соотношение для волн, распространяющихся в цепочке на подложке (1), имеет вид г акл
Т
Q2(k) = (ш2 -w2)sin2
+ Ш
1 >
(36)
ш2 = 4^2 + Ш!2.
Дисперсионное соотношение (36) содержит две частоты отсечки ш1 < ш2, таким образом цепочка на подложке является полосным фильтром. Волны с частотами, лежащими за пределами полосы пропускания, называются экспоненциальными (зигзагообразными) и не распространяются [58]. Для анализа дисперсионных зависимостей как при положительных, так и при отрицательных значениях жесткости С0 оказывается удобным перейти к параметрам ш и в, заданным таким образом, что
2 def
Ш1 = Ш cos
р+-
2 def
Ш2 = Ш Sin
/
V
Ш =Ш1 +Ш2.
(37)
Введенный параметр в зависит от жесткости подложки и жесткости связей между частицами в цепочке следующим образом:
tg
г\
р+-
def Ш2 Ш
C1 + 4Cq C
(38)
и1 \ С1
При отрицательных значениях жесткости С0 параметр в также принимает отрицательные значения. Значениям в из промежутка -я/4 <Р<я/4 соответствуют частоты ш1, ш2, меняющиеся в интервале (0, ш), или жесткости С1 и С2, изменяющиеся в пределах - С^4 < С0 < и 0 < С1 < га. Соответствующие значения величин представлены в таблице. Использование параметра в позволяет получить симметричные графики зависимостей величин, характеризующих распространение волн в системе (1) при любых допустимых значениях жесткости С0 и С1. При этом параметр в остается вещественным и не обращается в бесконечность на всем множестве своих значений.
Дисперсионная характеристика волн в цепочке на упругом основании (36) выражается через параметры в и ш следующим образом: ш
Q(k) = - sin (2P)cos(ak). v2
(39)
На рис. 2 изображено семейство дисперсионных кривых, соответствующих различным значениям параметра в, видно что при положительных и отрицательных в графики расположены симметрично. Таким образом, благодаря выбору способа параметризации, ясно просматривается аналогия в поведении системы (1) при отрицательных - С^4 < С0 < 0 и положительных значениях 0 < С0 < га жесткости цепочки. При в = П4 подложка отсутствует, при в = 0 жесткость С0 = 0 обращается в ноль, что соответствует системе несвязанных осцилляторов. Значению в = -п/ 4 соответствует предельная величина отрицательной жесткости С0 = = - С1 /4, при которой цепочка на подложке сохраняет устойчивость. Все графики дисперсионных зависимостей на рис. 2 ограничены снизу и сверху частотами
Таблица
Связь введенных параметров в и ш с параметрами ш1, ш2, ш0, с* и коэффициентами С0, С1. Величина ш входит в дисперсионные соотношения как множитель
e -П 4 0 -П 4
ш1 Ш V^/2 Ш 0
Ш2 0 42/2 ш Ш
Шо - i/ 2 ш 0 i/ 2 ш
c* a/2 Ш 0 a/2 ш
Cq - Q/4 0 ^
Ci -4Cq - 0
с0>о
О тг/4 ак/2
Рис. 2. Дисперсионная характеристика волн в цепочке на упругом основании при различных значениях параметра в, который меняется в пределах от 0 до п/4 с шагом п/16, для положительных С0 > 0 (а) и отрицательных С0 < 0 значений жесткости (б)
О 71/4 ак/2
Рис. 3. Групповая скорость волн в цепочке на упругом основании при различных значениях параметра в, который меняется в пределах от 0 до п/ 4 с шагом п/16, для положительных С0 > 0 (а) и отрицательных С0 < 0 значений жесткости (б)
отсечек:
Ш1 = w sin
р+-
w = w cos
(40)
Зависимость групповой скорости от волнового числа примет вид
def d^ a ш sin (2Р )sin (ak)
(41)
dk 2a/2^1 - sin (2P) cos (ak)' Групповая скорость волн в цепочке на подложке всегда меньше, чем в цепочке без подложки, которой соответствует график при Р = п/4 на рис. 3, а. При положительных и отрицательных значениях жесткости графики групповой скорости расположены антисимметрично. При С0 > 0 групповая скорость положительна, при С0 < 0 отрицательна и направлена противоположно фазовой. Максимальное по модулю значение групповой скорости в цепочке на подложке определяет максимальную скорость распространения энергии в системе, которая составляет
c = a |ш2-ш11 * 2 ' Величина c* пропорциональна разности верхней и нижней частот отсечек. Переходя к параметрам в и ш запишем выражение для скорости: ato . . П|
c = — |srnB |. (42)
V2
В случае отсутствия подложки, при = 0 и, следовательно, Р = п/ 4, скорость с* принимает значение с = = ato0.
Параметр to при различных допустимых значениях жесткости -С1/4<С0 и 0<С1 меняется в пределах - 14С1 /т < to < ^ и входит в выражения для дисперсионной зависимости (39), групповой скорости (41), а также максимальной по модулю групповой скорости волн (42) только как масштабный множитель. Таким образом интерес представляет только зависимость перечисленных величин от параметра в-
На рис. 4 показана зависимость максимальной по модулю групповой скорости (42) от параметра в. При отрицательных значениях параметра -я/4 <в< 0, соответствующих отрицательным значениям жесткости - ^/4 < С0 < 0, график зависимости ведет себя симметрично интервалу 0<в<П4, соответствующему значениям положительной жесткости 0 < С0 < Видно, что в крайних точках графика в = -п/4, п/4 (т.е. при С0 = -^/4 и С0/С1 ^ ^ соответственно) скорость с* принимает свое наибольшее значение, равное скорости распространения длинных волн в цепочке без подложки: с = а 0)0. В случае в = 0, соответствующем с0 = 0, волны в цепочке не распространяются: с* = 0.
Принятый в работе способ параметризации дисперсионных соотношений позволяет выделить две вели-
Рис. 4. Зависимость максимальной по модулю групповой скорости от параметра в
чины, одна из которых, ш, входит в соотношения дисперсии как множитель и является масштабным фактором, другая, в, определяет качественное поведение графиков зависимостей й(к, в), С^^, в) и с*(в). Обе величины единственным образом определяются из параметров модели (1), при этом выполняется соответствие между в, ш и другими использованными обозначениями, см. таблицу.
5. Заключение
В работе показано, что распространение короткого теплового возмущения в кристалле на упругом основании определяется уравнением баллистической теплопроводности того же вида, что и в кристалле без упругого основания [45]:
Т + -Т = с1т", П=0 = 50 (х), П=0 = 0,
единственным параметром данного уравнения является величина с* = а | ш2 - ш11 /2 — максимальная по модулю групповая скорость, т.е. максимальная скорость, с которой в кристалле на упругом основании может распространяться энергия. Данная величина пропорциональна модулю полуразности верхней ш2 и нижней ш1 частот отсечек и если жесткость подложки много больше жесткости цепочки (С1 >> С0), то с* ~ с^С0/С1, напротив, если С1 << С0, то с* - с. Скорость распространения тепловой волны с* в кристалле на упругом основании с положительной жесткостью С0 всегда меньше, чем скорость тепловой волны с в кристалле без упругого основания с* < с.
Установлено, что уравнение (31) справедливо как для положительных значений жесткости цепочки, так и для отрицательных, для которых выполняется соотношение С0 >-С^/4. Если жесткость цепочки отрицательна и выполняется условие С0 >-С^4, то скорость тепловых волн остается вещественной величиной и может быть найдена по формуле (42).
Для иллюстрации получено точное решение уравнения (31) для параболического начального профиля температуры, моделирующее нагрев одномерного кристалла на подложке коротким лазерным импульсом.
За счет наличия дисперсии механических волн в цепочке на подложке их групповая скорость зависит от волнового числа к и соотношения жесткостей связей в цепочке и упругом основании в. Тепловой фронт распространяется с максимально возможной в системе групповой скоростью, которая зависит только от в. Таким образом, скорости тепловых и акустических волн при фиксированном значении в могут отличаться. Скорость переноса тепла может превышать скорость акустических волн, если жесткость связи между массами в цепочке положительна С0 > 0, так и, наоборот, может быть меньше скорости акустических волн при С0 < 0 (рис. 3).
Авторы выражают благодарность за всестороннее обсуждение работы А.К. Беляеву, О.В. Гендельману, В.А. Кузькину. Работа выполнена при финансовой поддержке Российского научного фонда (грант № 1811-00201).
Приложение А. Представление решения уравнения переноса кинетической энергии с помощью функций Бесселя
Докажем, что выполняется следующее равенство:
/ Л
1п
-Jcos
aC0kt sin 8
d5 =
^ат (1 + 4 С0/ С^2 V 2)
= Jo(c*kt). (А1)
Пользуясь интегральным представлением функции Бесселя [59], преобразуем правую часть
п
J cOS
0
aCkt sin 8
+ 4 C0/ C1 cos2 (8/ 2))
d8 =
= J cos(c*kt sin (8/ 2))d8.
(A2)
Выполним интегральное преобразование Лапласа правой и левой части равенства (A2) по времени t ^ ш:
2 \-1
a (C0k sin 5)
тш(1 + 4 C0/C1cos2( 8/ 2)) ш
d8
= J
-d8.
(A3)
0 с*2 k 28Ш2(5/ 2) + ш2 Вычислим определенный интеграл по 5 в левой и правой части (А3). В предположении, что с* = а/2х х | ш2 -ш1|, полученное выражение сводится к алгебраическому тождеству, которое после упрощающих преобразований принимает вид
д/ц + 2ш2 - 2^¡c^k
4 2 4
4 + цш2 + ш4 +
+ 7 2(ц + 2ш2) - 2c2 k2 =
= д/ц + 2ш2 + 2^/c4k4 + цш2 + ш4 , ц = c2k(A4)
Справедливость данного утверждения при положительных значениях жесткости С0 > 0 и параметров с, к, to > 0, X > 2 доказывается путем последовательных тождественных преобразований. В случае отрицательных значений - ^ /4 < С0 < 0 и параметра < X < -2 доказательство проводится аналогично. В случае С0 = 0, или, что тоже самое, при X = 2, тепловые и акустические возмущения в системе (1) не распространяются. Таким образом левая и правая части равенства (А1) имеют одинаковое изображение по преобразованию Лапласа, следовательно, их оригиналы тоже совпадают, что и требовалось доказать.
Приложение Б. Решение задачи об эволюции теплового импульса с начальным параболическим профилем
Введем обезразмеренные переменные х* и г* такие,
что
х = с*тх*, г = тг*, (Б1)
где с* — максимальная по модулю групповая скорость; т — масштаб по оси времени. В новых переменных и с учетом начальных условий (35) задача (31) приобретает вид:
1
T + -T = T', T t =0 = Го(С - x2)H(x - 4 ),
(Б2)
T\, =о = 0
где /* = // (с*т). Выполняя интегрирование (33), получим следующее выражение для температуры (для краткости индексы в обозначениях опущены):
Т (х, г) = Т (х, г) + Т (-х, г)+Т2 (х, г)+Т3 (х, г), (Б3) где величины Т1, Т2 и Т3 задаются выражениями
Ti(x, t) =
( , ■ i+x Л . , чл
p(x)l 2arcsm-+n 1+ 2q(x)
4n
x H (l2 - (t + x)2), T2(x, t) = — I p(x)( 2arcsinl—x + arcsinl + x | +
2nl I t t
Л
+ q(-x) + q( x)
H((t -1)2 - x2),
T3 (x, t) = 2p(x)H((t -1)2 - x2)H(l -1),
q(x) = (l - 3x)yjt2 - (l + x)2, p(x) = 2l2 -12 - 2x2. В решении наблюдается диффузия волны: задний фронт отсутствует, существует только передний. Функции T2 (x, t) и T3 (x, t) в различные периоды времени определяют внутреннюю часть распространяющейся волны, Tj(x, t) и Tj(-x, t) — ее внешнюю часть при любых значениях времени.
Литература
1. Shi L., Rohringer P., Suenaga K., Niimi Y., Kotakosk J., Meyer J.C., Peterlik H., Wanko M., Cahangirov S., Rubio A., Lapin Z.J. Confined
linear carbon chains as a route to bulk carbyne // Nature Materials. -2016. - V. 15. - No. 6. - P. 634-639.
2. Senga R., Komsa H. P., Liu Z, Hirose-Takai K., Krasheninnikov A.V., Suenaga K. Atomic structure and dynamic behaviour of truly one-dimensional ionic chains inside carbon nanotubes // Nature Materials. - 2014. - V. 13. - No. 11. - P. 1050-1054.
3. Rieder Z., Lebowitz J.L., Lieb E. Properties of a harmonic crystal in a stationary nonequilibrium state // J. Math. Phys. - 1967. - V. 8. -No. 5. - P. 1073-1078.
4. Gendelman O.V., Savin A.V. Heat conduction in a chain of colliding particles with a stiff repulsive potential // Phys. Rev. E. - 2016. -V. 94. - No. 5. - C. 052137.
5. Savin A.V., Kosevich Y.A. Thermal conductivity of molecular chains with asymmetric potentials of pair interactions // Phys. Rev. E. -2014. - V. 89. - No. 3. - P. 032102.
6. Кривцов A.M., Морозов Н.Ф. О механических характеристиках наноразмерных объектов // ФТТ. - 2002. - Т. 44. - № 12. - C. 21582163.
7. Hoover W.G., Hoover C.G. Simulation and control of chaotic non-equilibrium systems // Advanced Series in Nonlinear Dynamics. -World Sci: 2015. - V. 27. - 324 p.
8. Indeitsev D.A., Osipova E.V. A two-temperature model of optical excitation of acoustic waves in conductors // Dokl. Phys. - 2017. -V. 62. - No. 3. - P. 136-140.
9. Еремеев B.A., Иванова E.A., Морозов Н.Ф. Некоторые задачи наномеханики // Физ. мезомех. - 2013. - Т. 16. - № 4. - C. 6773.- doi 10.24411/1683-805X-2013-00031.
10. Eremeyev V.A., Ivanova E.A., Indeitsev D.A. Wave processes in nanostructures formed by nanotube arrays or nanosize crystals // J. Appl. Mech. Tech. Phys. - 2010. - V. 51. - No. 4. - P. 569-578.
11. Kuzkin V.A. Comment on "negative thermal expansion in single-component systems with isotropic interactions" // J. Phys. Chem. -2014. - V. 118. - No. 41. - P. 9793-9794.
12. Kuzkin V.A., Krivtsov A.M. Nonlinear positive/negative thermal expansion and equations of state of a chain with longitudinal and transverse vibrations // Phys. Stat. Solid. B. - 2015. - V. 252. - No. 7. -P. 1664-1670.
13. Гольдштейн Р.В., Городцов B.A., Лисовенко Д.С. Мезомеханика многослойных углеродных нанотрубок и наноусов // Физ. мезо-мех. - 2008. - Т. 11. - № 6. - C. 25-42.
14. Podolskaya E.A., Panchenko A.Y., Freidin A.B., Krivtsov A.M. Loss of ellipticity and structural transformations in planar simple crystal lattices // Acta Mech. - 2015. - P. 1-17. - doi 10.1007/s00707-015-1424-1.
15. Chang C.W., Okawa D., Garcia H., Majumdar A., Zettl A. Breakdown of Fourier's law in nanotube thermal conductors // Phys. Rev. Lett. -2008. - V. 101. - P. 075903.
16. Xu X., Pereira L.F., Wang Y., Wu J., Zhang K., Zhao X., Bae S., Tinh Bui C., Xie R., Thong J.T., Hong B.H., Loh K.P., Donadio D., LiB., Ozyilmaz B. Length-dependent thermal conductivity in suspended single-layer graphene // Nature Commun. - 2014. - V. 5. -P. 3689.
17. Hsiao T.-K., Huang B.-W., Chang H.-K, Liou S.-C., Chu M.-W., Lee S.-C., Chang C.-W. Micron-scale ballistic thermal conduction and suppressed thermal conductivity in heterogeneously interfaced nano-wires // Phys. Rev. B. - 2015. - V. 91. - P. 035406.
18. Lepri S., Mejia-Monasterio C., Politi A. Nonequilibrium dynamics of a stochastic model of anomalous heat transport // J. Phys. A. Math. Theor. - 2010. - V. 43. - P. 065002.
19. Kannan V., Dhar A., Lebowitz J.L. Nonequilibrium stationary state of a harmonic crystal with alternating masses // Phys. Rev. E. - 2012.-V. 85. - P. 041118.
20. Dhar A., Dandekar R. Heat transport and current fluctuations in harmonic crystals // Physica A. - 2015. - V. 418. - P. 49-64.
21. Bonetto F., Lebowitz J.L., Rey-Bellet L. Fourier's law: A challenge to theorists // Mathematical Physics 2000 / Ed. by A. Fokas et al. -London: Imperial College Press, 2000. - P. 128-150.
22. Lepri S., Livi R., Politi A. Thermal conduction in classical low-dimensional lattices // Phys. Rev. P. - 2003. - V. 377. - P. 1-80.
23. Thermal Transport in Low Dimensions: From Statistical Physics to Nanoscale Heat Transfer / Ed. by S. Lepri. - V. 921. - Springer, 2016.
24. Dhar A. Heat transport in low-dimensional systems // Adv. Phys. -2008. - V. 57. - P. 457-537.
25. Aok K., Kusnezov D. Bulk properties of anharmonic chains in strong thermal gradients: non-equilibrium Ф4 theory // Phys. Lett. A. - 2000. -V. 265. - P. 250-256.
26. Gendelman O.V., Savin A.V. Normal heat conductivity of the one-dimensional lattice with periodic potential // Phys. Rev. Lett. - 2000. -V. 84. - P. 2381-2384.
27. Giardina C., Livi R., Politi A., Vassalli M. Finite thermal conductivity in 1D lattices // Phys. Rev. Lett. - 2000. - V. 84. - P. 2144-2147.
28. Gendelman O. V., Savin A. V. Normal heat conductivity in chains capable of dissociation // Europhys. Lett. - 2014. - V. 106. - P. 34004.
29. Bonetto F., Lebowitz J.L., Lukkarinen J. Fourier's law for a harmonic crystal with self-consistent stochastic reservoirs // J. Stat. Phys. -2004.- V. 116. - P. 783-813.
30. Le-Zakharov A.A., Krivtsov A.M. Molecular dynamics investigation of heat conduction in crystals with defects // Dokl. Phys. - 2008. -V. 53. - P. 261-264.
31. Cattaneo C. A form of heat conduction equation which eliminates the paradox of instantaneous propagation // Compte Rendus. - 1958. -V. 247. - P. 431-433.
32. Babenkov M.B., Ivanova E.A. Analysis of the wave propagation processes in heat transfer problems of the hyperbolic type // Contin. Mech. Thermodyn. - 2014. - V. 26. - No. 4. - P. 483-502. - doi 10.1007/s00161-013-0315-8.
33. Ivanova E.A., Vilchevskaya E.N. Description of Thermal and MicroStructural Processes in Generalized Continua: Zhilin's Method and Its Modifications // Generalized Continua as Models for Materials with Multi-Scale Effects or under Multi-Field Actions / Ed. by H. Altenbach, S. Forest, A.M. Krivtsov. - Berlin: Springer, 2013. -P. 179-197.
34. Ivanova E.A. Description of mechanism of thermal conduction and internal damping by means of two component Cosserat continuum // Acta Mech. - 2014. - V. 225. - No. 3. - P. 757-795.
35. Tzou D.Y. Macro- to Microscale Heat Transfer: The Lagging Behavior. - John Wiley and Sons, 2015. - 566 p.
36. Allen M.P., Tildesley A.K. Computer Simulation of Liquids. - Oxford: Clarendon Press, 1987. - 385 p.
37. Ландау Л.Д., Лифшиц EM. Механика. T. 1. Теоретическая физика. - М.: Физматлит, 2004. - 224 с.
38. Кривцов A.M. Колебания энергий в одномерном кристалле // ДАН. - 2014. - Т. 458. - № 3. - С. 279-281.
39. Klein G., Prigogine I. Sur la mécanique statistique des phénomènes irréversibles III // Physica. - 1953. - V. 19. - No. 1-12. - P. 10531071.
40. Гузев M.A., Дмитриев A.A. Различные формы представления решения одномерной гармонической модели кристалла // Дальневосточный математический журнал. - 2017. - Т. 17. - N° 1. - С. 3047.
41. Гузев M.A., Дмитриев A.A. Осцилляционно-затухающее поведение температуры в кристалле // Дальневосточный математический журнал. - 2017. - Т. 17. - № 2. - С. 170-179.
42. Бабенков М.Б., Кривцов A.M., Цветков Д.В. Колебания энергий в одномерном гармоническом кристалле на упругом основании // Физ. мезомех. - 2016. - Т. 19. - № 1. - С. 60-67.
43. Кузькин В.А., Кривцов A.M. Высокочастотные тепловые процессы в гармонических кристаллах // ДАН. - 2017. - Т. 472. - № 5. -С. 529-533.
44. Kuzkin V.A., Krivtsov A.M. Fast and slow thermal processes in harmonic scalar lattices // J. Phys. Condens. Matter. - 2017. - V. 29. -No. 50. - P. 505401.
45. Кривцов A.M. Распространение тепла в бесконечном одномерном гармоническом кристалле // ДАН. - 2015. - Т. 464. - № 2. - C. 162166.
46. Кривцов A.M. Динамика тепловых процессов в одномерных гармонических кристаллах. Вопросы математической физики и прикладной математики // Материалы семинара, приуроченного к 75-летию проф. Э.А. Троппа. 30 сентября 2015. - СПб.: Физико-технический институт им. А.Ф. Иоффе, 2016. - С. 63-81.
47. Krivtsov A.M. Dynamics of Energy Characteristics in One-Dimensional Crystal // Proc. of XXXIV Summer School "Advanced Problems in Mechanics", St.-Petersburg, Russia, 2006. - P. 274-208.
48. Krivtsov A.M. On unsteady heat conduction in a harmonic crystal // ArXiv:1509.02506, 2015.
49. Кривцов A.M. Колебания энергий в одномерном кристалле // ДАН. - 2014. - Т. 458. - № 3. - С. 279-281.
50. Indeitsev D.A., Sergeev A.D. Correlation between the properties of eigenfrequencies and eigenmodes in a chain of rigid bodies with torque connections // Vestnik St. Petersb. Univ. Math. - 2017. - V. 50. -No. 2. - P. 166-172. - doi 10.3103/S1063454117020066.
51. Слепян Л.И., Яковлев Ю.С. Интегральные преобразования в нестационарных задачах механики. - Л.: Судостроение, 1980. -343 с.
52. Gendelman O.V., Shvartsman R., Madar B., Savin A.V. Nonstationary heat conduction in one-dimensional models with substrate potential // Phys. Rev. E. - 2012. - V. 85. - No. 1. - P. 011105.
53. Babenkov M.B., Krivtsov A.M., Tsvetkov D.V. Unsteady heat conduction processes in a harmonic crystal with a substrate potential // arXiv preprint arXiv: 1802.02037. - 2017.
54. John F. Plane waves and spherical means applied to partial differential equations // Courier Corporation. - 2004 (Originally published in 1955).
55. Courant R., Hilbert D. Methods of Mathematical Physics: Partial Differential Equations. - New York: Interscience, 1962. - V. II. - 830 p.
56. Соколов A.A., Кривцов A.M., Muller W.H. Локализованные тепловые возмущения в одномерном гармоническом кристалле: решения уравнения аномальной теплопроводности // Физ. мезомех. -2017. - Т. 20. - № 3. - С. 63-68.
57. Gavrilov S.N., Krivtsov A.M., Tsvetkov D.V. Heat transfer in a one-dimensional harmonic crystal in a viscous environment subjected to an external heat supply // Continuum Mech. Thermodyn. - 2018. -doi 10.1007/s00161-018-0681-3.
58. Крауфорд Ф. Берклеевский курс физики. Т. 3. Волны. - М.: Наука, 1974.
59. Gradshteyn, Ryzhik's Table of Integrals, Series, and Products / Ed. by D. Zwillinger. - Elsevier, 2014.
Поступила в редакцию 12.11.2018 г., после доработки 12.11.2018 г., принята к публикации 12.12.2018 г.
Сведения об авторах
Кривцов Антон Мирославович, д.ф.-м.н., чл.-корр. РАН, проф. РАН, зав. каф. СПбПУ, зав. лаб. ИПМаш РАН, akrivtsov@mail.ru Бабенков Михаил Борисович, к.ф.-м.н., доц. СПбПУ, нс ИПМаш РАН, mikhail.babenkov@gmail.com Цветков Денис Валерьевич, ассист. СПбПУ, dvtsvetkov@ya.ru
