Математическое моделирование задач теории упругости с односторонним дискретным контактом Текст научной статьи по специальности «Физика»

Математика к Математическое

моделирование

Ссылка на статью: // Математика и Математическое моделирование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2015. № 04. С. 93-110.

Б01: 10.7463/шаШш.0415.0801840

Представлена в редакцию: 08.07.2015 Исправлена: 22.07.2015

© МГТУ им. Н.Э. Баумана

УДК 539.37

Математическое моделирование задач теории упругости с односторонним дискретным контактом

Станкевич И. В.1*

ар1тех@уапс1е;ш1

1МГТУ им. Н.Э. Баумана, Москва, Россия

В работе рассматривается алгоритм построения численного решения контактной задачи теории упругости применительно к телу, которое имеет выраженное одностороннее дискретное контактное взаимодействие с абсолютно упругим полупространством. Особенностью рассматриваемого алгоритма является специально разработанная процедура коррекции касательных сил в дискретных контактных точках, позволяющая добиться достаточно точного выполнения принятого закона трения. Алгоритм встроен в общую конечно-элементную технологию, с помощью которой создана прикладная программа. Исследована вычислительная эффективность алгоритма с помощью решения модельных задач.

Ключевые слова: дискретное контактное взаимодействие; контактная задача теории упругости; метод конечных элементов

Введение

Разработка и эксплуатация объектов современной техники и новейшей технологии требуют надежной оценки прочностных характеристик ответственных элементов конструкций и технологического оборудования, находящихся в условиях воздействия высокоинтенсивного термомеханического нагружения, сопровождаемого, как правило, контактным взаимодействием [1, 2, 6 - 11]. Математическое моделирование напряженно-деформированного состояния таких деталей и узлов в зонах контакта, базирующееся на адекватных математических моделях, современных численных методах и эффективных прикладных алгоритмах, реализующих непосредственное определение собственно самих полей перемещений, деформации и напряжений, является основным инструментом, обеспечивающим оперативное получение требуемых данных для расчетов на прочность и долговечность [1, 6, 7, 10, 11].

В данной работе приведена постановка задачи одностороннего дискретного контактного взаимодействия упругого тела и абсолютно жесткого полупространства, описан чис-

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

1. Постановка задачи теории упругости с односторонним дискретным

контактом

Рассмотрим в двухмерном пространстве М2 с декартовой системой координат Ох:х2 тело, занимающее область О, ограниченную кусочно-гладкой границей д О (рис.1). Математическая формулировка квазистатической задачи теории упругости включает следующие соотношения: уравнения равновесия

ал( и ) + 0 ( х) — 0, х еО, (1)

граничные условия (кинематические и силовые)

и (х )|5 = и0 (х), х е ^ ^д О, (2)

(и) п1 = Р (*), х е ^д О (3)

2

соотношения Коши для компонент тензора деформации

Еч(х) = \(и^ П),хеО, (4)

и определяющие уравнения (в данном случае закон Гука) для компонент тензора напряжений а

,0

— СуЫ£и — СуЫ (еи еы ) , (5)

где и (х) - вектор перемещения точки, определяемой радиусом-вектором х — х1е^ здесь и далее . —1,2; и° (х) - вектор перемещения точки, расположенной на поверхности ^ ; Q (х) — 0 (х) е - вектор массовых сил; р (х) — р (х) е - вектор внешней нагрузки, действующей на поверхности ; еек1 - компоненты тензора упругой деформации; е°ы -компоненты тензора начальной деформации (для термоупругого тела таковыми являются температурные деформации); СуЫ - компоненты тензора коэффициентов упругости.

При решении задачи с односторонним дискретным контактом в заданных точках контакта Ат ( т — 1,М, где М - конечное число точек контакта) дополнительно должны быть выполнены два условия контактного взаимодействия с ограничивающей контактной поверхностью ^. Для этого зададим на контактной поверхности ^ в каждой контактной

точке Аот внешнюю относительно области О нормаль пт и касательный вектор тт (рис. 1). Тогда условия контактного взаимодействия можно записать в виде следующих двух соотношений:

кинематическое условие

< < 5тп, (6)

силовое условие

атп < 0, (7)

где и™ - перемещение точки дискретного контакта Ат в направлении внешней нормали пк ; - начальное расстояние (зазор) по нормали пт между точкой дискретного контакта Аот и некоторой сходственной точкой Ап , расположенной на контактной поверхности ^, (здесь д^ > 0); а^ = ат ■ пт - проекция вектора напряжений ат на нормаль пт .

Условие (7) означает, что в точках фактического контакта растягивающие напряжения отсутствуют. Если условие (7) не выполняется (в точках нарушения контакта &т =0 ), то данная точка Аот перестает быть контактной. Необходимо отметить, что при стандартном рассмотрении контактного взаимодействия, например, двух тел, для каждого тела строится своя внешняя нормаль, поэтому нормальные проекции напряжений ат на соответствующие нормали являются отрицательными, то есть контактные нормальные напряжения являются сжимающими. При отсутствии трения касательные контактные напряжения на поверхности контакта ^ равны нулю, т.е. в данном случае в точках дискретного

контакта т = 1,М имеем: ат =0, где ат =ат • тт (ъ -касательный вектор к контактной границе ^, рис 1).

Р

Рис. 1. Схема одностороннего дискретного контакта

Многие прикладные задачи требуют учета сил трения на контактных поверхностях. В данной работе трение учитывалось в рамках закона Амонтона-Кулона [1]. Касательные

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

т т

II

если имело место скольжение контактного узла по контактной поверхности ^, здесь / - коэффициент трения (трения скольжения). А в том случае, когда наблюдается прилипание точки контакта к контактной поверхности , должно выполняться строгое неравенство

(9)

Совокупность соотношений (1) - (9) составляет математическую формулировку задачи теории упругости с односторонним дискретным контактом.

Для численного решения сформулированной выше задачи и некоторых ее вариантов был использован метод конечных элементов (МКЭ) [1, 3 - 5]. Процедура его применения описана во многих указанных руководствах. В частности, в данной работе был использован вариант МКЭ, описанный в [5].

2. Основные матричные соотношения МКЭ

Приведем основные матричные соотношения, которые были использованы при разработке программного кода для решения рассматриваемых в данной работе контактных

задач. Введем в рассмотрение следующие векторы: & = {&} =

'12

вектор напряжений;

в

'и

= {г} = <! в22 I - вектор деформации; в0 = {в} =

- вектор начальной деформации;

< 2 (* X о/ ч I < (*)

2 (*) =} 7 I - вектор массовой (объемной) силы; и (*) = \ ^ - вектор заданных

I 22 (*)1

I и0 (*

\Р1 (*)

перемещений точек поверхности ^; р (*) = < ^ - вектор заданной распределенной

[Р2 (*)

нагрузки, действующей на поверхности .

Определяющие соотношения (5) запишем в матричном виде

& = Н(в-в0), (10)

где компоненты матрицы Гука Н зависят от вида рассматриваемого напряжено-деформированного состояния.

Решение задачи (1) - (5), то есть без учета условий на границе контакта ^, эквивалентно решению соответствующей вариационной задачи, т. е. минимизации функционала полной потенциальной энергии. Выражение функционала полной потенциальной энергии

п

линеино упругого тела, размещённого в пространстве Q (х) и поверхностными р (х) имеет вид [1, 3 - 5]

и нагруженного массовыми

П = 1 — г0)т тёу uтQdv - | ит,

(11)

и

(х)!

здесь и (х) = ! 7 ^ - вектор перемещении точек тела. При минимизации функционала

1у ( х Л

(11) должны быть учтены кинематические граничные условия (2) и аналогичные кинематические условия на контактноИ поверхности ^, например, и (х= и°к (х) ,

х е ^ С , где и°к (х) - заданный вектор перемещения точки, расположенной на поверхности ^.

Минимизацию функционала (11) выполним с помощью МКЭ. В качестве конечного элемента выберем простейший 3-х узловой конечный элемент. Разобьем область С на конечные элементы. Пусть в результате разбиения получилось К - число конечных элементов и К - общее число узлов конечно-элементной модели. Область, которую занимает

конечный элемент (е), обозначим У(е) .

Глобальный вектор перемещений {и} состоит из перемещений всех узлов конечно-

элементной модели, а так как каждый узел имеет два перемещения, то размерность глобального вектора перемещений равна 2 х Ки, таким образом, имеем

и

м=

и2 и

и

(12)

здесь все перемещения узлов в направлении оси О х1 имеют нечетный индекс, а в направлении оси Ох2 - четный.

Между компонентами локального вектора узловых перемещений {ие | и компонентами глобального вектора перемещений {и} имеет место связь

{и «}-

а

(е)'

и,

и2

и

и

а

(е)

{и},

(13)

2

е

е

здесь

а

(е)

- матрица геометрических связей конечного элемента (е) [5]. Вектор перемещения произвольной точки (у(х,X) С(е)) внутри рассматриваемого элемента (е) аппроксимируем с помощью выражения

и(е) =■

,( е)

,,(е )

'М(е ) О

о

Ы(;е)

Ы(.е)

О

Ке)

О

Н(е )

о N

(е )

= [ N(е) ] {и(е)} = [ N(е) ] [а(е) ] {и}, (14)

здесь [ N(е) ] - матрица-строка функций формы конечного элемента (е).

Связь между компонентами вектора деформации {В е)} и компонентами вектора узловых перемещений {и(е )} конечного элемента (е ) определяется соотношениями Коши (4) и в матричном виде их можно записать так

в

(е) _

Г }

-(е) '11

( )

22 .(е)

[ в'-]{и« }=-!

(е)

Ъ О bJ о ък о

О о, О О; О Ск

О ъ О; ъ; Ок Ък

(15)

где £(е) - площадь конечного элемента, а учитывая соотношение (14), получим

в

(е) _

Г } =

,(е)

( )

'22

,(е) '12

> = [в(е) ]{и(е) } = [в(е) ] а(е) {и} .

(16)

Компоненты вектора напряжений внутри конечного элемента (е ) вычисляются с помощью закона Гука (10)

&

Функционал (9) достигает экстремума, если выполняется условие

аП

Н(в(е) -в0 )).

сли вы

= {О}.

а {и}

(17)

(18)

Окончательно матричное уравнение (глобальная система линейных алгебраических уравнений), эквивалентное условию (18), имеет вид

[ к]{и} = {я}, (19)

здесь [ К ] - глобальная матрица жесткости; {Я} - глобальный вектор нагрузки.

Глобальная матрица жесткости [ К ] и глобальный вектор нагрузки {Я} всей конечно-элементной модели строятся путём суммирования соответствующих локальных матриц с учётом геометрических связей, определяемых конечно-элементной моделью [5]. При решении системы (19) не обходимо учесть кинематические граничные условия (2), а также дополнительные кинематические условия на поверхности . В данной работе был использован метод сопряженных градиентов, что позволило учесть указанные условия непосредственно при реализации процедур метода сопряженных градиентов, не перестраивая для этого матрицу жесткости [ К] и вектор нагрузки {Я} .

3. Численное решение задачи теории упругости с односторонним

дискретным контактом

В качестве примера рассмотрено плоское напряженное состояние пластинки, занимающей в двухмерном пространстве М2 с декартовой системой координат Ох1х2 область С с размерами в плане 0,11м х 0,05 м . Пластинка выполнена из материала с модулем упругости Е = 200 ГПа и коэффициентом Пуассона 0,3. На нижней поверхности £ имеются выступы треугольной формы (рис. 2). Внешними вершинами выступов , т = 1,М, где М - число вершин выступов (М = 5) , пластинка опирается на абсолютно жесткое полупространство. Роль контактной поверхности выполняет поверхность полупространства. Поверхность £ , на которой заданы кинематические граничные условия (2), состоит из трех участков £п, £12 и £13: на участке £п заданы отрицательные перемещения и2 = -0,00004м в направлении оси Ох2, а на участке £12 положительные перемещения и = + 0,0001м в направлении оси Ох. Поверхность £13 зафиксирована от горизонтальных перемещений (и = 0), но вертикальные перемещения и2 точек поверхности допустимы. На поверхности заданы условия не проникновения в полупространство, то есть принято, что для вершин выступов ит = 0, т = 1,М, а ограничений на горизонтальные перемещения ит нет. Таким образом, пластинка испытывает сжатие по оси Ох и растяжение по оси Ох .

Данная задача без учета трения между выступами и контактной поверхностью является задачей теории упругости для пластинки с заданными кинематическими граничными условиями.

Рис. 2. Схема одностороннего дискретного контакта

Эта же задача была решена с учетом трения. Решение имело итерационный характер и реализовывалось следующим образом. На первом шаге выполнялось решение обычной упругой задачи без учета трения, т.е. все контактные узлы Ат имели возможность при деформировании свободно перемещаться горизонтально вдоль поверхности контакта ^. Для удобства вычислений направления осей локальных систем координат АПк^к (см. рис. 1) выбирались так, чтобы положительные направления касательной и™ и нормальной и^ компонент вектора перемещения ит в контактных узлах Ат совпадали с положительными направлениями соответствующих координатных осей Ох1 и Ох2.

После вычисления касательных перемещений ит {ит = ит), нормальных сгт

( т т \

Сп = С22)

и касательных с

( т т \

С = С12)

компонент вектора напряжений в контактных

С

С

то

узлах Ат проверялось условие прилипания (9). Если устанавливалось, что на следующей итерации в рассматриваемом контактном узле А фиксировалось значение касательной компоненты ит = 0 вектора перемещения ит, то есть было запрещено горизонтальное перемещение рассматриваемого контактного узла Аот. А в том случае, когда условие прилипания (9) не выполнялось, проводилась коррекция касательной компоненты узловой силы Рт с помощью соотношения

Р

!(к+1) = рт {к)

_ тК{к) _ тТ {к)

тах

_тК{ к)

пТ {к)

-ар:

<к)

(20)

тлт (к) ~ ~ г»т л

где р у - значение касательной компоненты узловой силы Р в контактном узле Ат (т = 1, М) после выполнения к — ой итерации; Дрт (к) - заданное приращение касательной

компоненты узловой силы Рт в контактном узле Аот, например, Дрт(к^ =а

<к)

, где

О <а< 1;

_тК(к)

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

результатам к — ой итерации,

яГ (к)

теоретическое значение модуля касательной ком-

поненты вектора напряжений, здесь

>гГ (к)

= м

а

,(к)

а

,(к)

модуль нормальной компо-

ненты вектора напряжений, вычисленной по результатам к — ой итерации.

В выражении (20) при выборе знака необходимо руководствоваться следующим:

если Рт(к)> 0 и

яГ (к)

— 8<

а

тЯ(к)

<

а

_тЯ(к)

,Г (к)

>

яГ (к)

то выбирается знак «-»; если Р^(к)> 0 и

то выбирается знак «+»; если Р^(к) <0 и

тЯ(к)

>

пГ (к)

то

выбирается знак «+»; если Р^(к) <0 и «-».

Учитывая выбор направления осей систем координат Аппк^к , в данной работе было

7-»т (к) т-»т (к) т^т (к) т^т (к)

принято, что РТ ( ) = Р ( ) и Рп ( ) = Р2 ( ) .

Так как численно добиться точного выполнения равенства (8), учитывающего скольжения, не удается, то вводится некоторый параметр 8, который позволяет организовать устойчивый численный итерационный процесс коррекции значений касательной компоненты узловой силы Рт и в том случае, когда расчетные значения касательных напряжений в рассматриваемом контактном узле приближаются (по модулю) к теоретическим значениям атГ снизу. В процессе выполнения итераций модуль разности значе-

яГ (к)

— 8<

а

тЯ(к)

<

а

чГ (к)

то выбирается знак

ний касательных напряжений Да^ =

а

1Я(к)

а

чГ (к)

стремится к нулю.

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

Процесс вычислений заканчивается после оценки сходимости результатов итерации, например, для этого может быть использовано соотношение

<е, (21)

т (к+1) т (к)

и ( ) — и ( )

т (к+1) т а л -, *

где и у - вектор перемещения и контактного узла Аот, т = 1,М, после выполнения (к +1) — ой итерации; еи - требуемая точность по перемещениям контактных узлов.

На рис. 3 показаны конечно-элементная модель пластинки (рис. 3, а) и в увеличенном масштабе ее деформированное представление без учета (рис. 3, б) и с учетом (рис. 3, в) трения.

На рис. 4 - 8 представлены поля компонент компонент тензора напряжений и вектора перемещений, вычисленные без учета и с учетом трения.

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

Заметное количественное различие можно обнаружить при сравнении результатов численных исследований (см. табл. 1) для контактных узлов Ат, т = 1,М, здесь М = 5 (рис. 3, а), полученных при решении рассматриваемой задачи без учета трения (числитель) и с учетом трения (знаменатель).

в

Рис. 3. Конечно-элементная модель пластинки: а - до деформирования; б - после деформирования без учета

трения; в - после деформирования с учетом трения

б

Рис. 4. Поле компоненты тензора напряжений ои (МПа): а - без учета трения; б - с учетом трения

б

Рис. 5. Поле компоненты тензора напряжений а22 (МПа): а - без учета трения; б - с учетом трения

б

Рис. 6. Поле компоненты тензора напряжений о12 (МПа):а - без учета трения; б - с учетом трения

б

Рис. 7. Поле компоненты вектора перемещений Щ (м): а - без учета трения; б - с учетом трения

Рис. 8. Поле компоненты вектора перемещений и2 (м): а - без учета трения; б - с учетом трения

Таблица 1. Результаты сравнительных численных исследований

Номер контактного узла и х 10+4, м р, МН р, МН с, МПа а22, МПа аи, МПа

1 0,15579 0,0 1,657 41,985 -331,409 0,0

0,01727 -0,480 1,601 14,580 -320,238 96,071

2 0,32427 0,0 1,534 51,107 -306,823 0,0

0,17218 -0,463 1,545 36,611 -309,034 92,710

3 0,50213 0,0 1,520 52,961 -304,018 0,0

0,34588 -0,458 1,528 47,092 -305,603 91,681

4 0,68031 0,0 1,535 51,310 -307,086 0,0

0,52793 -0,462 1,538 53,661 -307,704 92,311

5 0,85141 0,0 1,657 44,012 -331,341 0,0

0,70284 -0,521 1,736 59,491 -347,379 104,213

На рис. 9 и 10 представлены зависимости, характеризующие процесс сходимости при решении задачи теории упругости с дискретным односторонним контактом с учетом трения. Цифрами обозначены номера контактных узлов (см. рис. 3, а).

0.5 -I -К5 -2 ■2.1

■г!

\

л 4 5

// \ 2

/

\ . 5

. ч

1 ■ / 1 /

о и» 20 зп «I 51> л1 (| | о ли за р ¡о IV

а б

Рис. 9. Изменения значений компонент вектора узловых сил Рт контактных узлов Ат в зависимости от числа итераций N : а - касательные компоненты; б - нормальные компоненты

Рис. 10. Изменения значений компонент вектора напряжений &т (а - в) и касательной компоненты ит = ит вектора перемещения ит (г) контактных узлов Ат в зависимости от числа итераций N

Из графиков на рис. 9-10 видно, что использование соотношения (20) для коррекции касательных сил, учитывающих эффект трения в контактных узлах конечно -элементной модели, позволяет получить устойчивое сходящееся решение задачи теории упругости с односторонним контактом.

Заключение

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

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

Работа выполнена при поддержке гранта государственной поддержки ведущих научных школ НШ-1432.2014.8.

Список литературы

1. Зарубин В.С., Станкевич И.В. Расчет теплонапряженных конструкций. - М.: Машиностроение, 2005. - 352 с.

2. Галин Л.А. Развитие теории контактных задач М.: Наука, 1976. -494.

3. Зенкевич О. Метод конечных элементов в технике. - М.: Мир, 1975. - 540 с.

4. Сегерлинд Л. Применение метода конечных элементов. - М.: Мир, 1979. - 392 с.

5. Котович А.В., Станкевич И.В. Решение задач теории упругости методом конечных элементов. - М.: Изд-во МГТУ им. Н.Э. Баумана, 2012. - 106 с.

6. Галанин М.П., Крупкин А.В., Кузнецов В.И. Лукин В.В., Новиков В.В., Родин А.С., Станкевич И.В. Математическое моделирование термоупругогопластического контактного взаимодействия системы тел // Mathematica Montisnigri, 2014. - Т. 30 . - С. 99

- 114.

7. Станкевич И.В., Яковлев М.Е., Си Ту Хтет. Разработка алгоритма контактного взаимодействия на основе альтернирующего метода Шварца // Вестник МГТУ им. Н.Э. Баумана. Серия "Естественные науки". Спецвыпуск: Прикладная математика, 2011 г.

- С. 134 - 141.

8. Лукашевич А.А., Розин Л.А. О решении контактных задач строительной механики с односторонними связями и трением методом пошагового анализа // Инженерно-строительный журнал. 2013. №1. С. 75-81.

9. Silveira R.A., Gon9alves P.B. Analysis of slender structural elements under unilateral contact constraints. Structural Engineering and Mechanics, 2001, v. 12, no. 1, pp.1-16.

10. Elabbasi N., Hong J.W., Bathe K.J.R. On the reliable solution of contact problems in engineering design. Int. J. of Mechanics and Materials in Design, 2004, no 1, pp. 3-16.

11. Fernandez J.R., Sofonea M. Variational and numerical analysis of the signorini's contact problem in viscoplasticity with damage. Journal of Applied Mathematics, 2003, no 2, pp. 87-114.

Mathematics I Mathematical Modelling

Mathematics and Mathematical Madelling of the Bauman MSTU, 2015, no. 04, pp. 93-110.

DOI: 10.7463/mathm.0415.0801840

Received: Revised:

08.07.2015 22.07.2015

© Bauman Moscow State Technical Unversity

Mathematical Modeling of Contact Problems of Elasticity Theory with Unilateral Discrete Contact

I.V. Stankevich1' ''aplmexjgyandexiu

:Bauman Moscow State Technical University, Moscow, Russia

Keywords: discrete contact interaction; a contact problem of elasticity theory; finite element method

Development and operation of modern machinery and latest technology require reliable estimates of the strength characteristics of the critical elements of structures and technological equipment under the impact of high-intensity thermomechanical loading, accompanied, as a rule, by complex contact interaction. Mathematical modeling of stress-strain state of such parts and components in the contact area, based on adequate mathematical models, modern numerical methods and efficient algorithms that implement the direct determination of displacement fields, strains and stresses, is the main tool that allows fast acquisition of data required for the calculations of strength and durability. The paper considers an algorithm for constructing the numerical solution of the contact problem of elasticity theory in relation to the body, which has an obvious one-sided discrete contact interaction with an elastic half-space. The proposed algorithm is specially designed to have a correction of the tangential forces at discrete contact points, allowing us to achieve sufficiently accurate implementation of the adopted law of friction. The algorithm is embedded in a general finite element technology, with which the application code is generated. Numerical study of discrete unilateral contact interaction of an elastic plate and a rigid half-space showed a high efficiency of the developed algorithm and the application code that implements it.

References

1. Zarubin V. S., Stankevich I. V. Rachet teplonapryazhenykh konstruktsiy [Calculation of heat-stressed structures]. Moscow, Mashinostroenie Publ., 2005. 352 p.

2. Galin, L. A. Razvitie teorii kontaktnykh zadach [Development of the theory of contact problems]. Moscow, Nauka Publ., 1976. 494 p.

3. Zenkevich O. Metod konechnix elementov v tehnike [The Finite Element Method in Engineering Science]. Moscow, Mir Publ., 1975. 540 p.

4. Segerlind L. Primenenie metoda konechnix elementov[Applied Finite Element Analysis]. Moscow, Mir Publ., 1979. 392 p.

5. Kotovich, A. V., Stankevich I. V. Reshenie zadach teorii uprugosti metodom konechnix elementov [Solution of elasticity problems by the finite element method]. Moscow, Bauman MSTU Publ., 2012. 106 p.

6. Galanin M. P., Krupkin A. V., Kuznetsov V. I. Lukin V. V., Novikov V. V., Rodin A. S., Stankevich I. V. Mathematical modeling thermoplastische contact interaction of a system of bodies. Mathematica Montisnigri, 2014, vol. 30, pp. 99-114.

7. Stankevich I. V., Yakovlev, M. E., Si Tu Htet. Development of contact interaction algorithm based on the alternating Schwarz method. Vestnik Bauman MSTU. Series "Natural Sciences". Special issue: Applied mathematics, 2011, pp. 134-141(in Russian).

8. Lukashevich A.A., Rozin L.A. On the decision of contact problems of structural mechanics with unilateral constraints and friction by step-by-step analysis. Magazine of Civil Engineering, 2013, no. 1, pp. 75-81 (in Russian).

9. Silveira R.A., Gon9alves P.B. Analysis of slender structural elements under unilateral contact constraints. Structural Engineering and Mechanics, 2001, vol. 12, no. 1, pp.1-16.

10. Elabbasi N., Hong J.W., Bathe K.J.R. On the reliable solution of contact problems in engineering design. Int. J. of Mechanics and Materials in Design, 2004, no 1, pp. 3-16.

11. Fernandez J.R., Sofonea M. Variational and numerical analysis of the signorini's contact problem in viscoplasticity with damage. Journal of Applied Mathematics, 2003, no 2, pp. 87-114.