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

УДК 539.3

М. К. Сагдатуллин

ПОСТАНОВКА ЗАДАЧИ ИЗГИБА ПЛАСТИН НЕПРЯМЫМ МЕТОДОМ

ГРАНИЧНЫХ ЭЛЕМЕНТОВ

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

В данной работе описывается применение НМГЭ к решению задачи изгиба пластин произвольного контура в линейной постановке. Приводятся интегральные уравнения метода компенсирующих нагрузок для различных способов закрепления пластин. Получены численные формулы интегрирования по элементам контура пластины. На примерах решений тестовых задач показана эффективность используемого метода.

Keywords: a boundary element method, deflection plate, loads balancing method, system of linear equations.

This paper describes the use of IBEM (indirect boundary element method) algorithm to solve the problem of bending of plates of arbitrary contour in the linear formulation. We give integral equations method of compensating loads of different ways of fixing plates. Numerical integration formulas on elements of plate contour. For examples of test solutions of the problems shows the effectiveness of the method used.

Введение

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

Огромный вклад в изучение и развитие методов граничных интегральных уравнений был сделан Купрадзе В.Д., Мусхелишвили Н.И., Михлиным С.Г., Смирновым В.И. и др.

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

Численной реализации методов потенциала в задачах теории упругости посвящены работы Партона В.З., Крауча С.Бреббия К., Угодчикова А.Г., Перлина П.И., Бенерджи П, Уокера С., Баттерфилда Р., Старфилда А., Круза Т., Вроубела Л., Риццо Ф., Громадки Т., Лей Ч., Кузнецова С.В., Лившица И.М., Розенцвейга Л.Н. и др. [1-3].

Метод граничных элементов (МГЭ) - это метод численного решения граничных интегральных уравне-

ний с дискретизацией границы области и заданием аппроксимации функций плотности на границе. Интерес учёных к применению МГЭ связан с достоинствами этого метода: отсутствием ограничений на геометрию контура, снижением на единицу размерности рассматриваемой задачи, аналитическим описанием особенностей решения, высокой точностью результатов решения.

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

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

Непрямой метод граничных элементов в задачах изгиба пластин известен как метод компенсирующих нагрузок. Функциям плотности придается смысл нагрузок, приложенных к бесконечной пластине и распределенных по границе области, или по некоторому контуру, внутри которого находится область. В задачах изгиба пластин первые работы в данном направлении выполнены Кореневым Б.Г. и дальнейшее развитие этот метод получил в работах Артюхина Ю.П., Грибова А.П., Венцеля Э.С.,Толкачева В.М., и др.

Исходные соотношения и гипотезы

В тех случаях, когда прогибы w пластины малы по сравнению с её толщиной И, имеется возможность построить вполне удовлетворительную приближенную

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

1. В срединной' плоскости пластина не испытывает никаких деформаций. При изгибе эта плоскость остается нейтральной.

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

3. Нормальными напряжениями в направлении, поперечном к срединной плоскости пластины, допустимо пренебрегать.

Основываясь на этих гипотезах, мы сможем все компоненты напряжений выразить через прогиб V пластины, являющийся функцией двух координат в плоскости пластины. Эта функция должна удовлетворять линейному дифференциальному уравнению в частных производных, которое, вместе с граничными условиями, полностью определяет V. Расположим координатные оси х и у в срединной плоскости пластины, ось г направим перпендикулярно к этой плоскости. Дифференциальное уравнение строится при рассмотрении элемента, выделенного из пластины двумя парами плоскостей, параллельных плоскостям хг и ух. Через Мх обозначим отнесенный к единице длины изгибающий момент, действующий по краям, параллельным оси у, через Му - момент, приложенный по краям, параллельным оси х, Мху -крутящий момент, ((х и Оу - вертикальные поперечные силы, приложенные по боковым граням элемента. Используя закон Гука, запишем выражения моментов через прогиб пластины V:

Mx = -D

(д 2w д 2w^

-— + V-—

дх£

My = -D

(д 2w

+V

dyÁ

ду2

dW

дх2

(1)

M

ху

■■ -D(1 - V)

д 2 w дхду

где D =

Eh"

- жесткость пластины при изгибе, Е

12(1 - V2)

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

д 2Mx

д 2M

y

дх

2

ду

2

+ 2-

д 2M

ху

дхду

= -q.

(2)

Из условий по моментам так же можно получить выражения для поперечных сил

°х = -^дх дх

(а 2

2... \

д 2w д 2w +

дх

2

ду

°у =-Dду

2... \

( д 2w + д2 w дх2 ду2

(3)

Подставляя выражения (1) в уравнение равновесия (2), найдем

д 4w

44 „ д w д w q ,.ч

+ 2—--- +-- = 4г . (4)

дх4 " дх2ду2 ду4 D Уравнение (4) представляет собой уравнение в частных производных 4-го порядка. Его можно записать также в символической форме

AAw = — , А = D

д2 д2

дх2 ду2

Чтобы краевая задача для уравнения 4-го порядка имела определенное решение, на контуре должны быть заданы граничные условия. Если край пластины заделан, то граничные условия на нем имеют вид

п дw

V = 0, — = 0; (5)

дп

при шарнирно опертом

V = 0, М = 0; (6)

при свободном

М = 0, О * = 0; (7)

здесь М - нормальный изгибающий момент, О * -

приведенная поперечная сила на контуре

* дМ

Q* =Q ,

д s

n , т - нормальное и касательное направления к контуру пластины, s - дуговая координата контура. Силовые факторы Q , М , М можно выразить через

Т 1 ¿L^ft 7 n 7 пт А А

Qx, Qy, Мх, Му, My по следующим формулам Q = Q cosf + Q sinf ,

^ n ^ x J У

M = M cos2 f + M sin2 f + M sin2f , (8)

n x J У J ХУ J J v '

Mn = Mху cos

2ф +1 (mу - Mх )sin

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

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

Mn = -D

д 2 w

+ V

дп2 д 2 w

+ V

д 2 w д г2

д 2 w

д г2

2

Mn =-D(l - v)

дп

д (дш д т\дп

Qn =-D — дп

22 д w д w

-Г- +-Г-

д г¿

Qn

-Di — I дп

дп2

( д 2 w д 2w^

-Г- +-г-

дп2

д г2

+

Мт = -D

+

+(1 -V)

д 3\

дпдг

(д 2\ дп2

д V

' д г2

где р - радиус кривизны контура.

Метод компенсирующих нагрузок при изгибе пластин

Рассмотрим тонкую изотропную линейно-упругую пластину, срединная поверхность которой занимает область О, ограниченную гладким контуром Г. Дополним область О областью О- до бесконечной области. По контуру Г к бесконечной пластине приложим компенсирующие нагрузки р(£) и т (9 . Нагрузка р(£) -распределенное по контуру Г усилие, нормальное поверхности пластины, т (9 - распределенный по контуру Г момент вокруг касательной к контуру Г (рис. 1).

Рис. 1

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

АА\ = ^ D

(10)

Решение уравнения (10) ищется в интегральном виде

ад №«ю

дп

оГ( 0

(11)

\кг{)) = \\в{),{,Ш9 ЛХZ)+XG^()Zi)F(/ ■ (12) о /

Здесь G (^9 - фундаментальное решение для задачи изгиба бесконечной пластины, t (х, у) - точка области О, П) - точка контура Г:

1 г2 |

в),0=-г2 1п г

8кО

(13)

г = ^х-92 +(у-П2 ленной нормальной нагрузки, Fi - модуль сосредоточенной силы, приложенной в точке £ i 0' = 1,2,...).

Компенсирующие нагрузки определяются из граничных условий (5)-(7) на контуре пластины. Граничные условия записываются на контуре Г для области О.

Отметим, что выражение (13) удовлетворяет дифференциальному соотношению (10). Чтобы (13) было решением краевой задачи, необходимо из системы сингулярных интегральных уравнений, которая получается при подстановке (13) в краевые условия (5)-(7)

д(9 - интенсивность распреде-

на контуре пластины, определить функции р(9 и т (9. При этом надо произвести предельный переход точки t(х,у) из внутренней области О на границу Г.

В области О изгибающие и крутящий моменты по известным р(С), т (9 определяются следующим образом

Мх() = -й\

д 2в д Щ

—— + V-

дх2

ду2

А 0-

( дЩ

- + V-

дЩ ^

дх2дп ду2дп

«9

0Г( 0+ М/))

Му()) = ^

д2Щ дЩ

- + V——

ду 2

дх

А 9

( д3С

- + V-

д3С ^

ду2дп дх2дп

«9

Му(/)=-Ц1 - V) J

Л о +Муг() А 9 -

д Щ

дхду'

д3Щ !-

дхдудп

«о

оГ( о + Мху^()

(14)

где

МГО = -0

М/)) = -Я

^ д 2\г() д 2\г() ^

дх

ду

д2\г() ^ -+ V-

ду£

дх'

М хуГ) = -0(1 - V)

д 2\г)) дхду

(15)

Граничноэлементное представление интегральных уравнений метода компенсирующих нагрузок

Полученная после подстановки (10) в краевые усло вия (5)-(7) на границе пластины система интегральных уравнений в общем случае не имеет ана-литиче ского решения. Реализация дискретного подхода решения состоит из следующих шагов.

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

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

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

4) Интегральные уравнения, полученные из граничных условий, сводятся к системе линейных алгебраических уравнений.

5) Значения функции V в области W и на контуре Г определяются соотношением (11).

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

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

шгО

виор о-^Я то

дпПО

г

дЩОд

дпо

дп{ 0

Р тХО

дп© дПО

С 0 = 0, (16)

С 0 = 0,

(17)

где t (х , у ) - точки на границе пластины Г, п (£) и п ^) - внешние нормальные направления к границе пластины в точках £ и / соответственно. Следует отметить, что (16) и (17) являются точным представлением решения задачи. При численном решении любая погрешность в конечном результате получается исключительно из-за дискретизации интегралов и последующего решения алгебраических уравнений.

Разбив контур Г на N граничных элементов, для отдельных узловых точек каждого элемента уравнения (16) и (17) записываются в дискретном виде

N

w

"(О ) + £ ¡рР 0 с )0Т( О-

У=!г,

N

"ХИО

У=1г,

дП 0

сТ( 0 = 0 :

(18)

дП,

N

У=1г,

N

■и Р О

У=1г/

дп( о дпо,-)

дп,

сг( 0 = 0.

С 0-

(19)

Здесь интегрирование ведется в пределах у - того элемента, t. - точка наблюдения (/ = 1,2,...Интегралы

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

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

Рис. 2

Таким образом, функции р(С) и т (£) могут быть вынесены из-под знака интеграла, что дает следующую систему линейных алгебраических уравнений

N

N

А, + £руву - £ туНу = 0 , / = 1,2.....N , (20)

У=1

N

У=1

N

В/ + £руЯу - £туМ,

У

0

/ = 1,2.....N,

(21)

1=1 У=1 где неизвестными являются дискретные значения компенсирующих нагрузок на граничных элементах р. = р(^ ), т . = т), коэффициенты в уравнениях имеют вид

в У = ¡вО/, ^)<сг( о

НУ = 1

дв(у, 0

дп{ О

С ():

Я/

= 1

двО/, £) дп(У)

С о

Му = \-ЙпС О-

1 г ПОдпО;)

А/ = wr{ti)

В/

дшгО/)

дЦО)

(22)

(23)

(24)

(25)

(26)

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

-+ V-—

дп (О

дг2( О

'д2в(),0 , V д2в(),о Л

дп2 (О

дг2( О

РР о-

' д3в()0 д 3в(), ¡3

-+ V-т(0

дп( 0дп2()) дп( Одг^ (О у

ссг = 0.

Подынтегральное выражение последнего уравнения является потенциалом двойного слоя и содержит в себе особенность t = £. Выделив некоторую окре-

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

д2^))

d2wr{f) mtt

dn2 )t)

дт2 )

diG)t,Q

+ v-

- mw + f

2 J

Г

d3G[t,Q

Г d20)t, Q

d2G)) Z Л

дп2))

дт2 W

дгфдп2)) дфдт2))

mKQ

d - 0,

где интегрирование ведется по контуру с выколотой особенностью и понимается в смысле главного значения по Коши. В данной работе вычисление интегралов с особенностями будет проводиться на дискретной модели контура.

В дискретном виде условие равенства нулю нормального изгибающего момента на границе примет вид N N

С/ + XР]Ти - Е т1Би = 0 ' / = 12.....М , (27)

У=1 У=1

где

Tj-i

д G{tj, Q дп2( j)

. д2G(t/,Q ^

дт2( t/)

d,

Г

д G{t/, Q

■ + v-

д G{t/,Z

С/ -

д/лодп^о) д/ 0дт2(/)

д2w'(t/) д w' (t¡)

+ v-

оГ,

д/1/)

дт2 (t/)

(28)

(29)

(30)

Учет действия равномерно распределенной нагрузки заключается в представлении ^) в (11) как частного решения уравнения (10):

(31)

48D

Числовые примеры

Настоящий раздел посвящен примерам численной реализации решения задач. Рассмотрено несколько тестовых примеров для апробации методики, представленной в предыдущих разделах, приводятся сравнения с аналитическими решениями задач. Задача 1. В качестве примера была выбрана квадратная пластина.

q - равномерно-распределенная нагрузка, ц- 2 106, а -сторона квадрата, а = 1, D - жесткость пластины при

изгибе, D - Е •

12(1-V2)

Результаты прогиба в центре квадратной пластины приведены в таблице 1. Таблица 1

Пост. элемент Линейный элемент

8 узлов 0.0015722422347459 0.0015585241670949

24 узла 0.0011151149982027 0.0011949410724905

40 узлов 0.0011104006543812 0.0011419081694110

аналит. решение 0.00110074 0.00110074

Аналитическое решение [7] для прогиба в центре квадратной пластины следующее:

,4

W - 0.00126 •

q • a D

Задача 2. В качестве примера рассмотрим круглую пластину.

q - равномерно распределенная нагрузка, ц - 2 • 106 , а - радиус круга, а = 1, Б - жесткость пластины при изгибе, й - -

E • h

12(1 -V2 ) '

Аналитическое решение [7] для прогиба круглой пластины в центре следующее:

4

W-

q • a

64 D

Далее в таблице 2 получим прогиб в центре пластины для круга (случай постоянного и линейного элемента) при разном количестве узлов и сравним с аналитическим решением.

Таблица 2

Пост. элемент Линейный элемент

8 узлов 0.0101531248198796 0.0114368921000835

24 узла 0.013348671920581 0.0134961333868765

40 узлов 0.0135029600389354 0.0135702357495802

аналит. решение 0.01365 0.01365

Заключение

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

Литература

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

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

3. Крауч С., Старфилд А. Методы граничных элементов в механике твердого тела. М.: Мир, 1987. -328 с.

4. Голованов А.И., Бережной Д.В. Метод конечных элементов в механике деформируемых твердых тел. Казань, 2001, 300 с.

5. Бережной Д.В., Сагдатуллин М.К., Саченков А.А. Универсальный конечный элемент для расчета комбинированных конструкций // Вестник Казанского технологического университета. - 2012. - №17. - С.150-157.

6. Бережной Д.В., Сагдатуллин М.К. Трехмерный конечный элемент для расчета оболочек средней толщины // Вестник Казанского технологического университета. -2013. - №9. - С.256-261.

7. Тимошенко С.П., Войновский - Кригер С. Пластинки и оболочки. М.: Наука. - 1966. - 636 с.

© М. К. Сагдатуллин, канд. физ.-мат. наук, доц. каф. теоретической механики и сопротивления материалов КНИТУ, [email protected]. © M. K. SagdatuIIin, PhD, associate professor, KNRTU, department of theoretical mechanics and strength of materials, [email protected].