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

УДК 539.3

В. П. Федотов, Л. Ф. Спевак, В. Б. Трухин

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

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

Для решения статических задач теории упругости авторами был предложен численно-аналитический метод решения [1, 2], основанный на методе граничных элементов [3, 4]. Метод позволяет построить алгоритмы решения, которые сочетают в себе преимущества параллельных вычислений перед последовательными и аналитических расчётов перед численными. В работе [5] был представлен новый подход к получению аналитических формул для вычисления интегралов от компонентов функций влияния, позволяющий существенно ускорить процесс решения. Данная работа посвящена аналитическому интегрированию производных от функций влияния, необходимому для вычисления напряжений внутри деформируемой области.

Статическая задача теории упругости. Задача теории упругости для плоской области Б состоит в нахождении вектора перемещений щ, тензора деформаций £{у и тензора напряжений , которые удовлетворяют в этой области системе уравнений:

= Ь1, (1)

= - [иц + щ,у), (2)

£Ч = 2 ищ’у 2u,v

а{у = 2^£гу + --— £кк 8{у, (3)

1 — 2v

и заданным граничным условиям

на поверхности Sf : апп1 + ^12^2 = /1 = /* , ^21«1 + ^22«2 = /2 = /2*, (4)

на поверхности Su : и1 = и* , и2 = и**. ( )

Здесь Ь{ — известные объёмные силы (по повторяющемуся индексу производится суммирование от 1 до 2); щ,у = дх1, Л — модуль упругости при сдвиге; V — коэффициент Пуассона; б(у — единичный тензор; П{ — вектор нормали к поверхности. Обратная к (3) зависимость имеет следующий вид:

1 ( ^ \

£"' = — аа ■ (5) Согласно методу граничных элементов [3, 4] решение системы (1)-(4) можно выразить через поверхностные перемещения щ (х) и поверхностные напряжения /• (х):

щ ({) = ^ и *■($, х)/* (х) — /* х) и у (х) dS (х) + ^ и * у ($, х) /у (х) — /* ($, х) и* (х) dS (х). (6)

Sf ^

Здесь х е S — граничная точка области; $ — внутренняя точка области, звёздочкой ( *) обозначены известные из граничных условий для группы поверхностей Su перемещения и* (х) и для группы поверхностей Sf — поверхностные напряжения /* (х). Функции влияния и*у ($, х) и /* ($, х) для двумерной задачи имеют вид:

и*у ($, х) = ——-1—— [(3 — 4v)ln(r)81у — Дг-гДуг], (7)

4 8п (1 — V) л

1 ( дг )

/*у($’х) = — 4л(l-v)r\^(1 — 2V)б{у+2Д{ГДу^ дп — (1 — 2V)(Д{ГПу — Дугщп , (8)

где г = г (х, $) — расстояние между точками х и $, Дг-г = ■

Для произвольной точки хо гладкой границы граничное интегральное уравнение имеет такой вид:

1 щ (хо) = ^ и * у (хо, х)/* (х) — /*у (хо, х) и у (х) dS (х) + ^ и * у (хо, х) /у (х) — /*у (хо, х) и* (х) dS (х). (9)

Это уравнение используется для определения неизвестных поверхностных перемещений иу (х) и поверхностных напряжений /у (х). Разбив границу области на элементы, мы приходим к численной схеме, в основе которой лежит граничный элемент (отрезок прямой либо дуга окружности) и произвольная точка, на которую оказывают влияние напряжения и перемещения, действующие на этом элементе. Если сделать разбиение поверхности деформируемой области достаточно мелким, то все поверхностные перемещения и напряжения, включая неизвестные, можно считать постоянными на элементе и помещёнными в его середину. В этом приближении, записав уравнение (9) для середины каждого граничного элемента, получим следующую систему уравнений:

1 (в) — и.

2 1

N

= Е /<

а=1

(а) *

^ и*у{х(в), xjdS (х) — ^ир^ fiУу[x(в), xjdS (х)+

N+М

+ I /у

а=N +1

/( \ N+М п . .

и*:^, х^ (х) — £ и(а) * I /*:(х(в\ х^ (х), • = 1,2, в = 1,2,..., М + N. (10)

-1 ' ' а= N +1 1 3 у ' >

а=N +1

о( а) эи

Здесь М и N — количество элементов на каждой из групп поверхностей. Соотношения представляют собой систему линейных алгебраических уравнений относительно иа и /!а.

эффициентами системы являются интегралы / и*Ах(в\ x)dS(х) и / /*Ах(в), x)dS(х) от компо-

_(а) V ) ,а) V )

(10)

Ко-

з(а /

и, /

( ) "и,/

нентов функций и*у и /"*у по граничным элементам для точек влияния, помещённых в узловые

точки х(в. В работе [5] были получены формулы для вычисления таких интегралов по произвольному отрезку прямой и произвольной дуге окружности для произвольной точки влияния.

( а) ( а)

После решения системы и получения значений иу и /у перемещение в любой точке внутренней области вычисляется по формуле (6). Деформации определяются по формуле

цу (О) = ^ ш* ук (О, х) /к (х) — gу (О, х) ик (х)

dS (х),

полученной из (6) путём дифференцирования, с учётом (2). Здесь

ш*у к (о х)=2

ди *

ди * , ___к +______¿к

дОу дО

1

gг у к (О х) = 2

(д £

д /*, Iк + ■)ук

(11)

(12)

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

Аналитическое вычисление интегралов по отрезку прямой. Как и в работе [5], для упрощения процедуры вычисления необходимых интегралов поставим в соответствие расчётному блоку из произвольного граничного элемента и произвольной точки влияния блок из специального элемента и соответствующей точки влияния таким образом, чтобы интегралы для двух блоков были связаны.

Рассмотрим на плоскости отрезок АВ, где А (А1, А2) и В (В1, В2) — произвольные точки, и точку влияния О (О1, О2) (рис. 1). Действующие на отрезке АВ перемещения и = (щ, и2) и поверхностные напряжения / = (/ь /2) вызывают в точке О некоторое перемещение и (О) = (и1(О), и2(О)). Осуществим преобразование координат, сохраняющее расстояния, и отображающее точку А в начало координат О (о, о), а точку В — в точку С (Ь,о), где Ь = \/(В1 — А1)2 + (В2 — А2)2 — длина отрезка АВ. Такое преобразование является комбинацией параллельного переноса и поворота на угол ф (рис. 1).

S

Произвольная точка х(х^ х^ на плоскости отображается в точку х(х 1, хО, связанную с ней соотношениями:

х = Кх + А, х = К—1(х — А), (13)

где х = | х^ |, х = | х, |, А ^ ^ |, К — матрица поворота:

Х2 А

К = [ I11 112\, К-1= Кт, кц = к22 = В-А-

к21 к22 ’ ’ 11 22 Ь

к12 = -к21 = -

В2 - А2

Ь

= - 81П р.

= 008 р, (14)

О Ь х\

Рис. 1. Преобразование координат

Произвольный вектор ш = [ш1, ш2) на плоскости отображается в вектор ш = (ші, ш2):

ш = Кш, ш = К 1 ш.

(15)

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

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

'К,) = /

ЛБ (х), I

АВ

:Ы=/

дЪ;

I

-Л- -

ОС ОС

АВ

д{ V-’-” I д<

ЛБ (х), I

(¡’к.-)=!

АВ

(Л..!)=!

д¡к ((• х) КI д¡"к (Ъ х)

ІБ (х), ІБ (х).

(16)

Поскольку точки О и О связаны преобразованием координат (13), справедливы следующие равенства: _ _

ди*к = ди± д_О1+ д_ик дО2= к . д_ик + к . Щк

дОу дО1 дОу дО2 дОу 1 у дО1 2у дО2 ’ (17)

д/*к д/* дО1 + д/*к дО2 , д/к +, д/*к ( )

----= —=-------I =-------------= к1 ; —=—+ ко у —=—.

дОу дО1 дОу дО2 дОу дО1 дО2

Учитывая последнее равенство, а также полученную в работе [5] связь между интегралами по отрезку АВ для точки влияния О и соответствующими интегралами по отрезку ОС для точки О, можно получить следующие формулы для вычисления интегралов от производных функций и * ■ (О, х) и /■*■ (О, х) по произвольному отрезку АВ для произвольной точки влияния О:

у у

ди * (Ъ, х)

( ) г ди2 «, х) г ди.7 (

'К,) = кііІ -*=-------------------ІБ(х) + Ы) ІБ(х),

АВ ^ АВ Ъ2

( ) Г д¡2 (Ъ, х) г д¡2 (Ъ, х)

і(Пч) = к1 і) -ид=--------------¿Б(х) + к211 -------ЛБ(х)

АВ АВ ^2

(18)

где

АВ

АВ

АВ

АВ

АВ

АВ

ди{1 ((, х)

дЪ

дик2 (Ъ, х)

дЬі

ди*2 (Ъ, х)

дїі

ЛБ (х) = Г11 (и И,і) + 2г2~І[иІ2іі) + Г31 (и22,і) ,

ЛБ (х) = Г2 11[и*22, і) - _І(и п, і)) + (г1 - г3) ~І(и*и, і), ЛБ (х) = Г31 (и И,і) - 2г2 1 [и*22,і) + Г11 [и*22,і) ,

В

АВ

/

АВ

dS (х) = и 1 (л*1,0 + г2 (7 [/*„) + + Г31 (г*2Л),

^ (х) = Г2 (7 (/*2Л) -1 (/*1',)) + и1 (/*2Л) - Г31 (/*1Л),

дА\($, X)

д$г

д /*($, х)

д$г

I С^ (х) = Г2 (1 (/*2.) -1 (/*1.)) - ^1 (/*2.) + п1 (/*1'.) '

АВ $1

/

АВ

г 1 = к^ = к^2 = к11к22 = СО&2 р, г 2 = кцкц = -кцкц = к12к22 = -к21к22 = - 2$т2р, г з = -к12к21 = к12 = к^ = 8Ш2 р.

д/*2«, X) д$г

Ж (х) = гз ¡¡и *1,.) - г2 (/(/Ц,,-) + I |/>*1,г'11 + г11(и*22,г) ,

(19)

Для вычисления интегралов от производных функций влияния по отрезку ОС были получены следующие формулы:

Ци^д] = С^1 C2Q2 + ^2^,

1 и

11,2)

С1 ((С2 + 1) Qз + ¡,2а^,

*12,1) = С1$2а1,

где

1 К„)

1 (И 22,1) = С1 (2 С2 ^2 - $ 2а2^ ,

1 (Л1 11) = С3 {^С4а2 +

1 (/1 2,1) = С3 (С4 а1 +

1 (/2 1,1) = С3 {<-С4а1 +

1 (/22,1) = С3 ((с4 + 2) а2 - а3), 1 [/22,2) = С^С4а1 - а4), 11

^ и12,2) = С1 (-2 ^2 - $2 а2^>

11 и22>2) = С1 \[С2 - 1) Qз - ^1] , 1 (Ли) = С3{(С4+2) а1 +

1 (/\2,2) = С3{,- (С4 + 2) а2 + а5),

1 [Л*!. 2) = С3{(С4 - 2) а2+а5),

(20)

С1 =

8л^ (1 - V)

е2 = 3 - е3 = -

4п (1 - V)

С4 = 1 — IV,

Ql=ln(($1 -ь) + Q2=ln($\ + $2) -($1 -ь) + Qз = arctg

1±

¿2,

аг^

$1- ь

*. $ ,

. $1 - ь $1 . $

а1 = (----------а2= $2

(- \2 -2 -р2 -

($1 - ^ + $2 $1 + $2

а4 = 2$2

/ ч 1 1 , а3 = 2$2 г I2! ($1- ь)2

-2 -2 )2 -2 ¿1+ $2 ($1 - ь) + $2 ^ ($2+$2 )2 ( Ъ - ь)2+~$2)2

^ ’ 1 > > )

\ г \

$1

$1 - ь

($1+$2) |[$1-ь)2+$2|

, а5 = 2$2

($1+$2) |[$1-^2+$2|

(21)

Приведённые формулы не могут быть использованы в случае, если $2 = 0, то есть когда точка влияния $ лежит на прямой, содержащей отрезок АВ. Для этого случая получены следующие специальные соотношения:

1 (и2и) = с^* , 1 (и1 1,2) =0, ^{/*1,1) = 0,

1 (/22,1)=0,

1( /2*1,2)=0,

1(и 12,1) = 0,

1(и *2,2) = -С^2,

1( /22,1) = ^3=4^

1_(и*2А) = С1С20* , 1(И2 2,2) = °’

7^ 1 ) _ _ С3С4ь

1[/21'1>= $1($1-ь)

$1($1-ь)

1(/112 = ыЦ ’ 1/22=0,

Т ( f 1 ) = С3С4ь

1[122’2> = $1-ь) •

1

1

2

2

1011 Па, V = 0,33.

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

Задача об однородном нагружении квадратной пластины с круговым отверстием. Рассмотрим упругую квадратную пластину с круговым отверстием в центре, подвергнутую однородному растяжению (рис. 2). Длина стороны пластины а = 20 м, радиус отверстия Я = 1 м. Параметры упругости были приняты следующими: Е = 2 •

Предполагаем, что деформация протекает в условиях плоского напряжённого состояния, поэтому во всех приведённых выше соотношениях, соответствующих плоскому деформированному состоянию, нужно заменить коэффициент Пуассона V на величину V' = 1Vv [3]. Для бесконечной пластины, в соответствии с решением, приведённым в [6], окружное напряжение вдоль оси х1 равно

/ ( Я2 Я4

°в =2 Г+ Т2 +374

(23)

где / — однородное напряжение на границе; г — расстояние до начала координат.

Ввиду симметрии задачи нет необходимости рассматривать всю деформируемую область. Сначала была рассмотрена четверть пластины (рис. 2 а). В этом случае граничные условия выглядят следующим образом:

х1 = 0: щ = 0, /2 = 0;

Х2 = 0 : Ы2 = 0, /1 = 0,

Рис. 2. Квадратная пластина с круговым отверстием: а — модель для четверти пластины; б — модель для половины пластины

хі = 2 Х2 = 2

/1 = 0, /1=/,

/2 = 0; /2 = 0.

(24)

Полученные решения системы линейных уравнений вида (10), а именно, значения (-/2) на границе Х2 = 0, сравнивались с аналитическим решением (23). Вычисления проводились для 200 элементов, при этом в одном решении круговая часть границы аппроксимировалась ломаной, и использовались только формулы для вычисления интегралов по отрезку, а во втором круговая часть разбивалась на дуги, для которых использовались соответствующие формулы интегрирования. Результаты расчётов при / =1 Н/м2 показаны на рис. 3.

Далее была решена задача для половины пластины (рис. 2 б). Граничные условия в этом случае выглядят следующим образом:

1 2345678 хъ м

Рис. 3. Модель для четверти пластины. Сравнение численных и аналитических решений: сплошная линия — аналитическое решение; точки — четверть пластины, аппроксимация отрезками; пунктирная линия — четверть пластины, аппроксимация дугами

Х1 = 0: ^1 = 0, /2 = 0; Х1 = |: /1=0, /2 = 0;

2 : /1 = -/, /2 = 0; Х2 = § : /1= /, /2 = 0.

Х2 = -

(25)

Для этой модели сравнивалось напряжение 022, вычисленное для внутренних точек вдоль оси Х1 в соответствии с соотношениями (6), (11) и (12), с аналитическим решением (23). Результаты сравнения показаны на рис. 4. Расчёты вновь проводились для 200 элементов с аппроксимацией круговой части границы отрезками и дугами.

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

Чтобы показать различия всех численных решений, которые были получены для одного и того же числа граничных элементов (200), сравнилось их относительные отклонения от аналитического решения, (рис. 5). Сравнение показало, что в модели для четверти пластины аппроксимация круговой части границы ломаной линией дает худшее решение по сравнению с полученным при аппроксимации дугами, особенно вблизи границы области. Это говорит о том, что система уравнений (10) более корректно составляется именно при аппроксимации дугами. Два решения, полученные для половины пластины, практически не отличаются, что показывает устойчивость предложенной процедуры вычисления напряжений внутри области к погрешностям, допущенным при аппроксимации границы.

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

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Гасилов, В. Л. Численное моделирование упругой задачи на многопроцессорных вычислительных системах [Текст] / В. Л. Гасилов, Е.С. Думшева, Е.С. Зенкова, В. П. Федотов / В сб.: Алгоритмы и программные средства параллельных вычислений. — Екатеринбург: УрО РАН, 2002. — Вып. 6. — С. 104-124.

2. Думшева, Е.С. Численно-аналитический алгоритм для решения задач упругости, теплопроводности, диффузии [Текст] / Е. С. Думшева, Е.С. Зенкова, В.П. Федотов и др./ В сб.: Алгоритмы и программные средства параллельных вычислений. — Екатеринбург: УрО РАН, 2003. — Вып. 7. — С. 70-86.

3. Бреббия, К. Методы граничных элементов [Текст] / К. Бреббия, Ж. Теллес, Л.М. Вроубел. — М.: Мир, 1987.— 524 с.

4. Бенерджи, П. Методы граничных элементов [Текст] / П. Бенерджи, Р. Баттерфилд. — М.: Мир, 1984. — 494 с.

5. Федотов, В. П. К аналитическому вычислению интегралов в численно-аналитическом методе решения задач математической физики [Текст] / В. П. Федотов, Л.Ф. Спевак // Вестн. Сам. госуд. техн. ун-та. Сер.: Физ.-мат. науки. — 2006. — № 43. — С. 92-99.

6. Тимошенко, С. П. Теория упругости [Текст] / С. П. Тимошенко, Дж. Гудьер. —М.: Наука, 1975. —576 с.

Институт машиноведения Уральского отделения РАН, г. Екатеринбург Поступила 29.03.2007

fedotov@imach.игап.ги

Рис. 4. Модель для половины пластины. Сравнение численных и аналитических решений: сплошная линия — аналитическое решение; точки — половина пластины, аппроксимация отрезками; пунктирная линия — половина пластины, аппроксимация дугами

1 2345678 II, к

Рис. 5. Относительное отклонение Д численных решений от аналитического: 1—четверть пластины, аппроксимация отрезками; 2— четверть пластины, аппроксимация дугами; 3— половина пластины, аппроксимация отрезками; 4— половина пластины, аппроксимация дугами