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

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

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

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

Б01: 10.7463/шаШш.0515.0812348

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

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

УДК 539.37

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

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

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

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

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

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

Введение

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

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

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

контактом

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

ар , ( и ) + £г ( х) = 0, х еС , (1)

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

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

аЛ (И) = Р (*) ' х Е ^ д О (3)

2

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

еЛх ) = (и * +UJ *)' х еО , (4)

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

= СуИ£И = Сук1 {£Ы ) , (5)

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

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

ип < , (6)

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

< 0, (7)

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

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

ат =Иап , (8)

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

<Ма„

, (9)

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

2. Алгоритм численного решения контактной задачи теории упругости с односторонним непрерывным контактом

Процесс решения имеет итерационный характер и реализуются следующим образом. Исходная задач рассматривается как задача теории упругости, при этом на контактной поверхности 8к задаются кинематические (по нормали п ) и силовые (по касательной т ) условия. Если силовые условия отсутствуют, то данная задача становится стандартной задачей теории упругости (см. рис. 1, а). Силовые условия необходимо учитывать в том случае, когда задача решается с учетом трения (см. рис. 1, б). Так как рассматривается односторонний контакт, то кинематические условия заранее известны, например, и\ = 0 (см.

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

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

гии линейно упругого тела, размещённого в пространстве М2 и нагруженного массовыми Q (х) и поверхностными р (х), получается матричное уравнение (глобальная система линейных алгебраических уравнений), имеющее вид [4],

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

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

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

б

Рис. 1. Схема одностороннего непрерывного контактного взаимодействия: а - без учета трения; б - с учетом

трения

а

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

5, 6, 7]. Алгоритм коррекции касательных сил реализуется следующим образом. На первом шаге выполняется решение обычной упругой задачи без учета трения, то есть все узлы конечных элементов, расположенные на контактной поверхности ^ имели возможность при деформировании свободно перемещаться вдоль ограничивающей поверхности

т—г т т ».»»..

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

m »-» m »-» m

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

рт

Лк+1) = рт (к)

_тй(к) _тТ (к)

тах

_тЯ(к)

_тТ (к)

-АР'

(к)

(11)

Т>т (к) ^ ^ пт л

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

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

рт

г(к)

, где

0 <а< 1;

_тК(к)

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

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

а

чТ (к)

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

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

а

чТ (к)

= \

а

г (к)

а

г(к)

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

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

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

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

_тК(к)

>

пТ (к)

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

пТ (к)

-8<

а

тК(к)

<

пТ(к)

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

_шЯ(к)

>

пТ (к)

то

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

пТ (к)

-8<

а

тК(к)

<

пТ (к)

, то выбирается

знак «-».

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

пт

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

« "Л _т _шК( к) _тТ (к)

ний касательных напряжений Д< = аТ ' - аТ ' стремится к нулю.

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

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

8{к+1) =

тах

т=1,М

т (к+1) т (к)

и v ; - и у '

< £„

(12)

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

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

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

О с размерами в плане 0,11м х 0,05 м (см. рис. 1). Пластинка выполнена из материала с модулем упругости Е = 200 ГПа и коэффициентом Пуассона 0,3. Поверхность ^ зафиксирована от горизонтальных перемещений (щ = 0), но вертикальные перемещения щ точек поверхности допустимы. На поверхности задана равномерно распределенная нагрузка с интенсивностью р = 350МПа . На контактной поверхности ^ заданы условия не проникновения в полупространство, т.е. принято, что для всех узлов Аот ^ = 0, а специальных кинематических ограничений на касательные перемещения щт нет. Таким обра-

зом, пластинка в результате приложения внешней нагрузки р на поверхности $2 испытывает сжатие по оси Ох2 и растяжение по оси Охх.

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

Эта же задача была решена с учетом трения (^ = 0,05) . Решение имело итерационный характер и реализовывалось так, что в процессе итераций уточнялись значения касательных напряжений ак, возникающих на контактной поверхности 8к, но при этом сохранялись условия не проникновения в полупространство, т.е. и™ = 0 (рис. 1, б).

Обе задачи были решены с использованием конечно-элементной технологии [4]. На рис. 2 показана конечно-элементная модель пластинки, состоящая из 480 простейших 3-х узловых конечных элементов с общим числом узлов 275. Число контактных узлов в конечно-элементной модели было равно 25.

Рис. 2. Конечно-элементная модель пластинки

При решении данной задачи без учета трения вертикальное перемещение узлов, расположенных на поверхности , было максимальным по абсолютной величине и составило значение щ = —0,875 х10—4 м. Максимальные горизонтальные перемещения имели узлы на поверхности £3: щ = 0,630*10 4м. Напряженное состояние пластинки было полностью однородным: ап =^12 = 0,0 МПа и а22 = —350,0 МПа.

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

^ Л1 у" гЬ ч-^-"--г-С- щ э ж; ад ~>к ¿Л

.я* .А"^ -ос й

"'1

'""Г ^"",1 ^

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

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

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

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

— <И. 72

Рис. 6. Поле компоненты тензора напряжений <ги с учетом трения (МПа)

б

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

=

б

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

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

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

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

Номер контактного узла и , м Р, МН Р, МН аи, МПа а22, МПа аи, МПа

1 0.10500е-04 0,0 1.750 0,0 -350,00 0,0

0.3021167е-05 -0.10 1.71 -47.69 -342.56 17.13

2 0.23625е-04 0,0 1.750 0,0 -350,000 0,0

0.1111654е-04 -0.09 1.73 -31.55 -346.05 17.30

3 0.36750е-04 0,0 1.750 0,0 -350,00 0,0

0.2081970Е-04 -0.09 1.74 -22.21 -348.40 17.42

4 0.49875е-04 0,0 1.75 0,0 -350,00 0,0

0.3159375е-04 -0.09 1.77 -16.22 -354.18 17.71

5 0.63000Е-04 0,0 0.88 0,0 -350,00 0,0

0.4282310е-04 -0.15 1.06 -25.94 -394.39 19.72

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

£„41 [

V 1

/у

ч< Г4*

1ч*

■

а б

Рис. 9. Изменения значений компонент вектора узловых сил Р контактных узлов в зависимости от числа

итераций N

п .. Л1Г 1.1

-

¡Я— 4 >

Г

1

й, .VI]Ь

\ /'

д е

Рис. 10. Изменения значений компонент тензора напряжений & и компоненты И] вектора перемещения и

контактных узлов в зависимости от числа итераций N

в

г

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

Заключение

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

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

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

1. Станкевич И. В. Математическое моделирование задач теории упругости с односторонним дискретным контактом // Математика и математическое моделирование. 2015. № 04. http://mathmjournal.ru/doc/801840.html. DOI: 10.7463/mathm.0415.0801840.

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

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

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

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

- 114.

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

- С. 134 - 141.

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

8. 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.

9. 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. 10. 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. 05, pp. 83-96.

DOI: 10.7463/mathm.0515.0812348

Received: Revised:

10.09.2015 24.09.2015

© Bauman Moscow State Technical Unversity

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

I.V. Stankevich1' ''aplmexigyandexju

1Bauman Moscow State Technical University, Moscow, Russia

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

The work [1] presents the formulation and numerical solution of the problem concerning the unilateral discrete contact interaction of an elastic body and a rigid half-space. However, many parts and components of engineering structures have a pronounced continuous contact within a given surface [2, 3]. In this paper we consider a special case of this option of contact interaction when, the elastic body of finite size, subjected to external forces, is based on a rigid half-space. Contact occurs through a dedicated contact surface, which in general can change their sizes.

Developed to solve this problem, a numerical algorithm is a further adaptation and development of the approaches described in [1]. The paper shows results of solving the model problem of the elasticity theory with and without taking friction into account. In the latter case, were additionally obtained numerical data characterizing the convergence of the solution

References

1. Stankevich I. V., Mathematical modeling of problems of elasticity with unilateral discrete contact. Matematika i matematicheskoe modelirovanie MGTU im. N.E. Baumana =Mathematics and mathematical modeling of the Bauman MSTU, 2015, no. 04. DOI: 10.7463/mathm.0415.0801840 (in Russian).

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

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

4. 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.

5. 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.

6. 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).

7. 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).

8. 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.

9. 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.

10. 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.