Численно-аналитический метод решения задач упругости с особенностями Текст научной статьи по специальности «Физика»

Механика деформируемого твердого тела

УДК 531/534

В. П. Федотов, А. В. Горшков

ЧИСЛЕННО-АНАЛИТИЧЕСКИЙ МЕТОД РЕШЕНИЯ ЗАДАЧ УПРУГОСТИ С ОСОБЕННОСТЯМИ

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

При решении практических задач математической физики (упругости, теплопроводности, диффузии, пластичности и т.д.), как правило, используются алгоритмы и программистская технология последовательного счёта. В последнее время наблюдается прогресс в распараллеливании программ для решения таких задач, однако возможности распараллеливания в классических пакетах программ остаются ограниченными. Разработка алгоритмов, в которых изначально была бы заложена идеология распараллеливания, может существенно сократить время решения реальных задач. Вторым резервом для роста быстродействия может быть увеличение аналитической части предварительных расчётов в рассматриваемой задаче. В работе [1] была рассмотрена плоская задача теории упругости и предложен подход, основанный на методе граничных элементов, позволяющий ввести распараллеливание на уровне алгоритма, и существенную долю расчетов провести аналитически. В целях упрощения вычислений в ней принималось, что как известные, так и искомые значения поверхностных напряжений и перемещений были постоянными в рамках каждого граничного элемента. Этого допущения оказалось вполне достаточно для задач с регулярной границей, что показала быстрая сходимость численного решения к аналитическому для тестовых задач [1]. Однако для задач с особенностями вблизи границы наблюдалась некоторая осцилляция по напряжениям, что объясняется принятым допущением. В настоящей работе этот недостаток был устранен.

Как известно, метод граничных элементов сводит решение упругой задачи для тела V с границей Г к решению интегрального уравнения [2]:

и, (Х) = | Ы1 (Х Х) РУ (Х) аГх - | Р* (Х Х) 4у (Х) аГх + | иI (Х Х) Ьу (Х) ^х • (1)

Здесь х = (х1, х2) — координаты точки интегрирования, X = (X, X) — координаты точки «наблюдения», ц. (х) — . -тая компонента перемещения точки тела, р1 (х) — . -тая компонента усилий, иу (X,х) — компоненты перемещений, а Ру (X,х) — компоненты усилий фундаментального решения. Для плоской задачи функции и*. (X, х) и р*. (X, х) имеют вид

*. (X, X )= -(3 - П)1п г 8.у +

дг дг дх, дх,

/(лЕ (1 -п)).

Р, =■

((

(

\\

2

дг дг дх. дх.

дг

дп

(1 - 2п

дг

дг

-------п.---------------п2

ч дх2 дх. у

/(4л (1 -п)г) .

дг

Здесь г2 =(1 - х1) +(2 - х2) ,-производная по нормали к границе; V — коэффициент

дп

Пуассона; Е — модуль упругости.

V.

На части границы Ц заданы перемещения u. (x), на части Г2 — усилия. Согласно методу граничных элементов выберем на границе области N узлов с координатами X” = (x”5X2), n = 1, 2, ..., N, которые разобьют границу на N элементов Г”. В данном случае в качестве граничных были взяты отрезки прямых. Элементы и узлы нумеруются последовательно, против часовой стрелки. При этом область получается слева, а внешняя часть — справа. Каждому эле-

¿r1n ¿Т 2”

менту ставится в соответствие два узла — начало элемента точка x , и его конец точка x —

в направлении обхода нумерации. Первая цифра индекса показывает, что эти точки относятся к элементу с номером n .

На элементе перемещения аппроксимируются линейными функциями, усилия — кусочнопостоянными:

un (s) = uJX (s) + u2 j (s),

p” (s) = p”.

Здесь p” — значение i -той компоненты усилия на элементе с номером n , uk — значение i -той компоненты перемещения в узле k элемента с номером n , jkn (s) — базисные функции элемента n , j (s) = 1 - s / Ln, qj (s) = s / Ln, Ln - длина элемента. Функции подобраны так, что в одном из узлов она равна 1, в другом — 0. Такая система функций называется нормализованной. Вне элемента n на границе функции j” (s) обращаются в ноль. Заметим, что функции

(s) определены только на границе и здесь s не вектор, а скалярный параметр, меняющийся

вдоль границы.

После подстановки в (1) получим:

u. (X) = I1 j u* (X, x) d г x - j p* (X, x) (jr (x) + j (x У; )d г x |,

n=i [г” г” J

где u kj n — перемещение k -того узла элемента n в направлении j -того. Так как перемещения непрерывны, то uj = u‘n+1, j = 1,2. Для n +1 > N равенство преобразуется в u jN = u'1 — последний узел последнего элемента совпадает с первым узлом первого. Тогда соотношение для определения перемещений в произвольной точке области примет вид

( ^J

(2)

, V" Г+1 0J

Здесь u” — перемещение в направлении j -того узла с глобальным номером n . Можно считать, что глобальный номер узла равен номеру элемента, которому он принадлежит, или

u” = u 2 .

Записав (2) для всех узлов Xm, m = 1,2,..., N, получим систему линейных алгебраических уравнений, u'm = ut (Xm):

u.m = I j pn J u* ((, x )d Г x - u” fj p* (n, xj (x )d Г x + j p* (n, x)j”+1 (x)drx 11, m = 1,2,..., N .(3)

n=1 [г” у г” г*1 )J

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

(х)= j u* (х. x)d Г x;

Г”

P” (X) = J p* (X,x) j” (x)x + j+1 p* (X,x j+1 (x)drx .

Интегралы вычислены аналитически. Для вычислений используется локальная система координат. Начало системы совпадает с началом (первым локальным) узлом элемента n . Ось Z

N I

и< (X) = I p”n j u* (X, x) Г x - u” j p* (X, x) jj (x )d Г x + j p* (X, x) j”+1 (x )d Г x

направлена вдоль элемента, ось Z2 — по внешней нормали элемента. Тогда координаты точки X в новых осях примут вид:

X = px + Z cosa -Z2 sina, X2 = P2 + Z sina + Z2 cosa .

Здесь a — угол наклона элемента к оси üxx . По аналогичным формулам определяются координаты точки на элементе: xY = px + s cosa, x2 = p2 + s sin a, s — скалярный параметр.

В новых осях криволинейные интегралы превратятся в обычные.

Рассматриваются нижеследующие случаи взаимного положения точки X и элемента, по которому проводится интегрирование.

1. Точка X не лежит на элементе или его продолжении. В этом случае подынтегральные функции непрерывны, и интегралы — непрерывные функции параметров-координат точки X .

2. Точка X лежит на продолжении элемента. Локальная координата особой точки Z2 = 0 . Особенностей на отрезке интегрирования нет. Интегралы от функций ptí(X, x) равны нулю. А интегралы функций pij (X, x), i Ф j равны

Pkz , P„z (->)[L, + Z.ln(L-Ci)-ZilnZ, ]

P2(Zl) P2l(Zl) 4 (1 -n)pLn '

3. Точка X лежит на элементе, по которому ведется интегрирование. В этом случае подынтегральные функции ptj (X, x) имеют особенности. Интегралы вычисляются по схеме, изло-

женной в [1]. Интеграл по отрезку прямой заменим суммой трех интегралов:

Zi -e Ln 0

j (i e) = J pj ) )(x(s))ds + J pjj (X sj(x(s))ds + J pj(X, x(a)) p¡ ( ) e da.

0 Z\ +e -p

Третье слагаемое — интеграл по дуге малого радиуса e с центром в точке X, a — угол интегрирования. При интегрировании область, содержащая тело, остается слева и интеграл по дуге берется в отрицательном направлении. Так как радиус дуги мал, а функция (р, (x) ограничена,

то ее значение на дуге приближенно заменим ее значением в центре. Предел суммы при e — 0 и будет значением несобственного интеграла.

Для функций ptí (X, x) интегралы по прямым участкам обращаются в ноль, а предел интеграла по дуге равен

P, (Z, ) = -рр.

Для функций ptj (X, x), i Ф j интегралы по дуге обращаются в ноль. Пределы интегралов

Zi -e I,

Pjnk (ie)= J pj (X р (x(s))ds + J pjj (X, Р (x(s))ds

0 Zi +e

при e —— 0 принимаются значения

Pn, (z ) Pn, (z ) (1 - 2")[L- -(Ln -Z ) |n (L„ - Z, ) + (L„ - Z, )lní, ]

P21 P 4(1 -n)pL ’

n (4)

/ ч (1 -2v)[L + Z ln(L -Z )-ClnC ]

Pn2 Z \ = -Pn2 Z \ =' >L n b1 V n J

P21 P2 4(1 -n)pLn '

Заметим, что пределы lim f pj (Z1 ,Z2, (s)ds, i Ф j

Z ®00 j

совпадают со значениями интегралов,

вычисленных при ¡Х2 = 0 .

4. Точка X совпадает с одним из концов элемента. При й,\ ® 0 или й,\ ® ^п выражения (4) стремятся к бесконечности. Но сумма интегралов по соседним элементам (точка X лежит между ними) сходится. При вычислении будем использовать ту же схему, что и выше. Отличие

состоит только в том, что так как угол между элементами может быть любым, то интегралы по дуге будут зависеть от угла между элементами:

Ьп_е А.+1 ап+1

I Р* (X,5) $ (х0+ | р* (X, 5)^1п+1 (х(5))^у + | р* (X, х(а)) е ¿а, (5)

0 е /

где ап — угол положительного направления элемента п с осью х1 , а /п — угол обратного направления (от точки qn к точке рп) элемента п с осью х1 . Функция фпг (Ьп) = 1 и $ (0) = 1, поэтому в окрестности особой точки базисные функции приближенно заменили единицей. Интегралы

Ьп -е Ьп

I Р* (Х 5Ж (хО)№ , I Р* (X, 5Ж (х(5))й5

0е

имеют конечные пределы при е ® 0, так как $ (0) = 0 и ср[ (Ьп) = 0 .

Для функций ри (X, х) интегралы по прямым равны нулю, а интеграл по дуге равен

_ 411 _п)(ап+1 - ап + л) + 51п2ап+1 - «¿п2ап

8 (1 _п)л ’

если область выпуклая, и

_ 4 11 _П)(ап+1 _ ап ) + «¿п2ап+1 _ ^п2ап

8 (1 _п)л ’

если область не выпукла.

Для функций р12(Х, х) интеграл (5) равен:

оо8(2ап+1) _ оо8(2ап)

(1 _2'> ш-Ь 8(1 _п)р !„

а для функции р 21 (X, х) —

1п-

8л (1 _п) оо8(2ап+1) _ оо8(2ап)

8 (1 _П)Л Ьп+1 8л( _П)

Значения интегралов как функций параметра при аргументе X на границе не определены. При операциях с рациональными числами точные равенства, как правило, не выполняются. Поэтому для расчета значений интегралов вводится дополнительный малый параметр е . Считается, что если |^21 < е , то точка X лежит на элементе или его продолжении, и интеграл вычисляется по соответствующим формулам. Аналогично определяется совпадение с концом элемента: если выполняется неравенство |^1 | < е, то считается, что особая точка совпадает с началом

элемента, а если \^1 _ Ьп| < е — то

с концом. При расчетах е выбрано равным 10_7.

Подставив вычисленные интегралы в (3), получим систему линейных алгебраических уравнений относительно неизвестных узловых перемещений и усилий. Эта система решается численно.

Для анализа точности решения с помощью рассмотренного алгоритма была рассмотрена классическая плоская задача о растяжении прямоугольной полосы Ь х Н заданной силой (см. рис. 1). Кромка полосы х1 = 0 зафиксирована, к кромке х1 = Ь приложена равномерно распределенная

сила, направленная параллельно оси х 1. Кромки х2 = 0 и х2 = Н свободны. В середине полосы

имеется эллиптическое отверстие с полуосями а и Ь. На полосу не действуют массовые силы. Граница Г1 состоит из отрезка х1 = 0, 0 < х2 < Н , и (х) = 0, і = 1,2. Граница Г2 состоит из трех отрезков: х2 = 0, 0 < х1 < Ь, х1 = Ь,0 < х2 < Н, х2 = Н,0 < х1 < Ь . На первом и третьем отрезках сила равна нулю рі (х) = 0, і = 1,2, на втором задана сила р1 (х) = /, р2 (х) = 0 . На внутреннем контуре эллипса заданы нулевые силы, перемещения неизвестны.

При вычислениях напряжения и усилия масштабировались. За единицу выбрано значение модуля упругости. Поэтому значение относительного модуля упругости равно 1. Расчеты проводились для полосы с Ь = 100, Н = 1. Центр эллипса находился в точке с координатами

(50,0.5). Наибольшая полуось параллельна по оси Ох2 и равна а = 0.2.

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

Полученные решения сравнивались с теоретическими, полученными на основе упругих потенциалов [3,4]:

-{-|^4т3р2 - 4тр6 + р4 (3 + р4) + т2 (1 + 8р2 + 3р4)^00829 +

+ (1 + Р2)[2т3 + т4 -рб + тр1 (3 + р2) + т1 р1 (1 -р2) +

+тр2 (1 + р2)со84$Лу/ 2(т2 + р4 -2тр2оо82$) ; (6)

ар =-(р2 - 1){-2т3 -т4 + 3тр2 + тгрг -трА + тгрА --рб +[р4(р2 -3) + т2(3р2 -1)^00829-

-тр2 (р2 -1) 00Б 4^|у/|^2 (т2 + р4 - 2тр2 008 2$)

(7)

=

(р2 -і)sin2j|^2 (і-Р1)-2т3р2 + 3р4 + 2тр4 + р6 --Imp1 (і + р2)cos2jJ^т2 + р4 -2mp2cos2.9] . (8)

Здесь т = (a - b )/(a + b ) ; р — радиус вектор; J — полярный угол. Переменные р и J связаны с декартовыми координатами xi соотношениями

L a + b ( т | H a + b ( т | „

р----I cos J , x2 = —+- —І р+-I Sin J .

2 2 ^ р^ 2 2 ^ р

Точка р = 1 соответствует границе внутреннего контура. Сравнение проводилось в сечении, выделенном штриховой линией на рис. 1: х1 = 50, 0 < х2 < 0.5 - а. В этом сечении главные на-

пряжения в декартовых координатах совпадают с главными напряжениями в полярных.

Результаты для эллипса с отношением полуосей 1 представлены на рисунках. На рис. 2

2

представлена зависимость ар от р в сечении х1 = Ь/2. На рис. 3 — аналогичная зависимость для а 9. На рис. 4, 5 представлены разности расчетных напряжений и вычисленных по формулам (6), (7).

На рис. 6-9 представлены аналогичные результаты для эллипса с отношением полуосей 10 . Для эллипсов с большим соотношением осей погрешность достигает 5%.

Р и с. 2 Р и с. 3

Р и с. 4 Р и с. 5

Р и с. 6 Р и с. 7

Р и с. 8 Р и с. 9

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

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

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

3. Оден Дж. Конечные элементы в нелинейной механике сплошных сред. М.: Мир, 1976. 464 с.

4. Седов Л. И. Механика сплошной среды. Т 2, М.: Наука, 1983. 420 с.

5. Тимошенко С.П., Гудьер Дж. Теория упругости. М.: Наука, 1979. 560 с.

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