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

электронное научно-техническое издание

НАУКА и ОБРАЗОВАНИЕ

Эя № ФС 77 - 30569. Государственная регистрация №0421100025. ISSN 1994-0408_

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

77-30569/353180 # 04, апрель 2012

Станкевич И. В., Яковлев М. Е., Си Ту Хтет

УДК 519.3

МГТУ им. Н.Э. Баумана [email protected] [email protected] [email protected]

Введение

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

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

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

задач в рамках конечно-элементной технологии используют алгоритмы, основанные на метод штрафных функций, методе множителей Лагранжа, комбинировании методов штрафных функций и множителей Лагранжа [2, 4], имеются работы, в которых используется понятие "псевдосреды" [1]. В данной работе рассматривается алгоритм решения контактной задачи теории упругопластичности, базирующийся на альтернирующем методе Шварца [9], основанном на раздельном рассмотрении контактирующих тел и принципе поочередности. Этот метод позволяет отказаться от значительной части ограничений, накладываемых на используемые конечно-элементные модели (сетки); не возникает необходимости в переформировании глобальных матриц жесткости при изменении геометрии зон контакта; метод отличается относительно быстрой сходимостью.

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

Рассмотрим два двухмерных однородных и изотропных линейно-упругих контактирующих тела А и В, занимающих на плоскости (в пространстве Я2) области Ол и Ов и ограниченных кусочно-гладкими границами д0А и д ОВ . Введем на плоскости

декартовую систему координат Оx1x2 (рисунок 1).

Рисунок 1 - Схема контактного взаимодействия двух тел

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

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

Ъ,; (и,Т) + (х) = 0, х еаа , I,] = й; а е {А,В}; граничные условия (кинематические и силовые соответственно):

и (х)|^ = и0 (х), х е^сдО,; ае {А, В};

Ъ (и, Т) п\ = рг (х), х е^2аедаа, г, ] = \2;ае {А, В};

соотношения Коши:

% (х) = 1 (Иг, }■ (х) + и}■ ,г (х)) , х е °2 , г, ■ = 1, 2 ;2 е {Л В};

(1)

(2)

(3)

(4)

и определяющие уравнения - закон Гука в виде (начальные напряжения отсутствуют):

ъ (х) = Н (е (х) - е (х)), х е аа, а е {А, В}, (5)

где х = ■< 1 I - вектор координат, Н - матрица Гука, ъ = {ъ} =

'22

'12

вектор

напряжений, е = {е} =

"11

е22 Г12

- вектор деформации, где уп = 2е12, е0 = {е0} =

Г12

вектор начальной деформации, 0 - компоненты массовой (объемной) силы

«х) - {«(■ >}=й'«х )=К I!

вектор заданных перемещений точек

поверхности , pi (х) - компоненты заданной распределенной нагрузки р (х) =

ГР1 (х У

Р2 (х )\

на поверхности £2.

Компоненты матрицы Гука зависят от вида рассматриваемого напряжено-деформированного состояния. Если расчет проводить по схеме плоского напряженного состояния, то матрицу Гука следует записать в виде:

" 1 ¡л 0

Л 10, (5 а)

_ 0 0 (1 -¡)/2_

Н = [Н] = т

Е

¡

в случае плоского деформированного состояния имеем:

1 //(1 -/) О

н = [н]- е(1 -/)

(1 + /)(1- 2/)

/(1 -/) 1 О

О О (1 - 2/)/2 (1 -/)

(5 б)

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

< (х)- ивп (х) = дя (х), X е Бк; (6)

и по напряжениям (силовое условие):

сЯА (х)=-авп (х )< О, х е Бк; (7)

А В »-»

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

Совокупность соотношений (1) - (7) составляет математическую формулировку контактной задачи теории упругости. Для ее решения в данной работе был использован алгоритм, основанный на альтернирующем методе Шварца [9].

2. Основные процедуры альтернирующего метода Шварца

Альтернирующий метод Шварца является итерационным методом. Его суть состоит в следующем. На первом шаге на контактных поверхностях тел А и В

соответственно и задают начальные перемещения и (х)| = и°А (х) и

и (х) В = и°В (х) , которые имеют смысл дополнительных кинематических условий.

Значения компонент векторов и0А (х) и иВ (х) выбирают из диапазона ожидаемых значений для зоны контактного взаимодействия.

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

2(а)), но могут и не совпадать (рисунок 2(б)). Это зависит от особенности постановки решаемой контактной задачи. В конце решения контактной задачи геометрически скорректированные поверхности контакта совпадают = .

Далее решают независимо две подобные задачи теории упругости для тел А и В . Затем вычисляют поверхностные силы рА (х) и рвк (х) на контактных поверхностях и

Б В и их корректируют так, чтобы выполнялись силовые контактные условия (7).

На втором шаге на контактных поверхностях БА и БВ задают силовые контактные

условия, в качестве которых используют скорректированные поверхностные силы рА (х)

и рВ (х) :

а; (и, Т) пЧ\а= ра (х), х е^сЭб . а = А,В, /, ] = 1,2; (8)

Рисунок 2 - Начальная геометрия контактных поверхностей (а - контактные поверхности совпадают, б - контактные поверхности не совпадают)

Далее вновь решают независимо задачи теории упругости отдельно для тел А и В . По результатам полученных решений, выполняют коррекцию компонент векторов перемещений иА (х) и иВ (х) соответственно точек контактных поверхностей и Б^В с

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

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

Скорректированные перемещения контактных поверхностей иА (х) и иВ (х) рассматривают в качестве новых кинематических граничных условий на геометрически

измененных поверхностях контакта и 5кВ и вновь решают независимо задачи теории упругости отдельно для тел А и В :

и (х)| = < (х), х е Б2.2 = А,В ; (9)

Таким образом, данный алгоритм состоит в реализации итерационного процесса поочередного задания на поверхностях контакта и векторов перемещений ^ (х) и

икВ (х) и векторов поверхностных сил рА (х) и рВ (х), а также в их соответствующей

коррекции с тем, чтобы были выполнены либо силовые контактные условия, если в зоне контакта заданы перемещения, либо кинематические, если заданы поверхностные силы. Вопросы сходимости подобного типа алгоритмов рассмотрены в работах [8,9].

Процедуры коррекции векторов перемещений иА (х) и и^В (х) и векторов сил

ркА (х) и ркВ (х) рассмотрены ниже в разделе 4.

3. Вариационная постановка линеаризованных задач теории упругости и

построение СЛАУ

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

Выражение функционала полной потенциальной энергии линейно-упругого тела, размещённого в пространстве №2 и нагруженного массовыми «(х), поверхностными

р(х) и дискретными (сосредоточенными) Яп (п = 1 ,Кк) силами, имеет вид [2 - 3]:

1 К„

П =1 |(е-ео)То СУи- |итиX-р\итр^, (10)

2 О О 5, п=

( ) \и (х)1 " ( ) \и (хп )1 К1

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

[V (х)[V (хп)К

вектор перемещений фиксированных точек хп тела, в которых заданы дискретные силы

Яп (п = 1 ,Кк); в = 0, если на контактной поверхности 5к заданы перемещения ик (х),

в = 1, если на контактной поверхности 5к заданы поверхностные силы рк . Вариационная задача эквивалентная задаче (1) - (5) , (9) имеет вид

(11)

П ^ min;

u(x= u0 (x); x £ S1 e dG;

u (x )l S =uk (x);x £Sk edG ^ Sk

Вариационная задача эквивалентная задаче (1) - (5) , (8) может быть записана следующим образом:

П ^ min;

u (x= u0 (x); x £ S1 e dG.

Для решения вариационных задач (11) и (12) используем метод конечных элементов.

С учетом закона Гука (5) функционал (10) можно представить следующим образом:

(12)

1 nR

П = -j(s-so)T H(s-so)dv- juTQdv- juTPds- SukRk -ßjuPkds =

2 G G S2 k=' Sk

11 Пк = — jsTHsdv - jsTHs0dv + —js^Hsd - j uTQdv - j uTpds S uTRk - ßj uTpkds.

(13)

~ о о ~ о о я2 к=1 Бь

Введем обозначения для интегралов, стоящих в правой части последнего равенства:

Ц = 2 jsTHsdv; П^ = -jsT Hs0dv; Ц = 1 jsT HS)dv; Ц = -j uTQdv;

2 G G 2 G G

KR

Ц = - juT Pds; nR = -S uk Rk;Ц = - juT Pkds.

(14)

Тогда получим в виде:

п = ц. + п. + ц. +пв +Пр + П + вП. (15)

Минимизацию функционала (15) выполним с помощью метода конечных элементов. Необходимым условием достижения функционалом (15) экстремума является выполнение условия

¿а+^ + ^ + ^ + ^ + ^ + в-^ = {0}. (16)

d {U } d {U } d {U } d {U } d {U } d {U } d {U } d {U }

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

{и } =

Ul U 2

и

и

2хK„

(17)

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

рассматриваемого тела.

После выполнения всех конечно-элементных процедур соотношение (16) сводится к матричному уравнению (глобальной системе линейных алгебраических уравнений -СЛАУ), которая имеет вид:

[ к ] {и } = {*}.

(18)

4. Алгоритм численного решения контактной задачи теории упругости на основе альтернирующего метода Шварца

Рассмотрим подробно работу алгоритма на примере контакта двух тел А и В. Пронумеруем узлы контактных поверхностей (а = А, В), тогда векторы

соответственно перемещений контактных узлов {и^ } и узловых сил {^ } в них можно записать соответственно в виде:

К }а

К }(а),1 К }(а),2

К }

(а)Ма

(а),1

(а),2

(а)Ма

(19)

К }(а) =

1^к} 1

I к )(а),\ («к }

I к )(а),2

(«к }

(«)м*

(а),1

(а),2

(а)Ма

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

контактных сил в узле т (1 < т < Ма ) .

Последовательно для каждого тела в соответствии с алгоритмом выполняется коррекция компонент векторов (ик} ^ и («} ), а = А,В . Коррекция компонент

векторов (и}( ^ и }( ), а = А,В реализуется поочередно: на каждой четной итерации

корректируются компоненты вектора (ик } , а на нечетной итерации - (« } .

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

(и к }

(А)

выполняется в соответствии с соотношением:

(и}

2п Iй

к /(А),т

2п

(А),т

и

\и\ IV

п = 0;

(А),т 2п-1

(А),т

+ а

2 п—1 (А),т

(( Л 2п—1

и I

(В Ъ

и

IV

2п—1 Л

(А),т

(21)

п = 1, 2,...

? * ' * 5

где а2А" т - итерационный параметр, т (1 < т < М^ ) - номер текущего узла, лежащего на контактной поверхности Б А тела А , здесь М А - число контактных узлов на

2п—1

оА

поверхности Б ;

к ; л > - вектор перемещений сходственной точки ^, лежащей на

^( В),

контактной поверхности тела В .

Значения компонент и0 и V0 (п = 0) в узле т выбирают из ожидаемого диапазона значений на основании априорной информации о характере контактного взаимодействия.

Для коррекции компонент вектора контактных узловых сил («используется

формула:

2п+1 г л 2п ( г л 2п с л 2п Л

^=ЬI,= и !•• -ай"

К}...... „

(А),т ^ (А),т

1П [

Ч Я (В),* IЧ I,А),т J

п = 0,1,2, ..., (22)

где аАП)т - итерационный параметр, т ,1 < т < Ма ) - номер текущего узла, лежащего

2 п

на контактной поверхности Б^А тела А, \ > - вектор контактных узловых сил

УЧ(в

сходственной точки ^, лежащей на контактной поверхности Бк тела В .

Чтобы воспользоваться формулами (21) и (22) необходимо знать компоненты

вектора перемещений ив ={и}(_ =-11 и вектора сил ЯВ = {^}(в) = •! г в

(В ^ ^ 1( в(В ^ Уч 1( в

сходственной точке ^, лежащей на контактной поверхности БВ тела В. Существуют разные подходы к построению сходственных точек [9], например, в качестве сходственной точки ^ можно использовать точку пересечения перпендикуляра, опущенного из узла т , лежащего на контактной поверхности БА тела А, на контактную поверхность БВ тела В, состоящую из одномерных конечных элементов (граней). При этом точка пересечения должна находиться в пределах соответствующей грани. В данной работе в качестве сходственной точки ^ рассматривалась точка пересечения отрезка т0 т, соединяющего начальное т0 и конечное т положения узла, принадлежащего

контактной поверхности БА тела А, с отрезком г] , лежащем на контактной поверхности БВ тела В (рисунок 3).

Рисунок 3 - Построение сходственной точки ( £А0 - начальное положение контактной поверхности тела А )

После того, как сходственная точка ^ определена, то есть найдены ее глобальные координаты Х15 и Х2$ и грань ¡] контактной поверхности БкВ тела В, на которой лежит сходственная точка, на грани ¡] строится внешняя нормаль пВ и вектор касательной тВ к контактной поверхности БВ тела В (рисунок 4):

пВ = (Х2, — Х2г ) е1 +(Х. — Х1, ) е2 (23)

ТВ = (Х1г — Х1} ) е1 +(Х2 ] — Х2г ) в2 (24)

4 '

где - длина отрезка ¡] .

Рисунок 4 - Построение векторов нормали пВ и касательной тВ в сходственной точке

Затем в узле т вычисляется значение контактной силы ЯАВ в направлении

п ,т

внешней нормали nB , имеем:

RnB,m = PmC0SP + , (25)

здесь pm и qm - компоненты вектора контактных сил {Rk} = \ P I в узле m

V J( A\m q

Я J (A),m

(x ^ X )

(1 < m < Ma ) ; cos^ = cos(nB,e1 ) = -.

dj

Далее необходимо проверить выполнение условия:

RA >0. (26)

п т

Если условие (26) не выполняется, тогда узел m исключатся из описания контактной поверхности SA тела A на данной итерации, а при выполнении условия в сходственной точке s вычисляются компоненты вектора перемещений

UB = {U}( , = \ I (если итерация четная) или вектора сил RB = {R}(r) = <! P I (если

(B),s lvj(b),s (B),s Iq j(b),s

итерация нечетная). При этом используются следующие соотношения:

\us = N8) (Z) uB + Nf) (Z) uB, (27)

Vs = N<8) (Z) vB + Ng) (Z) vB, ( )

Ps = ^8) (Z) pB + N8) (z) PB, (28)

q = N8) (Z) qB + N(8) (Z) qB, ( )

где N(8)(Z) и N(8)(Z) - функции формы одномерного конечного элемента (грани), аппроксимирующего часть контактной поверхности SB тела B (рисунок 4); Z = dis/dj , dis - длина отрезка is, dtj - длина отрезка i j .

Вычисленные значения компонент вектора перемещений UB = {U}

(B ),s

V

используются в итерационных соотношениях (21) для коррекции компонент вектора перемещений -II , а компоненты вектора силы ЯВ = {Щ^ ^ =] г -соответственно в итерационных соотношениях (22) для коррекции компонент вектора

{л 2 n+1

> , узла m контактной поверхности SB тела A .

q J (A),m

Для компонент векторов перемещений и сил узлов контактной поверхности БкВ

тела В выполняется совершенно аналогичная процедура коррекции.

Заметим, что после выполнения коррекции компонент векторов перемещений узлов контактных поверхностей БкА и БкВ , то есть на четных итерациях, узлы контактных

поверхностей располагаются на общей контактной поверхности Бк .

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

2k

U ,Л = i

n,( A),m

u0 (/l) , k = 0;

/ ч (29)

2 k-1 . 2 n-1 / 2 k-1 2 k-1 с \ / 1 о

+ O.N u - U * A\ -o , k = 1,2, . ..,

(A),m X n,(B),s n,(A),m P •>•>•>

Un,( A),m

где a(2An)m - итерационный параметр; m (1 < m < Mд ) - номер текущего узла, лежащего на контактной поверхности SkA тела A , здесь M A - число контактных узлов на поверхности Sд; un (A) m - перемещение узла m , лежащего на контактной поверхности Sд тела A, в направлении внешней нормали n к поверхности тела B, построенной в сходственной точке s, лежащей на контактной поверхности SB тела B (рисунок 4);

OB

un(B)s - перемещение сходственной точки s, лежащей на контактной поверхности Sk тела B, в направлении внешней нормали n к поверхности тела B ; S - начальный зазор между контактными поверхностями SA и Sf. Коррекция контактных усилий и в этом случае проводится с помощью соотношения (22).

5. Вычисление итерационных параметров

Для эффективной работы рассмотренного алгоритма важным является вычисление численных значений итерационных параметров a^h и Оо)m ,а = A,B , 1 < m <Ma , здесь

Mа - число контактных узлов тела а .

В работе [8] показано, что итерационные параметры удовлетворяют условиям и между ними существует связь:

0 < а(2"-1, а(2" < 1 (30)

(a),m5 (a),m v ^

а?"- + а(2" = 1. (31)

( а ),m ( а ),m ^ '

Итерационные параметры можно принять постоянными в пределах контактных поверхностей SО ( а = A,B) и не изменять их значения при проведении итераций,

например, положить а^- т = а^ т = 0,5 . Однако, как показывают многочисленные

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

пары узлов (т, ') : узел т (1 < т < Ма ), рассматриваемой контактной поверхности ^^ ( а = А, В), и соответствующая ему сходственная точка ', лежащая на контактной

поверхности (в = А, В), а ^ в .

В данной работе с учетом результатов численных экспериментов был использован следующий подход. Пусть для определенности а = А и 1 < т < МА . На четных итерациях,

2 \п 12"

то есть когда выполняется коррекция компонент вектора перемещения и^, т = 1 г ,

IУ(А),т

итерационный параметр а2^ определяется с помощью выражения:

Ци2-1!

2 п-1 __|| т || /"ЗОЛ

а( А),т II 2 п-1 II II 2 п-111 , (32)

1 2 п-1 II , II 2п-1||

и + и

т

где ||итп ^ - норма вектора перемещения и2т" 1 узла т (1 < т < Ма ), лежащего на контактной поверхности 8А тела А, ||и2п-1|| - норма вектора перемещения и'п-1

сходственной точки ', лежащей на контактной поверхности тела В . Очевидно, что выполняются соотношения:

__2 п-1 __2 п-1 ^л __2 п-1 .__2п-1 1 /ллч

< а(А),т , а(В),1 < 1 и а(А),т + а(В- = 1 . (33)

Заметим, что алгоритм вычисления величины а(В- не требует.

Соотношение (32) можно записать для компонент вектора перемещения и

2п-1 : т

и2-1

а2п-1 = | п,т |

(А),т,и ~ I 2п-11 . 2п-11 '

и + и

П,т П, А'

I 2 11 (34)

V

2 п-1 = | Т,т |

(А),тV _ V2п-11 + V 2п-1| ,

Т,т т,а

2 п—1 2 п—1 2 п—1 2 п—1 2п—1

где ипт , vTm ,ип' и vTs - соответственно компоненты векторов перемещений ит и и2"-1 в направлении внешней нормали пВ и вектора касательной тВ в сходственной точке ', лежащей на грани г] на контактной поверхности 8В тела В (рисунок 4)

На нечетных итерациях, то есть когда выполняется коррекция компонент вектора контактных узловых сил {Як} , итерационный параметр 0С(А)т определяется с помощью

выражения:

2" Ят

II Т>2" \\Ят я2 "||

а(А),т = II —II II „ч„П' (35)

где \\Ят"\\ - норма вектора контактных сил Я^ в узле т (1 < т < М ^ ), лежащего на контактной поверхности 8^ тела А, Ця^Ц - норма вектора контактных сил Я2" в сходственной точке ^, лежащей на контактной поверхности 8В тела В .

Соотношение (35) можно записать для компонент вектора контактных сил Я^ :

(А ),т, р

а( А),т,д

|р2п 1

|Р2" 1 + \1г п,т и

q2n

#1'

(36)

2" 2" 2" 2" г%2"

где Рпт, ЧТт, Р„ц и - соответственно компоненты векторов контактных сил Ят и Я2" в направлении внешней нормали пВ и вектора касательной тВ в сходственной точке s, лежащей на грани ¡] на контактной поверхности 8В тела В (рисунок 4)

Как показали результаты численных исследований использование для вычислений итерационных коэффициентов а0—т и а(2а"т соотношений (32) и (35) является достаточно эффективным.

6. Учет трения при решении контактных задач

В данной работе было принято, что трение описывается законом Амонтона-Кулона. Отсюда учет трения между контактирующими поверхностями тел А и В (соответственно и 8В ) требует в каждом контактном узле т проверки неравенства:

< РРпт, 1 < т < ма, а = А, В, (37)

где / - коэффициент трения; рпт и qTm - компоненты вектора силы Ят в направлении

внешней нормали пв и вектора касательной тв в сходственной точке ^, лежащей на грани ¡] на контактной поверхности 8е тела в, в = А, В и в Ф а (рисунок 5).

Рисунок 5 - Учет сил трения

Векторы нормали пв и касательной тв определяются выражениями:

пв =

тв =

(Х2] — ^) , Х — Х] )

е1 +

2

(Х1 ] — Х1г ) ,(Х2] — Х2г )

(38)

^ +

где = ^Х1} — Х1г) + (Х2} — Х2г) - длина отрезка (рисунок 5).

Обозначим через рпт и цТт компоненты вектора силы и через ипт и vTm компоненты вектора перемещения в узле т , 1 < т < Ма , а = А, В .

Если имеет место строгое неравенство (37), то считаем, что в узле т, например, тела А наблюдается полное сцепление с контактной поверхностью 8В тела В в

сходственной точке ' . В этом случае следует принять:

V = V

тт та 5

(39)

где vTs - компонента перемещения сходственной точки а , расположенной на контактной поверхности 8В тела В, в направлении вектора касательной тв (рисунок 5, а = А,

в = В ).

Нарушение неравенства (37) вызывает взаимное проскальзывание контактирующих тел, поэтому необходимо использовать соотношения, накладывающие ограничения на касательную компоненту силы в узле т , 1 < т < Ма , а = А, В, имеем:

qтm =-МРпт^ёП (Т ) . (40)

Таким образом, в начале вычислений с учетом трения все узлы т , 1 < т < Ма ,

рассматриваемой контактной поверхности , а = А, В, считаются сцепленными. С этими граничными (кинематическими) условиями выполняется начальная (нулевая) итерация. После усреднения усилий в зоне контакта проверяется выполнение силового условия (37). В случае невыполнения условия (37) соответствующие узлы переводятся в разряд скользящих и в них ограничиваются касательные усилия по формуле (40). Далее итерационно уточняются компоненты векторов перемещений и сил в контактных узлах по формулам (21) и (22). Практика решения задач показала, что учет трения не увеличивает общего числа контактных итераций.

Также учитывается возможность упругопластического деформирования. Для этого используется метод переменных параметров упругости [7].

7. Анализ результатов численных исследований

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

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

При расширении тел начальный зазор был полностью выбран. На рисунке 7 представлено поле компоненты а22 тензора напряжений.

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

Рисунок 6 - Расчетная схема контактного взаимодействия равномерно нагретых тел

Рисунок 7 - Поле компоненты а22 тензора напряжений, МПа.

Исследовано контактное взаимодействие двух неравномерно нагретых прямоугольных тел (рисунок 8), выполненных из стали 40Х. Размер верхнего тела вдоль оси абсцисс 0.025 м, нижнего - 0.05 м. На верхнее тело действует распределенная горизонтальная нагрузка р2 = 100 МПа. Оба тела подвержены температурному

воздействию. В пределах верхнего тела температура меняются от оси закрепления Ох2 вдоль оси Ох1 до свободной поверхности справа по линейному закону. Максимальное значение температуры равно 600 К и относится к точкам тела, лежащим на оси Ох2, а точки, расположенные на свободной поверхности имеют минимальное значение температуры - 200 К. У точек нижнего тела, для которых 0 < х1 < 0.025, аналогичное

распределение температуры, а точки, у которых 0.025 < х1 < 0.050, имеют постоянную температуру 200 К. Обобщенный график распределения температуры по обоим телам показано на рисунке 9.

Рисунок 8 - Расчетная схема контактного взаимодействия неравномерно нагретых

прямоугольных тел

Рисунок 9 - Распределение температур в контактирующих телах

На рисунке 10 показана геометрия конечно-элементных моделей контактирующих тел после деформирования. Из характера деформирования видно, что температурное воздействие является определяющим. На рисунке 11 приведено поле компоненты а22 тензора напряжений. Данное поле характеризуется неоднородностью.

Рисунок 10 - Геометрия конечно-элементных моделей тел после деформирования

- 497

-184

i I A i ■ 128

Рисунок 11 - Поле компоненты а22 тензора напряжений, МПа.

Выводы.

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

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

1. Бабин А.П., Зернин М.В. Конечно-элементное моделирование контактного взаимодействия с использованием положений механики контактной псевдосреды // Механика твердого тела.2009. №4.с. 84-107.

2. Wriggers P., Nour-Omid B. A two-level iteration method for solution of contact problems // Computer methods in applied mechanics and engineering. 1986. №54. pp. 131-144

3. Kikuchi N., Oden J. T. Contact problem in elasticity: a study of variational inequalities and finite element methods. - Philadelphia: SIAM, 1988. 495 p.

4. Wriggers P. Computational contact mechanics. - England: John Wiley & Sons Ltd, 2002 -441 p.

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

6. Шабров Н.Н. Метод конечных элементов в расчетах деталей тепловых двигателей — Л.: Машиностроение, 1983. 213 с.

7. Малинин Н.Н. Прикладная теория пластичности и ползучести — М.: Машиностроение, 1975. 398 с.

8. Цвик Л. Б. Принцип поочередности в задачах о сопряжении и контакте твердых деформируемых тел // Прикладная механика. 1980. т. 16. № 1. с. 13-18.

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

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

electronic scientific and technical periodical

SCIENCE and EDUCATION

_EL № KS 77 - 3Ü56'». .V;II421100025, ISSN 1994-jMOg_

Mathematical simulation of contact interaction of elasto-plastic environments

77-30569/353180 # 04, April 2012

Stankevich I.V., Yakovlev M.E., Si Tu Htet

Bauman Moscow State Technical University [email protected] [email protected] [email protected]

The authors consider an algorithm of solving of a contact problem of the theory of elasto-plasticity on the basis of alternating Schwarz method. The article includes a numerical solution of the theory of elasticity contact problem with the use of implementation of an iterative process of the alternate appointment of displacement vectors and surface forces on surfaces of the contact, and also their corresponding correction. The authors consider the mechanism with account of friction and elasto-plastic deformations when solving contact problems with the use of the developed algorithm. Results of numerical simulation of contact interaction of heated bodies are included in the article. The geometry of finite-element models of contacting bodies after deformation is shown.

Publications with keywords: theory of elastity, tensor of tensions, alternating Schwarz's method, elasto-plasticity, iterative process, contact task

Publications with words: theory of elastity, tensor of tensions, alternating Schwarz's method, elasto-plasticity, iterative process, contact task

References

1. Babin A.P., Zernin M.V. Konechno-elementnoe modelirovanie kontaktnogo vzaimodeistviia s ispol'zovaniem polozhenii mekhaniki kontaktnoi psevdosredy [Finite element modeling of contact interaction with the use of the provisions of the contact mechanics pseudomedium]. Mekhanika tverdogo tela, 2009, no. 4, pp. 84-107.

2. Wriggers P., Nour-Omid B. A two-level iteration method for solution of contact problems. Computer methods in applied mechanics and engineering, 1986, no. 54, pp. 131-144.

3. Kikuchi N., Oden J.T. Contact problems in elasticity. Philadelphia, SIAM, 1988. 495 p.

4. Wriggers P. Computational contact mechanics. England, John Wiley & Sons Ltd, 2002. 441 p.

5. Zenkevich O. Metodkonechnykh elementov v tekhnike [Finite element method in engineering]. Moscow, Mir, 1975. 542 p.

6. Shabrov N.N. Metod konechnykh elementov v raschetakh detalei teplovykh dvigatelei [The finite element method in calculations of parts of heat engines]. Leningrad, Mashinostroenie, 1983. 213 p.

7. Malinin N.N. Prikladnaia teoriiaplastichnosti ipolzuchesti [Applied theory of plasticity and creep]. Moscow, Mashinostroenie, 1975. 398 p.

8. Tsvik L.B. Printsip poocherednosti v zadachakh o sopriazhenii i kontakte tverdykh deformiruemykh tel [The principle of sequential order in the tasks of the conjunction and contact of solid deformable bodies]. Prikladnaia mekhanika, 1980, vol. 16, no. 1, pp. 13-18.

9. Mozharovskii N.S., Kachalovskaia N.E. Prilozhenie metodov teoriiplastichnosti ipolzuchesti k resheniiu inzhenernykh zadach mashinostroeniia: V 2 ch. Ch. 2: Metody i algoritmy resheniia kraevykh zadach [The application of the theory of plasticity and creep to the solution of engineering problems of mechanical engineering: In 2 Pts. Pt. 2: Methods and algorithms for solution of boundary value problems]. Kiev, Vysshaia shkola, 1991. 287 p.

10. Stankevich I.V., Iakovlev M.E., Si Tu Khtet. Razrabotka algoritma kontaktnogo vzaimodeistviia na osnove al'terniruiushchego metoda Shvartsa [Development of algorithm of contact interaction on the basis of Schwarz alternating method]. VestnikMGTUim. N.E. Baumana. Ser. Estestvennye nauki [Bulletin of Bauman MSTU. Ser. Natural science], 2011, Spets. vyp. «Prikladnaia matematika» [Spec. iss. "Applied mathematics"], pp. 134-141.