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

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

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

Сетевое научное издание http://mathmelpub.ru

Ссылка на статью:

// Математика и математическое моделирование. 2017. № 05. С. 40-53

Б01: 10.24108/шаШш.0517.0000078

Представлена в редакцию: 30.09.2017

© НП «НЕИКОН»

УДК 539.37

Численное решение смешанных задач теории упругости с односторонними связями

Станкевич И.В.1' 'ар1тех@уаги1е;ии

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

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

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

Введение

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

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

ния таких задач с помощью конечно-элементной технологии обычно реализуются в рамках метода перемещений [2 - 6]. При этом из-за низкого уровня аппроксимации напряжений и деформации по сравнению с перемещениями возникают сложности при вычислении реакций связей в зонах контакта. В то же время контактные напряжения являются основными искомыми функциями в задачах теории упругости с односторонними связями, следовательно, должны быть определены с достаточно высокой степенью точности. Кроме того, при рассмотрении подобного типа задач ситуация с получением решения значительно усложняется, если априори неизвестны пределы границы контактного взаимодействия.

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

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

1. Математическая постановка задачи

Математическая формулировка квазистатической несвязанной задачи механики деформируемого твердого тела в рассматриваемой постановке включает следующие уравнения:

уравнения равновесия

(и,Т) + а (х) = 0, Xев, (1)

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

и (х)| = и0 (х), х е ^ сд в, (2)

(и,Т)п\ = р (х), х е ^ сдв, (3)

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

еа(х>1 ) =1 (и,; +и1 Ч), х ев, (4)

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

= Сцые« = Суы (е« ), (5)

где в - ограниченная область с границей дв; и (х) - вектор перемещения точки, определяемой радиус-вектором х = xiei, в общем случае / = 1,3; и0 (х) - вектор перемещения точки, расположенной на поверхности ; Q (х) = ^ (х) ei - вектор массовых сил, здесь

х еО; р (х) = р (х) в{ - вектор внешней нагрузки, действующей на поверхности , здесь х е сб О; геы - компоненты тензора упругой деформации; г°ы - компоненты тензора начальной деформации (для термоупругого тела таковыми являются температурные деформации); СуЫ - компоненты тензора коэффициентов упругости.

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

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

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

ип <Ъп, (6)

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

а <0, (7)

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

При отсутствии трения касательные напряжения на контактной поверхности ^ равны нулю, т.е. имеем: а = 0, где ах=а-т (т - касательный вектор к контактной поверхности ^, рис. 1). Многие прикладные задачи требуют учета сил трения на контактных поверхностях [4, 7]. В данной работе трение учитывалось в рамках закона Амонтона - Кулона [2]. Модуль касательного напряжения а в точках контактной поверхности вычислялся по формуле

Ы = ц|Ъп|, (8)

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

Н<цК|, (9)

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

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

Для численного решения задачи (1) - (9) применим смешанную формулировку метода конечных элементов, основанную на условии стационарности функционала Рейсснера [8].

Функционал Рейсснера запишем в следующем виде ф (а и) = |сттAudV -1 |атDстdV -1 и^ОёУ -1 (и - И )Т АаdS -1 ит, (10)

О О О ¿1

здесь А - матрица операций дифференцирования; Б - матрица податливости (обратная к матрице Гука Н = Б-1); А - матрица направляющих косинусов внешней нормали п к поверхности ^ .

Пусть конечно-элементная модель области О состоит из К однотипных конечных элементов (е), каждый из которых имеет М узлов и соответственно столько же функций

формы л), Р = 1, М , зависящих от локальных координат £ и л. Общее число уз-

лов конечно-элементной модели - N .

т

и

N

. С помощью

Введем в рассмотрение две матрицы функций формы

этих матриц функций формы аппроксимируем векторы напряжений а и перемещений и в пределах конечного элемента (е) . Имеем

(е) _

а

(е) И ' =

а ^х3М I )э

ы:

2х2М

и

2М

N

ы:

М"

3х3М

2х2М

а,

а.

М"

3М х3Ы

2М х2 N

{и}2N ■

(11) (12)

здесь {а(е)} - локальный вектор напряжений, состоящий из компонент тензора напряжений, относящихся только к М узлам фиксированного конечного элемента (е); {а} - глобальный вектор напряжений, состоящий из компонент тензора напряжений, относящихся

ко всем N узлам конечно-элементной модели;

а\

матрица геометрических связей

конечного элемента (е), используемая для связи компонент локального вектора напряжений {а(е)}и компонент глобального вектора напряжений {а} [9, 10]; {м(е)}- локальный

вектор перемещений, состоящий из компонент вектора перемещений, относящихся только к М узлам фиксированного конечного элемента (е); {и} - глобальный вектор перемещений, состоящий из компонент вектора перемещений, относящихся ко всем N узлам ко-

нечно-элементной модели;

а

матрица геометрических связей конечного элемента (е), используемая для связи компонент локального вектора перемещений {М е)} и компонент глобального вектора перемещений {и} [9, 10].

Окончательно функционал Рейсснера (10) в матричном виде для конечно-элементной модели можно записать так

е=1

К'

а),

{и}-1 {;}т±[4

е=1

I[

1 О

а

■{и}ТЕ[а,(*)]' {1«}-{и}ТЦ>18)]' {I<«)},

е=1 8=1

(13)

где

гМ"

Я Т [А]

N

М"

ёУ:

К')'

{131=1 [ N

У (')

N

М"

= ! [N1'']' [О]

У (')

ёУ {е(е)},

N

ёУ;

здесь

{о"1}

локальный вектор массовых сил, компоненты которого отнесены к узлам

элемента (е); {48)}= | Г

е)

N

(8)'

N

(8)"

{р1-8; Ь - число граней, в данном случае - од-

а..

- мат-

номерных конечных элементов, которые аппроксимируют поверхность $2; рица геометрических связей одномерного конечного элемента - грани (8), используемая для связи компонент локального вектора перемещений {^8и компонент глобального вектора перемещений {и}; {р(8 - локальный вектор внешней нагрузки, компоненты которого отнесены к узлам грани (8).

Так как заданные перемещения и0 на поверхности 5 можно учесть при решении системы линейных алгебраических уравнений (СЛАУ), то нет необходимости формиро-

вать матрицы, связанные с вычислением слагаемого

\Т

ис

К и - и0) А .

5,

Введем обозначения

К г

[ *п Ы [а

М"

е=1 Ке

[ К12 Ы [а

(е)"

е=1

К'

-Т " Г(е)] 1 2 Г (е)! «а

-Т Г К(')] К1 Г (е)! _ «и _

(14)

№ЫК]' {!:(е)М[а!8]]' {48)}.

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

равно . Для их вычисления необходимо построить СЛАУ, используя свойство стационарности функционала Рейсснера. Для построения СЛАУ надо продифференцировать

2

Т

Т

К

функционал ф({а},{и}) по компонентам глобальных векторов напряжений {а} и перемещений {и} и производные приравнять нулю, в результате получим

>}} {КЛ .{и0*2 }) '

K ] [K12 ]

[ K21 ] [ K22 ]

(15)

где [ K12 ] = [ K21 ]T; [ K22 ] = [0] - нулевая матрица; {^} = {0} - нулевой вектор.

3. Схема решения седловой задачи

Применим модифицированный метод SSOR (метод MSSOR) для численного решения системы (15), в основе которого лежит метод последовательной верхней релаксации (SOR - Successive Over Relaxation) [11]. Имеем

1 [ Kii ](Mk+^-Mk ) + [ K11 ]{a}k +[ K12 ]{U }k ={0};

B ]({U }k+1 -{U}k)+[ K12 ]T {a}k+12 = ft}; (16)

1 [ K11 ]({a}k+1 -{a}k+^ ) + [ K11 ]{a}k+^ +[ K12 ]{U }k+1 ={0},

здесь [B] - предобуславливатель; k - номер итерации; a и т - итерационные параметры. Применение схемы (16) требует в рамках одной итерации решения трех СЛАУ относительно глобальных векторов приращений напряжений {До} и перемещений {AU}

[ K„ ]{Ao}k+* =-т([ Ku ]{o}k +[ K12 ]{U }k);

[ b]{aU }k+1 = 2T([ K12 ]T {o}k+X-{R});

(17)

[ Kn ]{Ao}k1 =-т( [ Kn ]{o}k+12 +[ Ku ]{U }

k +1

\k+1

где

{а}к ^ = {а}к + {Аа}"/2; {а}k+1 = {а}"12 + {Aа}k+1 и ^}k+1 = {^ + {А^

Оценка численных значений итерационных параметров и сходимость вычислительной схемы (17) рассмотрены в работе [11]. Каждое СЛАУ, входящее в (17) можно решить с помощью метода сопряженных градиентов [12], некоторые особенности реализации которого в случае разреженных матриц изложены в [13].

ík+Уг. 1„ \k+>2

ik+1

4. Результаты численных исследований

В качестве примера рассмотрим плоское напряженное состояние пластинки, занимающей в двухмерном пространстве М2 с декартовой системой координат Oх:х2 область

О (см. рис. 1). Пластинка выполнена из материала с модулем упругости Е = 200 ГПа и коэффициентом Пуассона V = 0,3. Поверхность 5 зафиксирована от горизонтальных перемещений (и = 0), но вертикальные перемещения и2 точек поверхности допустимы. На поверхности £2 задана равномерно распределенная нагрузка с интенсивностью р . Снизу

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

Для решения были выбраны 8-ми узловые конечные элементы, допускающие квадратичную аппроксимацию как компонент вектора перемещения, так и тензора напряжений.

На рис. 1-5 представлены расчетная схема и поля компонент тензора напряжений ап, о22 и вектора перемещений их, и2, вычисленные с учетом трения и заданной равномерно распределенной нагрузкой с интенсивностью р2 = -10 МПа.

Как и следовало ожидать, поля компонент тензора напряжений и вектора перемещения имеют выраженный неоднородный характер. Наиболее нагруженной зоной является зона контакта пластинки и полупространства, а зоны периферийных областей менее напряженными (см. рис. 2 и рис. 3). Для области пластики вблизи поверхности приложения внешней нагрузки £2 в направлении оси Ох1 характерны растягивающие напряжения, в отличие от области в районе поверхности ^, в которой в указанном направлении получены сжимающие напряжения (см. рис. 2 и рис. 4).

Рис. 1. Расчетная схема

В направлении оси Oх2 зона сжатия локализована вблизи нижнего углового узла

(узел 1, см. рис. 1), а большая области анализа находится в состоянии растяжения (см. рис. 3 и 5).

Кроме того, особенностью данной задачи является необходимость определения в процессе решения тех узлов конечно-элементной модели, которые непосредственно находятся в контакте с жестким основанием. Этим можно объяснить заметное различие в характере напряженно-деформированного состояния, которое получено при решении данной задачи и задач, например, рассмотренных в работах [14, 15].

Рис. 2. Поле компоненты тензора напряжений C j (МПа)

Рис. 3. Поле компоненты тензора напряжений G22 (МПа)

Решение было получено в результате выполнения 6 итераций по схеме (17). Все исследования проводились на компьютере с процессором Intel Core 2 Duo T9550. Так как сетка содержала 121 узел (рис. 1), то число неизвестных было небольшим - около 600, поэтому время решения на данном компьютере было весьма незначительным.

Рис. 4. Поле компоненты вектора перемещений Щ (м)

Рис. 5. Поле компоненты вектора перемещений и2 (м)

Аналогичные работы проводили многие исследователи. Например, в работах [3, 4] сходимость была достигнута за 10 итераций. В работе [5] в результате применения метода множителей Лагранжа в рамках конечно-элементной технологии и реализации шаговой нагрузки сходимость была получена за 40 шагов. В работе [6] был использован альтернативный метод Шварца, что позволило получить приемлемое решение за 12-15 итераций. За такое же число итераций получены решения задач теории упругости с односторонними связями при рассмотрении соответственно дискретного и непрерывного контакта в работах [14, 15]. Приведенные данные позволяют надеяться, что предложенный алгоритм имеет перспективу практического применения.

Заключение

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

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

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

1. Чирков А.Ю. Применение смешанных вариационных формулировок метода конечных элементов к решению задач о собственных колебаниях упругих тел // Проблемы прочности. 2008. № 2. С. 121 - 140.

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

3. Можаровский Н.С., Качаловская Н.Е. Приложение методов теории пластичности и ползучести к решению инженерных задач машиностроения: учебник: В 2 ч. Ч. 2: Методы и алгоритмы решения краевых задач. Киев: Выща школа, 1991. 288 с.

4. Задачи контактного взаимодействия элементов конструкций / Подгорный А.Н., Гон-таровский П.П., Киркач Б.Н. и др.; отв. ред. Рвачев В.Л. Киев: Наукова думка, 1989. 232 с.

5. Zienkiewicz O.C., Taylor R.L. The finite element method for solid and structural mechanics. 6th ed. Amst.; Boston: Elsevier; Butterworth-Heinemann, 2005. 631 p.

6. Яковлев М.Е. Математическое моделирование контактного взаимодействия термовяз-коупругопластических сред: дис. ... канд. техн. наук. М.: МГТУ им. Н.Э. Баумана, 2014. 131 с.

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

С. 134-141.

8. Розин Л.А. Вариационные постановки задач для упругих систем. Л.: Изд-во ЛГУ, 1978. 223 с.

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

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

11. Быченков Ю.В., Чижонков Е.В. Итерационные методы решения седловых задач. М.: БИНОМ. Лаборатория знаний, 2010. 349 с.

12. Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. М.: Наука, 1978. 591 с.

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

14. Станкевич И.В. Математическое моделирование задач теории упругости с односторонним дискретным контактом // Математика и математическое моделирование. 2015. № 4. С. 93-110. DOI: 10.7463/mathm.0415.0801840

15. Станкевич И.В. Математическое моделирование контактных задач теории упругости с непрерывным односторонним контактом // Математика и математическое моделирование. 2015. № 5. С. 83-96. DOI: 10.7463/mathm.0515.0812348

Mathematics & Mathematical Modelling

Electronic journal http://mathmelpub.ru

Mathematics and Mathematical Modeling, 2017, no. 05, pp. 40-53.

DOI: 10.24108/mathm.0517.0000078

Received: 30.09.2017

© NP "NEICON"

Numerical Solution of Mixed Problems of the Theory of Elasticity with One-Sided Constraints

I.V. Stankevich1'* ''aplmexgjyandexju

1Bauman Moscow State Technical University, Moscow, Russia

Keywords: mixed problem of the theory of elasticity; the Reissner functional; finite element method

The paper deals with the application features of the finite element technologies to solve the problems of elasticity with one-sided constraints. On the one hand, the area of this study is determined by the fact that many critical parts and assemblies of mechanical and power engineering constructions have a significant contact within some given surface. To assess the strength and the life of these parts and assemblies, reliable stress-strain state data are demandable. Data on the stress-strain state can be obtained using the contemporary mathematical modeling means, e.g., finite element technology.

To solve the problems of the theory of elasticity with one-sided constraints, a method of finite elements in a traditional classical form can be used, but it is necessary to consider some of its shortcomings. The most significant one is an approximation of the tensile stress and strain, as well as a considerably lower order of convergence of the approximation for stresses and strains as compared to displacements. Improving the accuracy through increasing a density of the finite element models and/or the transition to more complex approximations is not always optimal, because increasing a dimension of the discrete problem leads to a significant computational cost and demand for expensive computing resources.

One of the alternatives in numerical analysis of contact problems of the elasticity theory is to use the mixed variational formulations of the finite element method in which stresses and/or strains appear in the resolving equations along with displacements as equal unknown. A major positive factor when using the mixed formulations of the finite element method is reduction of the approximation error of stress and strain, which leads to a more accurate assessment of the stress-strain state in comparison with the classical approach of the finite element method in the form of the method of displacements.

Besides, mixed schemes of the finite element method enable us to ensure continuous approximation of not only displacements, but also stresses and strains. Mixed schemes to solve the boundary value problems lead to the saddle-point problems. Their solutions use various iterative

techniques. One of the most effective techniques is a modified SSOR (MSSOR) method, based on the SOR (Successive Over Relaxation) one.

The paper considers one of the options of the finite element method in the framework of mixed scheme that uses a Reissner functional. The procedures of the algorithm proposed in the paper are used to solve the problem of contact interaction when an elastic body of the finite dimensions, being under a load of the external forces, relies on the absolutely rigid half-space. The contact occurs with the distinguished contact surface, which in the general case can change its size during thermo-mechanical loading. The algorithm is implemented as an application software complex. The numerical study of the one-sided contact interaction between the elastic plate and the perfectly rigid half-space has shown a fairly high efficiency of the developed algorithm and the code that implements it.

References

1. Chirkov A.Yu. Application of mixed variational formulations based on the finite element method to the solution of problems on natural vibrations of elastic bodies. Strength of Materials, 2008, vol. 40, no. 2, pp. 253-268. DOI: 10.1007/s11223-008-9005-3

2. 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. Inzhenerno-stroitel'nyj zhurnal [Magazine of Civil Engineering], 2013, no. 1(36), pp. 75-81.

DOI: 10.5862/MCE.36.9 (in Russian)

3. Mozharovskij N.S., Kachalovskaia N.E. Prilozhenie metodov teorii plastichnosti i polzuchesti k resheniyu inzhenernykh zadach mashinostroeniia [The application of methods of the theory of plasticity and creep to the solution of engineering tasks engineering]: a textbook: in two parts. Pt. 2: Metody i algoritmy resheniia kraevykh zadach [Methods and algorithms for solving boundary value problems]. Kiev: Vyshcha shkola Publ., 1991. 288 p. (in Russian).

4. Zadachi kontaktnogo vzaimodejstviia elementov konstruktsij [Problems of contact interaction of structural elements] / Podgornyj A.N., Gontarovskij P.P., Kirkach B.N. a.o.; ed. by V.L. Rvachev. Kiev: Naukova Dumka Publ., 1989. 232 p. (in Russian).

5. Zienkiewicz O.C., Taylor R.L. The finite element method for solid and structural mechanics. 6th ed. Amst.; Boston: Elsevier; Butterworth-Heinemann, 2005. 631 p.

6. Yakovlev M.E. Matematicheskoe modelirovanie kontaktnogo vzaimodejstviia termoviazkouprugoplasticheskikh sred [Mathematical simulation of contact interaction thermoviscoplasticity environments. Cand. diss]. Moscow: Bauman MSTU Publ., 2014. 131 p. (in Russian).

7. Stankevich I.V., Yakovlev M.E., Si Tu Khtet. Development of contact interaction algorithm on the basis of Schwarz alternating method. Vestnik MGTU im. N.E. Baumana. Ser. Estestvennye nauki [Herald of the Bauman MSTU. Natural Sciences], 2011, spec. iss. "Applied Mathematics", pp. 134-141 (in Russian).

8. Rozin L.A. Variatsionnye postanovki zadach dlia uprugikh sistem [Variational formulation of problems for elastic systems]. Leningrad: Leningrad Univ. Publ., 1978. 223 p. (in Russian).

9. Zarubin V.S., Stankevich I.V. Raschet teplonapriazhennykh konstruktsij [Calculation of heat-stressed structures]. Moscow: Mashinostroenie Publ., 2005. 351 p. (in Russian).

10. Kotovich A.V., Stankevich I.V. Reshenie zadach teorii uprugosti metodom konechnykh elementov [The solution of elasticity problems by the finite element method]. Moscow: Bauman MSTU Publ., 2012. 106 p. (in Russian).

11. Bychenkov Yu.V., Chizhonkov E.V. Iteratsionnye metody resheniia sedlovykh zadach [Iterative methods for solving saddle task]. Moscow: BINOM. Laboratoriia znanij Publ., 2010. 349 p. (in Russian).

12. Samarskij A.A., Nikolaev E.S. Metody recheniia setochnykh uravnenij [Methods of solving the grid equations]. Moscow: Nauka Publ., 1978. 591 p. (in Russian).

13. Kotovich A.V., Stankevich I.V. Reshenie zadach teploprovodnosti metodom konechnykh elementov [The solution of heat conduction problems by the finite element method]. Moscow: Bauman MSTU Publ., 2010. 84 p. (in Russian).

14. Stankevich I.V. Mathematical modeling of contact problems of elasticity theory with unilateral discrete contact. Matematika i matematicheskoe modelirovanie [Mathematics and Mathematical Modelling], 2015, no. 4, pp. 93-110. DOI: 10.7463/mathm.0415.0801840 (in Russian)

15. Stankevich I.V. Mathematical modeling of contact problems of elasticity theory with continuous unilateral contact. Matematika i matematicheskoe modelirovanie [Mathematics and Mathematical Modelling], 2015, no. 5, pp. 83-96. DOI: 10.7463/mathm.0515.0812348 (in Russian)