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

Р и с. 3: Кривая нагружения, связанная с линиями уровня

а

/ / РнН

- / / &

/ /2 4 , 6 8 10 Я

он

-6

Р и с. 4. Сечение поверхности деформационных состояний плоскостью симмет -рии, проходящей через ось С и п

Это позволяет постулировать: при пересечении кривой, соответствующей процессу деформирования критической линии (6.1) происходит разрушение материала.

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

Если предположить, что механические свойства материала соответствуют поверхности нагружения п.5, и второй и третий инварианты тензора деформации Альманси мало влияют на разрушение материала, то уравнения критической линии (6.1) можно записать в виде

Г Е + Е2 + Е3 = И (Б);

[(1 -2Е)(1 -2Е2)(1 -2Ез) = 1. ( . )

Критическая линия в виде (6.2) совпадает с линией максимально возможного упрочнения материала, которая определяется величиной

ИШах = И (Б0шах), где — диссипация энергии при ортогональном процессе деформирования. Естественно предположить, что дополнительная диссипация энергии, производимая элементом объема материала при неортогональном деформировании, снижает способность материала к упрочнению. Это влияние может зависеть от уровня деформаций (величины И). Уравнения (6.1), (6.2) означают, что критическая линия уровня, определяющая момент раз -рушения каждой частицы, приближается к неде-формированному состоянию в процессе пластического деформирования соответственно диссипации

энергии. Функция Н(В) должна определяться экспериментально.

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

1. Быковцев Г. И., Ивлев Д. Д. Теория пластичности. — Владивосток: Дальнаука, 1998. —529 с.

2. Кочеров Е. П., Хромов А. И. Деформационные состояния и разрушения идеальных жесткопластических тел // Вестн. Сам. гос. техн. ун-та. Сер.: «Физ .-мат. науки», 2006. — № 42. — С. 66-72.

3. Хромов А. И. Разрушение жесткопластических тел, константы разрушения // Известия РАН. МТТ, 2005. — № 3.— С. 137-152.

Поступила 5.07.2006 г.

УДК 539.3

В. П. Федотов, Л. Ф. Спевак

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

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

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

параллельных вычислении перед последовательными и аналитических расчетов перед численными.

1. Алгоритм решения упругой задачи. Задача теории упругости для плоскоИ области Б состоит в нахождении вектора перемещении иі, тензора деформации Еу и тензора напряжении

Оу, которые удовлетворяют в этоИ области системе уравнении

Оі і=ьу; (1)

%=2 и і+иі, у); (2)

Оу = 2теу + 1^% 8у (3)

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

на поверхности 5/ : СцЩ + °12п2 = /1 = /і*, О21П1 + О22П2 = /2 = /2*>

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

Здесь Ьі — известные объемные силы (здесь и далее по повторяющемуся индексу производится суммирование от 1 до 2); иі у = Ци.; т — модуль упругости при сдвиге; V — коэффициент

Пуассона; 8у — единичным тензор; пі — вектор нормали к поверхности. Обратная к (3) зависимость будет иметь вид:

£» =

1 2т

1

3п £

О1- 7ГП°И 8Ї

, , (5)

1 + у 1 /

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

и (^)= | [и] (х))(х) - х )(х)] ^ (х)+ \\и*] (х ))(х) - /(Кх )и*,-(х)] ^(х). (6)

Здесь х е 5 — граничная точка области; X — внутренняя точка области; звездочкой (*) обозначены известные из граничных условий для группы поверхностей Би перемещения и* (х) и для группы поверхностей Б/ поверхностные напряжения /](х). Функции влияния и* (X,х) и /■] (X, х) для двумерной задачи имеют вид

и](Х>х) = -8Р(1-П>т[(3-4п)1п(г)8« -А(гА/]; (7)

/](х) = - 4Р(1-П)Г{[(1 -2у)8« + 2А'гА/Ё- (1 - 2п))А.гп] -А]гП)}, (8)

где г = г (х, X) — расстояние между точками х и X; А(г = -|х- .

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

вид

2Щ(хо) =| [и](хо,х)/](х) - (хо , х)и](х)] ^(х) +

+ | [и](хо ,х)/](х) - //(хо , х)и*](х)] ^(х). (9)

Это уравнение используется для определения неизвестных поверхностных перемещений и] (х)

и поверхностных напряжений /■ (х). Разбив границу области на элементы, мы приходим

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

2uf} = f J u* (x(P),x)ds(x)-Xu(a) J f*(x(P),x)ds(x)

2 a=1 s(«) a=1 f)

sfa)

N+M

N+M

X //“) 1 и* ((Р)’*)(х)- X 1 //(*(Р^х)(х)’

а=N+1 $( а) а=N+1 $( а)

I = 1,2, р = 1...М + N . (10)

Здесь М и N — количество элементов на каждой из групп поверхностей. Соотношения (10)

представляют собой систему линейных алгебраических уравнений относительно и(а) и /7-а). Коэффициентами системы являются интегралы 1 и* (х(Ь),х)dS(х) и 1 /*(х(Ь),х)dS(х) от

s(a)

s(a)

Х2і •X •X /

О с / А /\Ф

L Хі

и,/ °и,/

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

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

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

Рассмотрим на плоскости отрезок АВ, где А (А1, А2) и В (В1, В2) — произвольные точки, и

точку влияния Х(Х1, Х2) (рис. 1). Действующие на

отрезке АВ перемещения и = ( и1, и2) и поверхностные напряжения / = (/1, /2) вызывают в точке X некоторое перемещение и(Х) = (и1 (X),и2 (X)). Осуществим преобразование координат, сохраняющее расстояния, и отображающее точку А в начало координат О(0,0), а точку

В — в точку С(Ь,0), где Ь = ^(В1 - А1 )2 +(2 - А2 )2 — длина отрезка АВ. Такое преобразование является комбинацией параллельного переноса и поворота на угол ф (см. рис. 1). Произвольная точка х (х1, х2) на плоскости отображается в точку х (х1, х 2), связанную с ней соотношениями

х = Кх + А; х = К-1 (х - А), (11)

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

где x =

/ \ /— \ A, ї

x1 — x1 , A =

, x = 12 AA

x2 x 2 V /

22

, K — матрица поворота: K-1 = Kт cos j, k12 = -k21

rk k Л

*11 *12

k21 k22

B2 - a2

-sin j.

(12)

І ' ^ I

Здесь и далее во всех соотношениях одномерные массивы трактуются как матрицы-столбцы.

Произвольный вектор ^ = (w1,w2) на плоскости отображается в вектор ^ = (W1, W2):

w =Kw , w = К~lw. (13)

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

ли на отрезке ОС действуют перемещения и = К 1и и поверхностные напряжения / = К 1 /, то они вызовут в точке Х = К_1 (X- А) перемещения

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

и (Х) = К 1и (^)-I (и** )= 1 и* (X, х )<® (х); I (/■*) = | /*(Х, х )dS (х);

АВ АВ

1 (и * )= 1и* (^ х )dS (х); 1 (/*)= 1 /* (^ х )dS (х )•

(14)

(15)

ОС

ОС

Используя соотношение (14) и выражение (6) для перемещений точек X и X, можно установить следующую связь между интегралами от компонентов функций и (X, х) и / (X, х) вдоль отрезка АВ и интегралами от компонентов функций и (X, х) и / (X, х) вдоль отрезка ОС :

I (и*! ) = ГI(и*! ) + 2Г21(и*2 ) + ГзI(^ ) ; I (и*2 ) = I («2х) = >2 (7(^2 ) -1(и*1)) + ( - > )(и*2 ) ;

I (и*2 ) = Гз I (и*1)-2Г21 (и*2 ) + ГI (и*2 ) ; I ( /11 ) = Г11 ( /11) + Г2 (I ( /1*2 )+1 ( /2*1)) + Г31 ( /2*2 ) ;

где

I (/1*2 ) = >2 (I (/2*2)-1 (/1*1)) + >11 (/1*2)-гз I (/2*1);

I (/2*1 ) = >2 (I (/2*2)-1 (/1*1))-гз I (/1*2) + >11 (/2*1);

I ( /2*2 ) = Гз I (и*1)-Г2 (I ( /1*2 )+1 ( /2*1)) + Г11 (и*2 ) ,

Г1 = кп = ^22 = к11к22 = 008 ф ;

г2 = к11к12 = -к11к21 = к12к22 = -к21к22 = -^28т2ф ;

(16)

(17)

гз = к12к21 = к12 = к21 = 8*п ф .

Таким образом, мы получили, что для вычисления интегралов по произвольному отрезку АВ для точки влияния X достаточно построить матрицу К (12), определить через соотношение (11) точку X и вычислить интегралы по отрезку ОС, что является более простой задачей. Для вычисления интегралов I (и*) и I(/*) для произвольной длины отрезка Ь и произвольной точки влияния X получены простые аналитические формулы, которые приведены ниже:

I (и*1) = С1 (2 + 1) ( - Ь) + ^ ( + XlQ2 )

V 2

(и12 ) = -2С1 ^^2 ; 1(и22 ) = С1 (с2 - 1)^бз + С2 ^((1 +X1Q2 )- Ь

V

I(/1*1) = сз (( + 1)бз +X2dl) ,1 (/1*2 ) = С

(18)

I (/*1)

= сз

2 С402 X2d 2

Г(/22) = Сз ((с4 + 1)0з )>

где

с =-

с2 = з - 4у, сз = —

с4 = 1 - 2у

8яц(1 -V) ’ 4я(1 -V)

01 = 1П(( - Ь)2 + § ], б2 = 1п (X2 +X2 )- 1п ((( - Ь)2 +X

^1 - ЬЛ

бз = агс*8

- аг^

X2

d- X1-L _Л]_ d -X

d1 _ ч 7 -2 -2 -2 ’ d2 _ ^2

1

1

tt2 tt2 /— \2 —;

X1 + X2 ( - L) +X

2 tt2

2

(-l)+X2 X + X2 ^ ^

Приведенные формулы не могут быть использованы в случае, если Х2 - 0, то есть, точка влияния X лежит на прямой, содержащей отрезок AB. Для этого случая получены следующие специальные формулы:

1 (и*)- Cl (С2 (LQ* -l) + Xi02*)-L),

1(u*2) - 0, 1 («2*2) - qC2 (L (q* -1) + XiQ2),

1 (/1*1)-1 (/2*2)-0, 1 (/1*2 )--1 (/2*1 )--C3C402*. (20)

где

6* - ln| L-X1, 62 - In I- ln IL-X1- (21)

Что касается неопределенностей при условиях Х2 - 0 , X1 - 0 или Х2 - 0 , X1 - L , соответствующих совпадению точки влияния X с одним их концов отрезка AB, то для выбранной ап-

проксимации они не реализуются, поскольку узловые точки х(b) находятся в серединах элементов.

3. Аналитическое вычисление интегралов по дуге окружности. Рассмотрим произвольно ориентированную дугу окружности M , которая определяется следующими параметрами: MC (х1с, х2с) — центр;

R — радиус; ja и jb (ja <jb ) — полярные углы, соответствующие началу и концу дуги при обходе против часовой стрелки, соответственно (рис. 2).

Учитывая связь декартовых координат х1, х2 с координатами р, j полярной системы с центром MC

х1 - х1с +р cos j, х2 - х2 C +р sin j, (22)

можно свести интегралы от компонентов функций влияния по дуге M к определенным интегралам:

jb

J «* (X, х )ds (х )- J и* ( ^2, х2 )ds (х )- R J U* ( х1с + R cos j, х2 с + R sin jJ j ,

Р и с. 2. Произвольно ориентированная дуга окружности

M

M

jb

J /* (Х х)dS (х) - J /* (Х2,хь%2 )ds (х) - R J /* (Х2,х1с + Rcos j,х2с + R sin j)) j . (23)

MM ja

Отметим, что интеграл по дуге окружности от функции ln (r), входящей как слагаемое в диагональные компоненты функции и* , не вычисляется аналитически, поэтому для вычисления интеграла от этого слагаемого использовалось разложение в ряд Тейлора по переменной j

j a +jb m

в окрестности точки j0 - 2 соответствующей середине дуги. Тесты показали, что при

величине дуги меньше 10° достаточно учесть отрезок ряда четвертой степени. В этом случае можно принять, что

j b г I-----

: Jln |v(

Нс '

R cos j) +(Х2

^2 с '

R sin j)

- ln (Q10)

0\(jb - j a ) „0 (jb - ja )

- 62

q0

1 - 9Q3 - 6

60

160

(24)

где

ві = ^ІП

Фа +ФбЛ

+ 52 008

в4 = 51 йІП

Фа +ФбЛ

+ 52 008

Фа +ФбЛ

Фа + Фъ

+ 53, в20 =-

вів0 +(в5°)

24 (в°)

вз° =

в°л2 ві°

в° = ^0081 Фа +ФЪ І-

Фа +ФЪЛ

51 = 2дЯ2, 52 = 2рЯ2, 53 = (р2 + д2 +1)Я2; р = Хіс Х, д = Х2с Х . (25)

' ' Я Я

Интегралы от остальных составляющих функций влияния были вычислены аналитически. В

результате были получены следующие формулы:

| М11 (Хх)(х) = 2с1Я {2с21 ° — 2в1 + 85в2 — 8бв4 — й2 ) ,

М

| М22 (Хх)й$ (х) = 2с1Я {2с21 ° — 2в1 — 85в2 + 8бв4 + й2 ) ’

М

| МІ2 (X, х)<К(х) = С1Я (— 48бв2 — §5в4 + ) ,

М

| /1І (X,х)<&(х) = С3 ((4 + 1в + ^ ) , | /122 (X,х)<&(х) = с

где

М

М

" с4в4 + 4йб Л 2 84

| /21 (х)&(х)

= с

М

с4в4 + 4й б ' 2 84

| /22 (X,х)<К(х) = С3 ((С4 + 1)в3 — ^ ) , (2б)

М

в1 = ar0tg (/ъ )— ar0tg (іа), в2 = в1 + ar0tg

^ 84/ъ + 2дЛ — arotg "84/а + 2д '

83 V 83 \ з /

в3 = в1 — ar0tg

ґ 84/ъ + 2дЛ + arotg "84/а + 2д '

83 V 83 \ з /

в4 = 1п|й4ъ| — ІП1й4а I — ІП (й1Ъ ) + ІП (й1а ) ; (27)

й1а = 1 + /о й1Ъ = 1 + Ь , й2 = —

2

81

ґд + р/ъ д + ріс Л

й

1Ъ

й,

1а

й3 =~

2

81

р — 9/ъ — р — діа

1Ъ

й

1а

й4а = 84^ + 4діа + 82 , й4ъ = 84^ + 4д/ъ + 82 ,

й-5 — —

5 84

1

й4Ъ \ V

89{Ъ + 81°/Ъ

й

—д88

1Ъ

1

йб =-

4Ъ

84/2Л

д82/ъ—д2 р + "4‘Ъ

1Ъ

89/а + 81°/а — д8;

1а

4а

ґ 2 2 8 /2Л д&А — д р +

//

Ма

й-

1а

(28)

81 = р2 + д2, 82 = (р +1)2 + д2, 83 = р2 + д2 — 1, 84 = (р — 1) + д

8 = д4 — р4 + р2 — д2 8 = (р + д 1)рд

85 = ^ 9 8б = "

„2 2

/а = tg

2 2 2 > 7

(р + д )

д2), 51° = (р —1)( р

"Фа ' 1 2 \ , /ъ = ^ Фъ " 1 2 0

(29)

(3°)

Очевидно, что выражения (30) предполагают, что фа ^ я и ф6 ^ я. В этих случаях соотношения (27), (28) следует принять в следующем виде: при фа = Я —

84?ь + 2д

в1 = Ръ ) — - , в2 = в1 + ar0tg

— (83)!

03 = 61 - аг^

84/Ь + 2д

83

3

+ (83), £4 = 1п |^4^ - 1п ^1Ь )- 1п (84);

(27а)

^ь = 1 + /2, d2 = 2 (! +РЬ ) , $3 = 2 ( р ^ ) , $4Ь = 84^! + 4^/ь +82,

81$1

1Ь

8^

$5 —

84$ 4Ь

89/3 +810?Ь - д88

$

1Ь

$6 = -

13

43

д^ь - д 2 Р + 8тЬ-

при Фь = Р

£1 = р - (а ) , £2 = 01 + Р «ип ( 83 ) -

-*1Ь

^а + 2?Л

(28а)

Р

03 = 01 - (83) + агс18

84?а + 2<?

83

V ^ /

$1а = 1 + 1а , $2 = -

2 ( 4 + Р^ ) $3 = -

04 = 1п (84 )- 1п| $ 4а| + 1п ($1а ) ;

2 ( Р - ?^а )

(27Ь)

8^

$5 —

1а

89?а + 810^а

81$

1а

84$,

4 4а

- ?88

1а

$6 = -

1

$4а = 84^а + 4д/а + 82 ,

2

2, „2„, 84'а

4а

д8^ - д р +

$

1а

(28Ь)

Формулы (24)-(30) не могут быть использованы в случаях, когда: 1) точка X лежит на окружности, содержащей дугу М (83 = р2 + д2 -1 = 0); 2) когда точка X лежит в центре этой окруж-ности (81 = р2 + д2 = 0). Для этих случаев получены специальные формулы.

В случае 1) интеграл /0 может быть вычислен по формулам

10 = 1п ((й*) (Фь - Фа ) + —£

1600

(Фь - Фа )3 , Г1 1 ] (ФЬ - Фа )5

3 V _ Л V 3 00 0 80 /

если ф^Фа +Фь , (31)

/ 0 =

1п (Я) + 1п Фь Ф -1

(Фь-Фа -(-У!, если фх = фаФ,. (31’)

288 230400

Здесь

00 = 1 - соэ

(32)

ФХ - полярный угол, соответствующий в системе координат (22) лежащей на окружности точке X. Итоговые формулы для случая 1) имеют следующий вид:

| и*1 (X, * )$£ (х) = с1Я

М

С2/0 + 01* - ФЬ 2 Фа |, |и12 (X, Х)$Б (х) = -с1Я02 =

М

М

| и22 (X, х)<® (х) = С1Я |с2/0 - 01* - ФЬ Ф а I, | /1*1 (X, х)$8 (х) = с3 (с4 + 1) Фь 2 Фа -

2 0 М V 2

| /1*2 (X, х(х) = с3 (с403 + д*2 ) , | /2*1 (X, х)&(х) = с3 (-с403* + 02 ) =

М

М

| /2*2 ( х№(х) = С3 (С4 + 1) ФЬ ■ Фа + 02

(33)

где

М

0» = (фь -Фx)- эт (фа -Фx ) 0* = СОв (фь -Фx)- с^ (фа - Фx )

01 = ;; , 02 =

2

0* = 1п

2

Фь -Фx

- п

Фа -ФX

(34)

В случае 2) необходимости приближенного интегрирования нет, поскольку для точки X, лежащей в центре окружности 1п (г) = 1п (Я). Формулы для вычисления интегралов имеют вид:

м

м

1

| м22 (X, х )<& (х ) = С1Л с2 1п (Л )-2 1(фъ-Ра ) + 01*

м \Ч 0

| /1*1 (х)8(х) = С3 (( + 1)(ф6 ф ) + 2°Г ) , | /1*2 ( х)&(х) = -2сз02*

м м

| /*1 ( х)&(х) = -2сз02*, | /22 (х)<&(х) = с3 ((С4 + ^(фъ -Ра )- 201* ) ,

м

м

где

о111 =

®^п (2Ръ )-81п ( 2Ра )

о* =-

сов (2Ръ )- сов (2Ра )

(36)

4 4

В заключение отметим еще раз, что фа < фъ. В случае, когда фа >Ръ, т. е., если обход вдоль окружности происходит по часовой стрелке (например, исследуемая область имеет внутреннее отверстие), следует вычислять интегралы по приведенным формулам, сменив знак у интегралов

от функций /у на противоположный.

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

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

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

2. Думшева Е. С., Зенкова Е. С., Федотов В. П. и др. Численно-аналитический алгоритм для решения задач упругости, теплопроводности, диффузии // Алгоритмы и программные средства параллельных вычислений. —

Екатеринбург: УрО РАН, 2003. — Вып. 7. — С. 70-86.

3. Федотов В. П., Спевак Л. Ф., Трухин В. Б. и др. Исследование сходимости численно-аналитического метода решения задач упругости, теплопроводности и диффузии // Вестн. Сам. гос. техн. ун-та. Сер.: «Физ.-мат. науки», 2004. — № 30. — С. 24-32.

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

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

Поступила 31.03.2006 г.