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

М Е Х А Н И К А

УДК 539.3

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

Докт. физ.-мат. наук, проф. ЧИГАРЕВ А. В., докт. физ.-мат. наук МЕЛЕШКО И. Н.,

ПРОНКЕВИЧ С. А.

Белорусский национальный технический университет

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

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

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

Эх Эv

Эv Эх

Рис. 1. Внутренний контакт областей с круговыми границами

Закон Гука имеет вид:

о" = 2G

1-2ц

о" =2 G

е™ +——е а

1-2 ц

а _ (-1 и „ и

\vv — — '7 £.w •

Формулы Коши записываются в виде:

Э и

р =- р = -

Эх ' v

а

ЭУ

Эу

а ди dv а а , а

Эу Эх

(2)

(3)

Наука итехника, № 6, 2013

где индекс а = 1 относится к цилиндрической полости, а = 2 - к цилиндру. Примем: Г\, г2 -радиусы отверстия полости и цилиндра соответственно; е = г2 - Т\ - радиальный зазор.

Переходя к комплексному потенциалу, представим перемещения и, V и напряжения в полярных координатах в виде [2]:

o¡- +аг =2

ф01 w + фа W

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

г\

£ + — Е

2d,

Oj Г) +Oj Г)

" G11+V1G21 Or

+

'эс

kx +1 4 щ

Ф, л -<Í\ Л

_ ' 2

- о° + 2i xarQ = 2е2,(- [w Ф' w + W w ]; (4) " £

2а.

Ф2 с +Ф2 С

- g12+v2g22 о,.

+

2|i ua+iva =kaфа w -wФа w w ,

где - коэффициенты Ляме для тела с полостью (а = I) и для цилиндра (а = 2); фа w ,

\|/а w — комплексные потенциалы Колосова -Мусхелишвили:

ф«' w =ф" w ; xf: W W .

Кроме того, имеют место выражения [3]:

Ф, w = -

kx iY 1 г а,, со da 2к \ + kx w 2ni I со-w

Ф, » =■

-i Y 1

iY w

2k l + k2 w к \ + k2 r2'

1 го,- % d%__ra,. % d^

2ra' l - w 4кг [

- +

+ Л

ЭС

k2+1

4ц,

Ф2 ti -Ф2 л

, (6)

где \/L, = rl/ r2ri ; £ = r,T| /rx—m.

Таким образом, из (5) и (6) получаем инте-гро-дифференциальное уравнение в виде [2]

Г| г о,, со ¿/со

ra¿ со-т)

= Yiar Л

/у --Y2

71

Т| i?2

Г ,

-Уз--Y46-Y5^

к

где

Yi ^

G12-v2G22 G,, -v,G2l

2 Tj^Gjj+Г2^С12

y2

l + v2 Elrlr2+kl 1 + Vj E2rf

4 r{'E2Gll+r2ElGl2

,2

,2

w w

Л.2 Л

w

V У

r~ ' -^Фа w w

где Y - главный вектор сил, приложенных к границе.

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

8 + г/j cos 'Q +1, sin 'С, = и2 cos 'С, +

+ v2 -8-е sin С, , (5)

где u1, Vi - компоненты перемещений для тела с цилиндрической полостью; u2, v2 - для цилиндра; 5 - осадка центра цилиндра.

Наука итехника, № 6, 2013

у =_G12£1¿i_

2 г, 1 + k2 r{íE2Gl j + r2E2Gx 2

у - о

4 KE2GU + Г2Е1012

ъ =

ЕхЕ2

1 ra

2 rlE1Gn +r1E]Gn ' rx

- f—с/со;

i j ril

2717 l СО

Ti = Re/'^.

Для дальнейшего рассмотрения удобно преобразовать полученное выше интегральное уравнение к виду [4]

Г х 1 ^Г' ^

----dt=f х,-1<х<1, (7)

ТГ '

В х к J, t - х

b

ri Y

где Г i =о w ; В х - —; fx---у2х

V

7ГГ)

X

+ —

Y

у

п

- известные

Г) Я2

ч

функции; Г(х) =ог - искомая функция.

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

Г(-1) = а; Г(1) = р. (8)

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

Преобразуем уравнения Прандтля к интегральному уравнению специального вида. Обозначим

1 'г Г" I

тг J

-их

71 ^ t - X

и применим формулу обращения сингулярного интеграла с ядром Коши [6], тогда получим

Г' х =

1

Г и t

где С - произвольная постоянная и

■Г

Г х = [ Г' х & = — \ и t Н х^ Л + \ л:-',

-dt + -

С

А ■ к

+С arcsinxH— 2

-с.

(9)

где Н x,t =1п

l-xi + ^j h-x2 1-t2

\t-x\

С -

произвольная постоянная.

С помощью условия (8) находим, что

С, =а: С =

к

Таким образом, краевая задача (7), (8) свелась к интегральному уравнению относительно функции и(х)

и х +-

ц

> v J

пВ х

и t Н х J dt - F х

где

F х - f x +

В x

г

a-p

к

arcsin x + —

V 2 У

-a

Обозначим q x — — X и перепишем

В x

последнее уравнение в виде

1

q х

и(х)--, \u(t)H(x,t)dt = F(x). (10)

ллД "2

-X2 I:

Введем оператор

К и -К и,х

¡/ х

тгл/i-

х2 Ч

Jw t Н x,t dt. (11)

Тогда интегральное уравнение (10) приводится к функциональному уравнению

Из УТ и < max

11 НС .те -1,1

и-К и — /'. | q х |

K<Jl-X2 -1

|я x,t dt

(12)

iml

следует, что q(x) и Ах) непрерывны на отрезке [-1, 1]. Тогда с учетом (11) получаем неравенство

|К и Ц^ тах^ х||Н|с.

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

Зададим на отрезке [-1, 1] систему точек

2

xk = kh\ к = -п, ..., -1, 0, 1, ..., и; h =

2п + \

и аппроксимируем функцию и(х) на этом отрезке по формуле

п

м(х) = й(х) = ^едДх)м(хД

А А

в которой 0д,(х) = 1, если хе

9Д, (х) = 0, если х I

h

Xl X*. "h 1

2 2

1 2 2

h

Наука итехника, № 6, 2013

1

-l

В результате получается квадратурная формула

Y 1

= \u(t)H(x,t)dt-1 J

JT-v/l - X2 -1

п

= ^At(x)u(xt:) + Е(и,х),

—n

где коэффициенты

h

Л—

1 г2

I , J H{x,t)dt;

Tly\ -X h

(13)

(14)

Е(и, х) - остаточный член формулы (13).

Очевидно, что все Л ¿(х) неотрицательны для всех А' е [-1, 1]

я 1 1

%Ак(х) = —г=\Н(х,0Л = 1. (15)

-п Пл] 1 - X -1

Вычислив интеграл в равенстве (14), находим

Ак(х) =

1

г

xk+-~x LV 1 J

\ f H

, л

f

Xk ~ x v Z J

Л f

H

XlXk + r.

v 1 J

, Л

x> XA- ,, v 1 J Л

+

(16)

1

+ — 71

( ЬЛ ( ЬЛ

arcsin хк+ — - arcsin хк —

2 2

- v v -

Если функция u(x) непрерывна на отрезке [-1, 1], то остаток приближенной формулы (15) можно оценить неравенством

\u(x) — ii(x)\<—h, М= max Iи х I. (17) 1 1 2 I

Остаточный член квадратурной формулы (13)

1 г

Е(и;х) =—___ n(t)-ii(t) H(x,t)dt.

7t\l - "2

х2

Оценивая его по модулю с учетом (11), получаем неравенство

Е и;х < max |м(х)-м(х)|.

(18)

Подставив вместо интеграла в уравнение (10) его представление квадратурной формулой (13), приходим к равенству

и(х) - q(x)^Ak(x)u(xk) — F(x)q(x)E(u,x). (19)

-и

Удовлетворяя это соотношение в узлах квадратурной формулы (13), получаем систему равенств

п

u(Xj) - q(Xj (х, )и{хк) - F(Xj) + q{x} )Е{и, х;),

j = -п, ..., -1, 0,1, ..., п.

(20)

Удовлетворяя этому уравнению в узлах квадратурной формулы хк и отбрасывая слагаемое q(x)E(u, х), будем иметь систему линейных алгебраических уравнений:

uj-q(xj)Y,Mxj)uk =F(xj)

—п

j = -n, ..., -1,0,1, ..., п,

(21)

где uk - приближенные значения u(xk) в узлах xk.

Приближенное решение интегрального уравнения (10) можно получить из (19), если отбросить в нем слагаемое q(x)E(u, x) и заменить значения u(xk) приближенными значениями uk, найденными из системы (21). Имеем

п

й(х) = F{x) + q(x)^Ak (х)ик. (22)

—п

С другой стороны, точное решение

п

и(х) = F(x) + q(x)^Ak (х)и(хк ) + q(x)E(u, х).

—п

Сравнивая два последних равенства, получим

п

I и(х) - й(х)\ < \q(x)\Y}Ak (х)| | и(хк) - ик | +

—п

+\q(x)\\E(u,x)\.

Из этого неравенства, равенства (15) следует равномерная по хе [—1; 1] оценка погрешности приближенного решения (22)

\и(х)-й(х)\<-^—Щи, К). (23)

1 -q

В качестве приближенного решения уравнения (7) можно взять

Г(х) = B(x)F(x) - В(х)й(х)

Наука итехника, № 6, 2013

п

2

2

или, вспоминая, что q(x) = ставляя u (х) из (22), имеем

B( x)

. п

Г X = J\-X^Ak(X)uk.

и под-

(24)

Нетрудно убедиться, что точное решение можно записать в виде

Г(х) = ф-х2^Ак {х)и{хк) + ф - X2 E(и.х). (25)

-и

Из равенств (24) и (25) находим разность

Г(х) - Г(х) = <J\-x2 X А (х) и(хк) ~ ик +

—п

х2 E(u, х). (26)

Отсюда с учетом равенства (11) получим оценку погрешности Г(х).

Численное моделирование. В качестве примера решения были рассмотрены тело с отверстием единичного радиуса Я = 0,1 м и цилиндр радиусом г = 0,099 м. Материалом для тела и цилиндра являлась сталь со следующими физическими характеристиками: модуль Юнга Е = 2-1011 Па, коэффициент Пуассона V = 0,23.

Изложенное выше решение сравнивали с решением, предложенным И. Штаерманом [7], а также с численным решением с использованием конечно-элементной программы ANSYS.

Распределение контактного давления р, полученное в системе ANSYS и на основе аналитического решения И. Штаермана с помощью вышеизложенного алгоритма, представлено на рис. 2.

25 -, р, МПа

15

20 ф, град. 2В

Рис. 2. Распределение контактного давления в области соприкосновения цилиндра, вложенного в цилиндрическую полость: 1 - аналитическое решение; 2 - решение в системе ANSYS

Таким образом, на основе численно-аналитического метода решения уравнения Прандтля построено решение задачи о внутреннем контакте цилиндрических тел.

В Ы В О Д Ы

1. Сравнение результатов, полученных на основе построения приближенного решения (1) и с использованием конечно-элементного моделирования, показывает, что результаты близки. В то же время использование систем компьютерной математики (Mathematica, Mapple, MathCad и др.) для решения такого рода задач значительно проще и не требует изучения таких громоздких систем, как ANSYS, и дает результаты, точность которых не хуже, чем при решении их с помощью данного конечно-элементного пакета.

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

Л И Т Е Р А Т У Р А

1. Горшков, А. Г. Теория упругости и пластичности / А. Г. Горшков, Э. И. Старовойтов, Д. В. Тарлаковский. -М.: ФИЗМАТЛИТ, 2002. - 416 с.

2. Кравчук, А. С. Механика контактного взаимодействия тел с круговыми границами / А. С. Кравчук, А. В. Чигарев. - Минск: Технопринт, 2000. - 196 с.

3. Каландия, А. И. Математические методы двумерной теории упругости / А. И. Каландия. - М.: Наука, 1973. - 304 с.

4. Prandtl, L. // Nach. Kgl. Ges. WissMath. Phys. -1918. - Р. 451-470.

5. Голубев, В. В. Лекции по теории крыла / В. В. Голубев. - М.; Л.: ГИТТЛ, 1949. - 480 с.

6. Габдулхаев, Б. Г. Прямые методы решения сингулярных интегральных уравнений первого рода / Б. Г. Габдулхаев. - Казань: Изд-во Казанского ун-та, 1994. - 288 с.

7. Штаерман, И. Я. Контактная задача теории упругости / И. Я. Штаерман. - Л.: Гостехиздат, 1949. - 270 с.

Поступила 02.09.2013

-п

Наука итехника, № 6, 2013