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

УЧЕНЫЕ ЗАПИСКИ Ц А Г И Том XXI Т 9 9 О

№ 5

УДК 624.07

519.63:512.643

ПРИМЕНЕНИЕ ПОЛУАНАЛИТИЧЕСКОГО МЕТОДА КОНЕЧНЫХ ЭЛЕМЕНТОВ К РЕШЕНИЮ КОНТАКТНОЙ ЗАДАЧИ ДЛЯ ОСЕСИММЕТРИЧНЫХ ТЕЛ

Л. П. Железное, А. Е. Колмагоров

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

Конечно-элементное решение трехмерных контактных задач с использованием трехмерных объемных конечных элементов (КЭ) очень громоздко и требует больших затрат времени и памяти ЭВМ. В известных решениях таких задач рассматривались, как правило, для контактирующих тела при довольно грубом разбиении. Общую формулировку и решение трехмерной контактной задачи можно найти, например, в работах [1—3]. При рассмотрении задачи упругого контакта осесимметричных тел при неосесимметричном нагружении в целях сокращения объема вычислений возможно применение полуаналитического метода конечных элементов, который заключается в использовании кольцевых конечных элементов и представлении перемещений узлов КЭ рядом Фурье по окружной координате. Именно такой подход был предложен в работе [4]. При этом в случае заранее неизвестной границы зоны контакта решение сводилось к задаче квадратического программирования, которая в свою очередь решалась модифицированным симплекс-методом. Затраты компьютерного времени оказались более чем в 10 раз меньшими по сравнению с традиционным методом конечных элементов.

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

Рассмотрим два упруго контактирующих осесимметричных тела при неосесимметричных условиях нагружения. Каждое тело дискрети-

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

и = 2' и1 сов I ? >

г=о

I

V = ^ С08 I ¥ »

2=0

I

ге> = т151п I <р , £ [0, 2тг],

1=1

а узловые силы от внешней нагрузки соответственно в виде £ £ ра=2 р1и^о%1^, «= 2со51 ?.

1—0

1=0

(1>

где м, и и до— радиальные, осевые и окружные перемещения соответственно.

Согласно [5] полная система уравнений равновесия метода конечных элементов для отдельного тела разбивается на Ь подсистем

#Г*8'=/« 1 = 0, ..., Ь,

Рис. 1

где £ — число гармоник, К1 — матрица жесткости, соответствующая 1-й гармонике,

Я-

о,-

/=■ =

и-,

VI

да,-Р1 ■

г ш р1 ■ р' .

где п — число узлов (узловых окружностей) в теле, Ь1, Я—векторы амплитуд узловых перемещений и обобщенных узловых сил для 1-й гармоники.

Запишем уравнения совместности полных перемещений для двух контактирующих узлов

й }(<?) —и] (?)=д „;■(?),

(?) — 0? (?) :

; Д®< (?) >

(?) — ®| (?) = Дш- (?) , ? £ 10, 2тс],

(2)

где г = 1, &, £— число контактирующих пар узлов, ДЦ1- (?), Д^-/?),

■Д«>/ (?)— функции зазоров (натягов) для г-го контактного узла.

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

«/ (?) = Щ (?) + их1 (?) ,

Щ (?) = (?) + (?) >

Щ (?) = да, (?) + тх! (?),

(3)

где и, (?), vl (?), (?) — местные упругие перемещения для данного

тела, их1 (?), чюх1 (?) — перемещения от жесткого смещения тела вдоль оси х, vгi(^f) — от жесткого смещения вдоль оси г. Выразим местные перемещения контактирующих узлов через узловые контактные усилия

= £ (8/ + &) V, (4)

1=0

где

1=1

= Х А‘1р'-

/=1

Здесь А\; (3X3) — компоненты матрицы А1 = К~и, <?/—компоненты разложения узловых контактных усилий

(5)

1=0

Ъ1Рш— амплитуда перемещений контактного узла /-й гармоники от действия внешней нагрузки.

Функции перемещений от жесткого смещения выразим через величины жестких перемещений <£>х по оси х и о)г по оси z. Очевидно (рис. 1,а), что

их1 (?) == “д; cos ? > )

Ъг1 (?) = «,, (6)

wxl (?) = — w* sin ? . J

Функции зазоров (натягов) представим в виде разложения в ряд Фурье по координате

Ki (?)==2 cos /?,

/ = 0 L

— ^ Лcos/ ? *

*=о

L

Дml (?) = 2 Sin 1 9

1=1

(?)

Подставим выражения (3) и (4) в уравнения совместности (2) и с учетом (6) и (7) приходим к системе уравнений

2 (2 {А}} + а%) я1! + а! и = 2 (л‘ + 8“ - « (8>

/ = 0 1/=1 I 1 = 0

1=1,..., k,

где

а' =

при 1 = 0,

Q\ = .

toje---------- 11),

при 1-І,

Q\= 0 при /> 1, А/:

A lui д ы

ДІ

Очевидно, что в случае, если контакт по окружной координате <р не нарушается во всех контактных узлах, система уравнений (8) разбивается на Ь независимых подсистем

2 ^ = + 8^, — 8^ ■ (9)

/=і

Таким образом, имеем ЗХпХЬ уравнений 3ХПХЬ неизвестных коэффициентов разложения контактных усилий плюс формально 2X2ХЬ неизвестных жестких смещений. Фактически от жестких смещений всего две неизвестных 0)г И При 1 = 0 и два смещения ш’с И при Ь — 1. Поэтому систему уравнений совместности перемещений необходимо дополнить фактически двумя уравнениями статического равновесия по оси х и двумя по оси г. Выпишем их для одного из тел:

Л тс пи

2 / (Яч (*) С055 (?) 8ІП ?)<*?= 2 1 (Раі С0§ 9'

<=10 • - -

/= 1 о

Рщ] (?) 8ІП?) Й? »

(10)

(И>

2 | Я VI (?) = X / Ро1 (?) .

/ = 10 /=10

где <7и1, <7^, ^ — компоненты вектора С}{.

Учитывая (1) и (5) и интегрируя уравнения (10) и (11) получаем

X [Яиг-?»*) = 2 (Я^-Р4г)=Р*. /-1,

«=1 /=1

/=0-

г=1 /=1 »

Запишем полную систему уравнений совместности перемещений с учетом (12) в виде

\

(12)

А1 С?

Я1, / = 0,1,...,/.,

(13)

где

/?' = {А[ + - 8& ..., А'* + §2Л - 8&, РІ, />1, ЯІ, Р;

‘Вї, СОг

7їг 1 2 2 2)Т

С? = (С?1, Як, <»х,...............

,2 \ Т х, г г]

их, "г, ш2)

Матрица Л* симметричная и имеет следующую структуру:

А‘ =

А(ЗкхЗМ)

йгт

"(4X3 к)

ВЇ

ЗАХ4

0

При 1=0 для г-го контактного узла

0 0 0

вГ=

при I = 1

В" =

при / > 1, В? =■ 0

0 1 0

0 0 0

0 1 0

1 0 -1

0 0 о

1 0 —1

0 0 о

При решении подсистем (13) нулевые строки и столбцы, а также строки и столбцы матрицы А1, соответствующие жестким смещениям закрепленного тела, автоматически исключаются из решения системы.

Перейдем из глобальной цилиндрической системы координат ггср в местную для каждого узла, связанную с контактной поверхностью, систему ттр (рис. 1,6). Тогда система (13) преобразуется к виду

А1$ = Я1, 1 = 0, ..., Ь, (14)

где

= Я‘ = Я13, Лг = Яг5.

Здесь 5 — матрица преобразования координат.

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

Пример 1. Одноосное растяжение на бесконечности пластины с круговым отверстием, в которое впрессовано упругое кольцо. Трение отсутствует. Для этой задачи дано точное аналитическое решение в книге [6]. Расчетная схема приведена на рис. 2, где размеры даны в мм. Используются 4-угольные конечные элементы с билинейной аппроксимацией по г и г. Узловая нагрузка для данного случая содержит лишь две гармоники 1 = 0 и 1 = 2:

Р,«1 = р л п

= рн О + С08 2(Р) >

1 = ^2=- — рн эшг?,

где к — толщина пластины. Упругие характеристики для стали Ест = = 2,06-105 МПа, Vcт = 0,3, для меди £м=1,13-105 МПа, ум = 0,3. Натяг А = —Я2 соответствует значению, при котором, согласно [6], контакт-

ное давление равно нулю при ф = 0. Задача сводится к решению двух подсистем (14) для 1 — 0 и 1 = 2. Результаты расчетов приведены ниже.

В числителе приведены результаты из [6], в знаменателе значения, полученные в данной работе методом конечных элементов. Сравнение результатов свидетельствует о хорошем совпадении.

Стальное кольцо в медной пластине Д=0,020926

<7/°о

*/2

0,000

0,0015

—0,116

—0,1164

Медное кольцо в стальной пластине Д=0,010948

Ф по

0,000

0,0008

—0,036

—0,0362

Ф20

Ф15

V, 1

11111111"! Г“

Рг

Р1

Пример 2. Штифт в круглой пластине при одноосном растяжении. Расчетная схема приведена на рис. 3, а. Силы трения не учитываются.

МПа

500

250

О а- У

#3°

0

1

МПа

500

250

2,5 5,0 О

2 ММ

б)

1=0

2,5 5,0 ,

г мм

В) Рис- 3

Нагрузка такая же как в примере 1. Размеры на рисунке указаны в мм. Упругие характеристики штифта £'8 = 2,1-105 МПа, ^=0,3, для пластины .Ер= 1,05• 105 МПа, 6Р = 0,3. Штифт установлен с равномерным натягом Д = 0,02 мм. На рис. 3,6 показано распределение контактных давлений по толщине при фиксированных значениях окружной координаты <р, на рис. 3, в составляющие давлений по гармоникам. Пунктирные линии на графиках соответствуют значению давления, вычисленному по формуле для прессового соединения [7]

__ ______2ЕР ______ /1ГЧ

я — е (1_^) + (1 + >) » и о;

где е — Ер!Е8, —диаметр отверстия пластины.

Если увеличивать внешнюю растягивающую нагрузку, то в районе <р = 0 произойдет нарушение контакта и условие совместности (2) для соответствующих узлов не будет выполняться на всем промежутке от

0 до 2я. В этом случае непосредственно в данной постановке на пря-

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

итераций.

Для построения итерационного процесса удобно представить перемещения и контактные давления тригонометрическими рядами Фурье на дискретном множестве точек по окружной координате. Поскольку условия ортогональности дискретного ряда выполняются [8], то все полученные выше соотношения справедливы.

Имеем значения некоторой функции £(ср) в точках <pi = ni/N, где

1 = 0,..., N, на промежутке [0, я]. Требуется найти разложение по косинусам

ёс (Ь) = 2 С08 (Ч- *) ■

/=0 ' '

Согласно [8] при N четном коэффициенты разложения определяются по формуле

£'=^{у [£Ы + (-1)г£Ы] +

{N—2)12

+ X & № + 1)Z g ^N-i) 1cos y 1 }•

i=i >

(16)

При этом если N=1, то £ 0?;) = (?/). Дискретизация задачи

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

В качестве начального приближения берем решение, удовлетворяющее системе уравнений (14) при условии сопряжения по ср на всем промежутке от 0 до я для всех узлов

Л,(?,(0)-Я,(0), 1 = 0,..., I. (17)

На (у+1)-й итерации выполняются следующие операции:

а) определяются значения контактных усилий

# (?!> = 2 5(/-1)хз+2 сое г,

1 = 0

(18)

і = 0, ..., N; j = 1, ..., я;

б) формируются соответствующие контактным узлам вектора из значений функции

я? (ь), q{i° (<рЛ > 0

о , ?iJ)(?,)<o,

г = 0, N; j = 1, ..., я;

в) проверяется условие

/!т) '(?,) < •,

(19)

(20)

если б — малое положительное заданное число, определяющее точность решения. Если условие (20) выполняется для всех I и /, то решение задачи получено и вычисление прекращается;

г) вычисляются коэффициенты разложения

ТС/

г? ы-2 со* -1 1 = 0

по формуле (16), в которой g (<рг) = /)Т> (срг) ;

д) пересчитывается вектор правой части системы совместности

Rl^+1) =R1^ -У В1

я* (/— 1)хЗ+2 О/

вш

т ■■

(г — 1) X З Ч- 2; i=l, ..., п + 2; 1 = 0, ..., L;

е) решаем систему совместности

I (т+1).

A1 Q

и возвращаемся к пункту а.

(21)

(22)

111111 і І Г І І

Рис. 4

Рис. 5

Физически описанный выше итерационный процесс означает, что в тех точках узловой окружности, где на соответствующей итерации получается растяжение {д} (<?() > 0), прикладываются внешние нормальные силы /"л/(<рг) = <7,-(<рг) и Р% (?,) = — (?{)■ Таким

образом, контактное взаимодействие в этих точках сводится к нулю.

Пример 3. Рассмотрим нагруженный силой гладкий упругий штифт, вставленный с натягом в бесконечную пластину. Расчетная схема представлена на рис. 4. Коэффициент % вычисляется по формуле Х= (йр—й5)/йр, где йр — диаметр отверстия, — диаметр штифта, £8 = 2,1-105 МПа, е=ЕР/Е8, тр=т8 = 0,3, Ер, Ев, \Р, Vs — модули упругости и коэффициенты Пуассона соответственно штифта и пластины. На рис. 4 показана зависимость нагрузки от размеров зоны отрыва при различных значениях параметра е. Сплошные кривые получены в работе [9] аналитическим методом в плоской постановке. Кружочками отмечены результаты, полученные авторами при Ь = 30. При е=1 имеет место хорошее совпадение. При е=1/3 и 3 очевидно расхождение результатов, которое объясняется тем, что в [9] неправильно взята формула, соответствующая формуле (15). Эта формула записана в [9] следующим образом:

На рис. 5 показано распределение контактного давления при различном числе гармоник. Из рисунка видно, что при увеличении числа гармоник решение сходится. Причем для <р<75° решение при Ь = 6 практически совпадает с решением, когда /. = 30. Значение 0Р = 95° по-

лучено в работе [9] для р = оо. Следует отметить, что относительная погрешность итерационного процесса ('16) — (21) составляла менее 0,1 % при десяти итерациях.

ЛИТЕРАТУРА

1. Хворостухин Л. А., Шишкин С. В. Общий метод решения трехмерных конструкционно-контактных задач. — Проблемы прочности,

1985, № 1.

2. L е е В. С., К w а к В. М. A computational method for elastoplas-tic contact problems. — Computs Structs, 1984, N 18.

3. Tanaka М., Mivazawa H., Asaba E., Hango К- Применение метода конечных элементов к болтовым соединениям. Фундаментальные исследования по анализу болтовых соединений с использованием метода конечных элементов. — Jap. Soc. Moch. Eng., 1980, с. 46, № 410.

4. Park S. G., Kwak В. M. A semi-analytical finite element method for three-dimensional contact problems with axisymmetric geometry. —

Proc. Inst. Mech. Eng., 1986, c. 200, N 6.

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

1975.

6. С а в и н Г. Н., Т у л ь ч и й В. И. Пластинки, подкрепленные составными кольцами и упругими накладками. — Киев: Наукова Думка,

1971.

7. Прочность, устойчивость, колебания. Том 2.—Справочник. — /Под ред. Биргера И. А., Пановко Я- Г. — М.: Машиностроение, 1968.

8. Хемминг Р. В. Численные методы. — М.: Наука, 1972.

9. Г х о ш С. П., Д а т т а г у р у Б., Р а о А. К. Передача нагрузки от гладкого упругого штифта к пластине больших размеров. — РТК,

1981, т. 19, № 7.

Рукопись поступила 20/1 1989 г.