Об одном методе решения задач о термоупругодинамической неустойчивости скользящего фрикционного контакта Текст научной статьи по специальности «Физика»

Об одном методе решения задач о термоупругодинамической неустойчивости скользящего фрикционного контакта

В.Б. Зеленцов, Б.И. Митрин Донской государственный технический университет, Ростов-на-Дону

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

Ключевые слова: термоупругодинамическая неустойчивость, скользящий фрикционный контакт, покрытие, трение скольжения, динамика, термоупругость, полюсы, математическое моделирование, параметрический анализ, анализ устойчивости

Введение

В машиностроении и обрабатывающей промышленности при расчёте узлов трения машин возникает необходимость расчёта скользящего фрикционного контакта взаимодействующих между собой рабочих поверхностей машин, станков и агрегатов. Это, в первую очередь, связано с температурным саморазогревом от трения взаимодействующих поверхностей и возникновением так называемой термоупругодинамической неустойчивости решения соответствующих контактных задач [1 - 8]. Исторически первыми методами исследования такого рода задач являются методы малых возмущений

[3 - 5], позволяющие установить устойчивость или неустойчивость решения задачи, а также установить параметрическую область устойчивости или неустойчивости решения задачи [9-11]. К другим подходам решения задач о термоупругодинамической неустойчивости скользящего контакта можно

отнести численные методы, разработке которых уделяется всё больше внимания в последнее время [4, 5].

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

В данной работе в качестве упрощенной модели скользящего фрикционного контакта рассматривается одномерная задача о скольжении с постоянной скоростью V жесткой полуплоскости В (к < х < да) по поверхности (х = к) упругого покрытия в виде бесконечной полосы шириной к (0 < х < к), нижняя сторона которого (х = 0) жестко соединена с недеформируемым основанием в виде полуплоскости А (-да< х < 0) (рис. 1) [5].

т

Рис. 1. - Постановка задачи

Скольжение недеформируемой полуплоскости В по поверхности покрытия происходит с учетом кулоновского трения, но без учета износа покрытия. Полуплоскость В деформирует упругое покрытие, смещаясь вдоль оси х по закону и(к, г) = -Л(г), г > 0, где и(к, г) - вертикальное смещение поверхности покрытия при х = к. В начальный момент температура покрытия нулевая: Т(х,0) = Т0 = 0 (0<х<к), где Т(х,г) - функция распределения температуры в покрытии. Движущаяся полуплоскость В теплоизолирована,

а поток тепла К дТг), образующийся за счет трения, направлен в упругое

дх

покрытие. Так как нижняя сторона покрытия лежит на недеформируемом основании (х = 0) в виде полуплоскости А, то на этой стороне покрытия упругие смещения и(0, г) = 0. На нижней стороне покрытия поддерживается нулевая температура Т(0, г) = 0. С учётом того, что до начального момента

< 0) (рис. 1) [5].

жесткая полуплоскость А

и

времени покрытие находилось в покое, начальные условия на и

ди ди (х 0)

и — нулевые: и( х,0) = —= 0. д* дг

Сформулированная одномерная задача о скользящем контакте приводит к следующим граничным условиям:

х = к : и(к,г) = -Л(г), К^ЕУЫ! = -уТо(А,г), г > 0; (1)

дх

х = 0: и(0, г) = 0, Т(0, г) = 0, г > 0, (2) где а(х, г) - напряжения сжатия в покрытии, К - коэффициент

теплопроводности материала покрытия, Л(г) - закон внедрения

полуплоскости В в упругое покрытие, который принимается в следующем виде

- ^ < г < 0,

-1 + вег, 0 < г < гЕ, в> 0, (3)

Л(г) = Л 0

8'

гв < г < да,

где в > 0, гв = в-11п 2, Л0 - глубина внедрения полуплоскости В в упругое покрытие (0 <Л0 < к). Скорость внедрения полуплоскости В определяется как производная Л (г) и выражается формулой

Л (г) = Л 0

0, - да < г < 0,

веег, 0 < г < гв, (4)

0, гв < г < да.

г

Максимальная скорость внедрения достигается при г = гв и составляет

Л(гв) = 2Л0в, при этом начальная скорость внедрения равна Л(0) = у0 =вЛ0, откуда в = у0/ Л 0 (при заданных у0 и Л 0).

Изменение и( х, г) и а( х, г) в полосе описывается одномерным уравнением теории упругости [12]

^-р^и- = 0, 0 < х < к, г > 0, (5)

дх дг2

1

и

а температура Т(х, г) - уравнением теплопроводности [13, 14] с учётом нулевой начальной температуры

д 2Т 1 дТ Л

—----= 0, 0 < х < к, г > 0,

дх к дг

(6)

где р, к - плотность и коэффициент температуропроводности материала покрытия соответственно. Связь между а( х, г), и( х, г) и Т (х, г) устанавливается соотношением Дюамеля - Неймана [14]

а(х,г) = М-^ * - аТ(х,г),

1 - 2у дх 1 - 2у

(7)

где |, V, а - модуль сдвига, коэффициент Пуассона, коэффициент линейного расширения материала покрытия соответственно.

Для удобства построения решения задачи (1) - (6) выражение о( х, г) (7) подставляется в дифференциальное уравнение движения упругой среды (5). В результате этого получим дифференциальное уравнение движения упругой среды относительно и( х, г) и Т (х, г)

д 2и 1 д 2и 1 + V дТ -а-

дх а дг 1 - V дх

0 < х < к, г > 0, а =

V

2|(1 -V)

р(1 - 2v)

(8)

где а - скорость упругой волны.

Начальные условия на и (х, г), Т (х, г) и Л(г) нулевые

и (х,0) = ди(х,0) = 0, Т (х,0) = 0, 0 < х < к, Л(0) = 0. дг

(9)

Таким образом, решение рассматриваемой задачи, с учетом начальных условий (9), сводится к решению дифференциальных уравнений (6), (8) с граничными условиями (1), (2), а а(х,г) в (1) определяется формулой (7).

Решение поставленной задачи

Поставленная нестационарная начально-краевая контактная задача (1)-(9) относится к разряду собственно-связанных задач термоупругости, так как а(х, г) и Т(х, г) связаны не только в формулах (5) и (7), но и во втором

и

граничном условии (1). Её решение может быть построено различными методами математической физики [13]. С помощью интегрального преобразования Лапласа [15] решения поставленной задачи (1)-(9) T(х, t), u( x, t), а( х, t) определяются с помощью контурных квадратур

T (х, t) =

1 - у уУ 1 1 + у ак 2га

- ГБ(г)Ыт (х, г)-1 (г)ехр(г~ )г,

t tк ,

и(х, t) = -2-ГДг)Ыи (х, г) 1 (г)ехр(г~ )г,

а( х, t) = -

2^(1 - у)

Е ^ |л(а-(х, t ))+^\(ш+(х, I ^

+ -

(1 - 2у)к - ГБ(г)Nс(х, г(г)ехр(г~)

ТГ 7 *

(10) (11)

(12)

2га

Лг.

справедливых для 0 < х < к, t > 0, где Г = {г: + , + Лк} - контур интегрирования в комплексной плоскости, представляющий прямую линию, параллельную мнимой оси и отстоящую от неё на величину Лх^, в которой Л выбирается таким образом, чтобы все полюса подынтегральных функций в

(10)-(12) были бы левее Лтк,

Я(г) = (1 - у2г)еЬ л/г бЬ уг - уУ(еЬ л/г еЬ уг - у л/г бЬ л/г бЬ уг -1)

Nт (х, г) = л/г (1 - у2 г )еЬ у г бЬ (л/гхк-), Nu (х, г) = (1 - у2 г)еЬ л/г бЬ (угхк- ) -

- уУ(еи(л/гхк- )еЬ уг - ул/г бЬ л/г 8и(угхк- ) - еЬ(уг(к - х)к-1)), Na (х, г) = у г[[ (х, г) бЬ уг - еЬ(угхк-1 )(г )], N° (х, г) = (1 - у2г)л/гсЬ(угхк-) +

+ уУ (ул/г (бИл/ТеЬ(угхк-)- бЬ (л/гхк- )еЬуг )- вЬ(уг(к - х)к- )), Яа(г )= 8ЬугЯ(г ),

г) = г-1 (2ехр(- * ек)- 1)+ (г - 8tк)-1 (1 - ехр(- (г - е^)«)),

+ / Ч ( - 1)к ± х ш;(х, t ) = t->--,

(13)

(14)

(15)

(16)

(17)

(18)

и

где а - постоянная из (7), а - из (8), безразмерные параметры у и V определяются формулами

-Л V-^а 2[д(1+v)h

У_ ак ' = К 1 - 2у ' Внеинтегральные слагаемые получены в формуле (12) в результате выделения обобщенных квадратур [16 - 18] с помощью выражения

В(z)ch(yzxk~1)/shyz.

Квадратуры в формулах (10)-(12) существуют при выполнении условия г >(к - х)/а и благодаря алгебраическому убыванию подынтегральных функций на бесконечности

В(г)Ыт (х, г)Яг) - О

С _ 1/ Л

: е к V J

при | г |^да, (19)

Ы.(х.-) Я _1( Ю]О Ы„ (х, г )Я»|

( Г _0 к_хЛ ( \ -2 I ¡-п_х 2 I -У г-

г ' I 1 " 2

V,-Iе к,

+ о

-V |е к

V ^ I J

при | г да . (20)

Исследование подынтегральных функций в формулах (10)-(12) показывает, что все они мероморфны в комплексной плоскости переменной интегрирования г -£ +щ, то есть имеют в качестве изолированных особых точек только полюсы, которые доставляются обращением в ноль знаменателей этих функций: Я(г), Яа(г). Кроме того, следует отметить, что

подынтегральные функции в (10)-(12) при г = 0 и г - у-2 обращаются в ноль:

Ыт (х,0) - Ыи (х,0) - Ы„ (х,0) - Я(0) - Яа (0) - 0, (21)

(21)

Ыт (х, у-2) - Ыи (х, у-2) - Ыа (х, у-2) - я(у-2) - Яа (у-2) - 0, причём Ыа (х, г) и Яа(г) имеют двукратный ноль при г = 0.

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

и

Нули функции Я(г) в комплексной плоскости.

Для вычисления контурных интегралов в формулах (10)-(12) методами теории функций комплексного переменного необходимо знание нулей функции Я(г) и их свойств в комплексной плоскости. Нули функции Я(г) из (13) определяются в комплексной плоскости г -£ + щ из решения трансцендентного уравнения

Я( г) - (1 -у2 г)еЬ л/2 бЬ уг - у V (еЬ л/7 еЬ уг - уу[7 бЬ бЬ уг -1) - 0 , (22)

которое совпадает с соответствующим характеристическим уравнением [5].

При численном определении корней уравнения (22) в комплексной плоскости г -£ +щ использовались итерационные численные методы определения корней, требующие хорошего начального приближения. Для получения начального приближения корней уравнения применяется параметрический анализ функции Я(7) по двум безразмерным параметрам у и V. Зафиксировав у и положив V - 0 в (22), получим начальные приближения для нулей (22)

2

п

-(1 + 2к) 2

к - 0,1,2,..., (23)

пп

7оп -±' —, п - 0,1,2,3,. (24)

У

Нули г0к (23) и 2~0п (24), определенные при V - 0 из уравнения (22), являются начальными приближениями для последовательного с увеличением V определения соответствующих по номеру нулей гк -гк(V), к = 0,1,2..., и -2±(у), п = 1,2., из численного решения уравнения (22) при фиксированных у.

На рис. 2 представлено расположение множества полюсов гк - гк(V), к = 0,1,2., и - (V), п = 1,2., определенные численными методами для произвольных значений V - 0 и у-1: множества полюсов гк (V), к = 1,2.,

20к ~

и

представляют собой отрезки на отрицательной части действительной оси, точки множества г0 (У) располагаются на луче, выходящем из точки на отрицательном части действительной оси и совпадающего с положительной частью действительной оси; множества г±(У), п = 1,2..., располагаются в правой полуплоскости в виде незаконченных эллипсов. Стрелками на графиках указаны направления расположения точек при увеличении У .

Следует заметить, что графики г+ и г- обладают симметрией относительно действительной оси 1т(г) = п = 0, как и подынтегральные функции в квадратурах (10)-(12), а для г± выполняется равенство

г±= ^. (25*)

и

На всех траекториях нулей z± (F), n = 1,2,3,..., указаны точками значения zk (Vj) и z± (Vj), где j = 1,2,3, причём Vx = 1, V2 = 2, V3 = 2,7. Значения z± (Vj), n = 1,2,3,..., при одинаковых Vj, j = 1,2,3, соединены тонкими линиями.

Графики нулей функции R(z): zk(F),z±(F), n,k = 1,2,3,..., z0(F), представленных на рис. 2, рассчитаны с помощью программы [19], написанной для системы Maple.

Анализ полученных решений

После вычисления контурных квадратур в комплексной плоскости в (10)-(12) [20, 21] для функций T(x, t), u(x, t) и a(x,t) получаются удобные для

Рис. 2. - Расположение нулей Я( г) в комплексной плоскости г при У е[0, го), у = 1. В левом верхнем углу - общий план расположения нулей, в нижнем

левом углу - расположение множества точек г1(У) в увеличенном формате

и

вычисления формулы в виде рядов по полюсам подынтегральных функций:

T (х,г) -¿^± Sl (х,~) H ((-1)" ( - О),

1 + V а к к-1

и( х, г) - - А 0 ]Г Sk (х, ~) Н ((-1) (~ - г К)),

к-1

0 < х < к г > 0

0 < х < к г > 0

а(х, г) - -

2[(1 - V) А0

1 - 2v к

± га (А (-(х, г ))+А ((х, г)))-

к-1

+

(х,~) + Sk0 (х, ~)) ((-1)(~ - гек)) ,

к-1 _

где Sk(х, г) (к -1,2) выражаются формулами:

^ (х, г) - 51 (х, г) + G(x, г),

0(х, г)-2Яе ¿0(2+): (, х, г )+]Г0(2к )К '(гк, х, г),

0 < х < к, г > 0,

к-1

к - 0

S2(х,г) - (х,г) + о(х,г)- 20(х,г - гек), ^ (х, г) - К(гкв, х, г) - К(0, х,0), (х, г) - К(0, х,0),

К(г, х, г) - Ы(х 7) ехр(гг), К'(г, х, г) - Ы(X, 7) ехр(гг), Я(г) Я'(г) 4 у

К(0, х,0) - 11ш, 0(2) - /кЕ ч, 7^0 Я(г) 2(2 - гкв)

S10 (х, г)-2Яе ¿е^) (+к, х, г),

к-1

S20(х,г)-2S10(х,г - О-^(х,г)

(26)

(27)

(28)

(29)

(30)

(31)

(32)

(33)

где Я'(г) есть производная по г, Н(г) - функция Хэвисайда.

В формуле (28) под А(г) понимается нормированная функция А(()- веегН(( -г)н(() из (4), где в - v0/А0. Функция G(x,г) (30) в формулах для Sk (х, г) (к -1,2) содержит два бесконечных ряда, один из которых по полюсам 2+, к -1,2,3,..., следующего вида:

S+-¿ Ь(х, ь(х, 7 )-0(7 )Щ ,

(34)

к-1

и

сходимость которого не очевидна, так как Яе(г+) > 0, t > 0. С учетом асимптотики г+ при к ^ да можно показать, что последовательность экспонент } является ограниченной для фиксированных значений ~

(0 < ~ < да):

<М (М > 0),

а ряд (34) является сходящимся [22], в том числе за счет убывания 9(г) при |г| ^да. С другой стороны, учитывая упорядоченность {яе(г+)}, где для двух

(35)

соседних членов выполняется условие

Re(z+)> Re(z++1), k = 1,2,3,..., сумму ряда (34) можно представить в виде

У = 4x,z+fl + ]Г в(х, z;)"'-leZ|+~, B(x, z ) =

^ k = 2 J b\x, z 1)

где ряд в скобках является сходящимся, причем со скоростью геометрической прогрессии, так как Re(z+-z')> 0, k = 2,3,... Из выражения (35) вытекает, что при неограниченном возрастании ~ (~ ^да), неограниченно возрастает S' из (34), а решение, его содержащее, является неустойчивым при любых значениях V > 0.

В формуле для G(x, t) (37) из Sk (x, t) (k = 1,2) , кроме вышеупомянутого ряда (34), содержится и другой бесконечный ряд по полюсам zk, k = 0,1,2,3,... Учитывая свойства z0(V), изученные выше, а именно то, что при 0 < V < 2 полюс z0 (V)< 0, при V = 2 полюс z0 (2) = 0, при V > 2 полюс z0 (V)> 0, можно утверждать, что ряд по полюсам zk

S-=£b(x,Zk) , (36)

k=0

и

где Ь(х, 7) из выражения (34) для фиксированных значений ~ сходится со скоростью геометрической прогрессии для любых у и V. Однако при V > 2, когда г0 (V)> 0, сумма ряда S- (36), перестроенная по формуле

S- - Ъ(х, 70(1 + В(х, гк)(к 70 )гV0*, (37)

V к-1 У

где В(х, г) из (35), при ~ ^да неограниченно возрастает, хотя ряд, стоящий в

скобках, сходится. Это означает, что при V > 2 решение, содержащее ряд (36), является неустойчивым. Таким образом, решения рассмотренной задачи (26)-(28), содержащие ряды (34) и (36), являются термоупругонеустой-чивыми, начиная со сколь угодно малой скорости скольжения V > 0.

Численный анализ полученных решений

Для детального изучения термоупругонеустойчивого решения рассматриваемой задачи используются формулы (26)-(28), реализованные в программе для ЭВМ [19]. В качестве материала покрытия рассматривается

9 2

алюминий со следующими свойствами: ц = 25.5 ■ 10 Н/м , V = 0,34, а = 1,04 ■ 10-5 К-1, а = 6,24 ■ 103 м/с, к = 8,74 ■ 10-5 м2/с, К = 209,3 Втм/К. Принимаются следующие значения параметров задачи: / = 0,15, Д0 = 0,1 к,

_3 _7

VI) = 0,01 м/с, к = 2 мм, V = 5 м/с, г8 = 6,9310 с, га = 3,210 с, в результате чего безразмерные параметры задачи у и V приобретают следующие значения: у = 710-6 и V = 0,173.

На рис. 3 представлены графики изменения смещений и(х,г) по толщине покрытия (0 < х < к), рассчитанные по формуле (27) при х = 0,5к. На рис. 3 видно, как со скоростью геометрической прогрессии развивается амплитуда смещений неустойчивого решения и( х, г) на собственной частоте 1т(г+), теряя физический смысл с некоторого момента времени. На врезках, которые показывают тот же график в более мелком масштабе по времени и по величине смещений, наглядно показан рост амплитуды колебаний.

и

Рис. 3. - График смещений и(х,0 посередине покрытия (х = 0,5 к)

На рис. 4 приведены графики изменения напряжений а(к,t) на контакте. Начиная с некоторого момента времени, амплитуда колебаний а(к, Г) на собственной частоте 1т(г+) нарастает по экспоненте и неустойчивое решение теряет физический смысл. Врезки на рис. 4 подробно показывают изменение величины напряжений а(к, t) на различных временных отрезках в более мелком масштабе: различается образование на контакте волны сжатия и запаздывание по времени прихода волны, отраженной от жесткого основания, при котором происходит удвоение амплитуды волны,

Рис. 4. - График напряжений а(к^) на контакте

Выводы

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

начиная с некоторого момента времени, полученное решение теряет физический смысл - амплитуда собственных колебаний неограниченно возрастает по времени. Использованный в работе метод позволил получить эффективные формулы для вычисления решения задачи, которые позволяют исследовать свойства решения с помощью программы для ЭВМ [19].

Литература

1. Barber J.R. Thermoelastic instabilities in the sliding of conforming solids // Proceedings of the Royal Society A. 1969, V. 312, pp. 381-394.

2. Dow T.A., Burton R.A. Thermoelastic instability of sliding contact in the absence of wear // Wear. 1972, V. 19, pp. 315-328.

3. Burton R.A., Nerlikar V., Kilaparti S.R. Thermoelastic instability in a seallike configuration // Wear. 1973, V. 24, pp. 177-188.

4. Barber J. R., Dundurs J., Comninou M. Stability considerations in thermoelastic contact // Transactions ASME. Journal of Applied Mechanics. 1980, V. 47, iss. 4, pp. 871-874.

5. Afferrante L., Ciavarella M., Barber J.R. Sliding thermoelastodynamic instability // Proceedings of the Royal Society A. 2006, V. 462, pp. 2161-2176.

6. Ciavarella M., Johansson L., Afferrante L., Klarbring A., Barber J.R. Interaction of thermal contact resistance and frictional heating in thermoelastic instability // International Journal of Solids and Structures. 2003, V. 40, iss. 21, pp.5583-5597.

7. Moirot F., Nguyen Q.S. Brake squeal: a problem of flutter instability of the steady sliding solution // Archives of Mechanics. 2006, V. 52, pp. 645-662.

8. Kinkaid N.M., O'Reilly O.M., Papadopoulos P. Automotive disk brake squeal // Journal of Sound and Vibration. 2003, V. 267, pp. 105-166.

9. Кривоносов В.А., Митин А.С. Наблюдаемость и управляемость системы стабилизации уровней расплавленного металла на МНЛЗ //

Инженерный

вестник

Дона. 2013, №3. URL:

ivdon.ru/ru/magazine/archive/n3y2013/1831.

10. Красильников А.Я., Кравченко К.Ю. Устойчивость линейных дифференциальных уравнений с постоянным запаздыванием, описывающих процесс фрезерования // Инженерный вестник Дона. 2014, №1. URL: ivdon.ru/ru/magazine/archive/n1y2014/2250.

11. Моров В.Н., Черский И.Н. Термоупругая неустойчивость фрикционного контакта штампов с полупространством // Трение и износ. 1985, Т. 6, № 1. С. 27-38.

12. Лурье А.И. Теория упругости. Москва: Наука, 1979. 979 с.

13. Тихонов А.Н., Самарский А.А. Уравнения математической физики. Москва: Наука, 1977. 735 с.

14. Новацкий В. Вопросы термоупругости : пер. с польского. Москва: Изд-во АН СССР, 1962. 363 с.

15. Диткин В.А., Прудников А.П. Интегральные преобразования и операционные исчисления. Москва: Физматлит, 1961. 524 с.

16. Виленкин Н.Я. и др. Функциональный анализ. Москва: Физматлит, 1961. 524 с.

17. Бейтмен Г., Эрдейи А. Высшие трансцендентные функции. Москва: Наука, 1965. Т. 1. 296 с.

18. Брычков Ю.А., Прудников А.П. Интегральные преобразования обобщенных функций. Москва: Наука, 1977. 288 с.

19. Зеленцов В.Б., Митрин Б.И., Айзикович С.М. Вычисление температуры, смещений, напряжений в задаче о термоупругодинамической неустойчивости упругой полосы при скользящем контакте. Св-во о гос. рег. программы для ЭВМ № 2014661451, заявл. 16.09.14, зарег. 30.10.14.

20. Гурвиц А., Курант Р. Теория функций. Москва: Наука, 1968. 648 с.

21. Титчмарш Е. Теория функций. Москва: Наука, 1980. 464 с.

22. Бриллинджер, Д. Временные ряды. Москва: Мир, 1980. 536 с.

References

1. Barber J.R. Proceedings of the Royal Society A, 1969, vol. 312. pp. 381394.

2. Dow T.A., Burton R.A. Wear, 1972, vol. 19, pp. 315-328.

3. Burton R.A., Nerlikar V., Kilaparti S.R. Wear, 1973, vol. 24. pp. 177-188.

4. Barber J. R., Dundurs J., Comninou M. Transactions ASME. Journal of Applied Mechanics, 1980, vol. 47, iss. 4, pp. 871-874.

5. Afferrante L., Ciavarella M., Barber J.R. Proceedings of the Royal Society A, 2006, vol. 462, pp. 2161-2176.

6. Ciavarella M., Johansson L., Afferrante L., Klarbring A., Barber J.R. International Journal of Solids and Structures, 2003, vol. 40, iss. 21, pp. 55835597.

7. Moirot F., Nguyen Q.S. Archives of Mechanics, 2006, vol. 52, pp. 645-662.

8. Kinkaid N.M., O'Reilly O.M., Papadopoulos P. Journal of Sound and Vibration, 2003, vol. 267, pp. 105-166.

9. Krivonosov V.A., Mitin A.S. Inzenernyj vestnik Dona (Rus), 2013, №3. URL: ivdon.ru/ru/magazine/archive/n3y2013/1831.

10. Krasil'nikov A.Ya., Kravchenko K.Yu. Inzenernyj vestnik Dona (Rus), 2014, №1. URL: ivdon.ru/ru/magazine/archive/n1y2014/2250.

11. Morov V.N., Cherskii I.N. Soviet Journal of Friction and Wear. 1985. Vol. 6, №1. C. 27-38.

12. Lourie A.I. Teoriya uprugosti [Theory of elasticity]. Moscow: Nauka, 1979. 979 p.

13. Tikhonov A.N., Samarskii A.A. Equations of Mathematical Physics. New York: Dover Books on Physics, 2011. 800 p.

14. Nowacki W. Thermoelasticity. Reading: Addison-Wesley Publishing Company, 1962. 628 p.

15. Ditkin V.A., Prudnikov A.P. Integral'nye preobrazovaniya i operatsionnye ischisleniya [Integral transforms and operational calculus]. Moscow: Fizmatlit, 1961. 524 p.

16. Vilenkin N.Ya., Flaherty R.E. Functional Analysis. Groningen: Wolters-Noordhoff B.V., 1972. 379 p.

17. Bateman H., Erdelyi A. Higher transcendental functions. Vol. 1. New York: McGraw-Hill, 1953.

18. Brychkov Yu.A., Prudnikov A.P. Integral Transforms of Generalized Functions. New York: Gordon and Breach Science Publishers, 1989. 343 p.

19. Zelenstov V.B., Mitrin B.I., Aizikovich S.M. Certificate of state registration of computer program № 2014661451, appl. 16.09.14, reg. 30.10.14.

20. Hurwitz A., Courant P. Teoriya funktsiy [The Theory of Functions]. Moscow: Nauka, 1968. 648 p.

21. Titchmarsh E.C. The Theory of Functions (2nd edition). New York: Oxford University Press, 1976. 464 p.

22. Brillinger D.R. Time series: Data analysis and theory. New York: Holt, Rinehart & Winston, 1975. 512 p.