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

Вестник СамГУ — Естественнонаучная серия. 2007. №2(52) 19

МАТЕМАТИКА

УДК 539.3

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

© 2007 А.А. Абдрахманова2

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

1. Постановка задачи

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

С2 / с^ш\

г?К?Нг)- (1)

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

Производные от функции перемещения ш = ш(х) определяют деформацию стержня и действующие в нем внутренние силовые факторы ([1]): первая производная от функции перемещения — это угол поворота поперечного сечения; вторая производная определяет изгибающий момент и возникающие в балке нормальные напряжения; третья производная характеризует поперечную силу и касательные напряжения; четвертая производная тесно

1 Представлена доктором физико-математических наук профессором Л.С. Пулькиной.

2Абдрахманова Алия Айдаровна , кафедра математики Уфимского государственного авиационного технического университета, 450000, Россия, Уфа, ул. К. Маркса, 12.

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

Такому требованию удовлетворяют сплайны пятой степени дефекта 1, методика построения и применения которых для решения одномерных задач механики деформируемых твердых тел изложена в работе [2].

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

В предлагаемой работе представлена новая схема применения метода сплайнов, обеспечивающая шестой порядок сходимости.

2. Основные положения метода сплайнов пятой степени

При построении сплайна пятой степени дефекта 1 на отрезке [а, Ь] формируется сетка А , имеющая N узлов

А : а = х1 < х2 < ... < xN = Ь.

На сетке А по методике, изложенной в [2], строится сплайн-функция ^5д(х) степени 5 дефекта 1, имеющая

N. = N + 4

степеней свободы, равных минимальному числу параметров, однозначно определяющих сплайн-функцию ^5д(х).

В пределах каждого отрезка [хг-,хг+1],i = 1,...^ - 1 функция ^^(х) является многочленом пятой степени

5

^5,1(х) = ^ а‘а (х - х,)а, х е [х,, х+1], i = 1,...^ - 1, (2)

а=0

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

й* Т^5 1(х, - 0) й* ^5 1(х, + 0)

-----^-------- =----^5 = 0,1,...,4, 7 = 2,...,Ж-1.

йх* йх*

В работе [2] параметры, определяющие сплайн, сведены в вектор-столбец параметров сплайна

Я = (дк ,к = 1,2 ,..., N + 4)т ,

где

41 =

42 =

сіх4 ’

¿2ТУ5,і(хі) ¿/х2 ’

4+2 = ^5,і(хг), і = 1,2,...,М,

с2 ^5,і(хм)

4м+ъ =

4м+4 =

Сх2

с?\У5Л{хн)

сіх4

Число компонент вектора Я точно соответствует числу степеней свободы сплайна N. = N + 4.

Согласно (2), в пределах каждого отрезка [х,-,хг+1], i = 1,...^- 1 функция ^5д(х) однозначно определяется коэффициентами

а®, а = 0,1,...,5, i = 1,...^ - 1,

которые сведены в векторы-столбцы

Аа = (а®, i = 1,...^ - 1)Т , а = 0,1,...,5.

Векторы Аа однозначно определяются [2] через вектор параметров сплайна Я матричными выражениями

Аа = СаЯ, а = 0,1,...,5. (3)

Матрицы Са, а = 0,1,..., 5

Са =

с(а),

г>}

і = 1,...,М - 1, } = 1,...,М + 4, а = 0,1,..., 5

зависят от вида сетки узлов А . Методика построения данных матриц изложена в работе [2].

На основе (2) определяются производные от сплайн-функции ^5д(х) в любой точке с координатой х из области определения [а, Ь]

^ УУ5Л(х) ¿/х5

а!

(а - 5)!

аа\х - хі)а 5,

(4)

х є [хі, хі+1], і = 1,...,М - 1, 5 = 0,...,5.

В пределах каждого отрезка [хі, хі+1] функция ^5д(х) является многочленом пятой степени с коэффициентами, определяемыми в соответствии с (3) следующими выражениями:

М+4

аа = 2 С

7(і) _

4а

,(а) а ■ і,} (i},

і = 1,...,М - 1, а = 0,1,...,5.

(5)

}=1

При подстановке (5) в (4) получено

^ \У5і1(х) ¿/х5

а!

=5 (а - 5)!

М+4

2 с а}

}=1

(х - х )

а - 5

х є [хг, хі+1], і = 1,...,М - 1, 5 = 0,...,5.

При изменении порядка суммирования в (6) получено

т5 ТТГ / \ М+4 л УУ5Л(х) = ^

}=1

Сх

а!

->(«) (а - 5)! 1>}

С}х - х,)

Ч},

х є [хі, хі+1], і = 1,...,М - 1, 5 = 0,...,5.

= І / а! М с(а)(х - х,)а - 5,

.Л' а=5 (а - 5)! ч к и

х є [хі, хі+1], і = 1,...,М - 1, 5 = 0,...,5, } = 1,...,М + 4.

С учетом (8) выражение (7) приняло вид С5 Иг5,1(х)

(7)

В соответствии с (7) введена вектор-строка

Д*(х) = (Д* (х), ] = 1,...^ + 4) , х е [хг, хг+1], * = 0,...,5, (8)

компоненты которой определяются следующими выражениями

(9)

= Я5(х)Я, х є [хі, хі+1], і = 1,...,М - 1, 5 = 0,...,5. (10)

Соотношения (10) позволяют определять значения сплайна ^5д(х) и его производных в произвольной точке с координатой х из области определения.

3. Дискретный аналог уравнения равновесия

Дифференциальное уравнение равновесия (1) записано в развернутой форме

с14со 1 ^{Е1у)сРсо 1 & {Ыу) ср-ш

= дг(х).

ЕЕ

+ 2-

+

у йх4 dx йх3 йх2 йх2

Заменой искомой функции перемещения осевой линии ш = ш(х) на ее сплайновую аппроксимацию ^5,1 (х) определена поперечная нагрузка 7[г = = 7[г(х) , при которой сплайн ^5д(х) является точным решением дифференциального уравнения (11). При этом подстановкой (10) в (11) получен сплайновый аналог дифференциального уравнения равновесия.

Чі (х) =

Е1у #4( х) + 2

С(Е1у) С2 (Е1у)

Сх

■Яз( х) +

Сх2

#2( х)

Я.

(12)

На основе (12) введена вектор-строка

К( х) = Е1уЯ4(х) + 2

СІЕІу) С2 (Е1у)

Лг^з (*) + -ЬИд2М,

Сх ъ Сх2

и уравнению (12) придан более лаконичный вид

%(х) = К( х)Я.

При построении дискретного аналога дифференциального уравнения (11) использован полученный в теоретической механике результат [3]: силы,

(13)

(14)

произвольно расположенные в пространстве, можно привести к одной силе, равной их главному вектору и приложенной в центре приведения, и к паре сил с моментом, 'равным главному моменту всех сил относительно центра приведения.

Базируясь на этом положении, была выбрана сетка узлов А с четным числом отрезков [хг-,хг+1],} = 1,...,N—1 при N = 2М+1, М = 1,2,... В пределах этой сетки рассмотрены узлы с номерами п:

п = 2т, т = 1, ...,М, при М = ^ — 1)/2.

Данные узлы, имеющие координаты хп , выбраны в качестве центров приведения для заданной дг = дг(х) и приближенной 1дг = ~дг(х) распределенных нагрузок.

Для определения сплайн-функции ^5д(х), близкой к искомому решению ш = ш(х), использовалось условие эквивалентности внешних действующих на балку нагрузок, заключающееся в том, что в пределах каждой пары отрезков [хп-1, хп] и [хп, хп+\] распределенные нагрузки дг = дг(х) и

& = 'Цг(х) приводились к центру приведения, имеющему координату хп. Получаемые при этом на основе точной дг = дг(х) и приближенной 1дг = ~дг(х) нагрузок главный вектор сил и главный момент приравнивались. В итоге получалась система из N — 1 уравнений, имеющая следующий вид:

хп+1 хп+1

f ~дг(х)^х = f дг(x)dx,

хп— 1 хп— 1

хп+1 хп+1 (15)

f ~дг(х)(х — хп^х = f дг(х)(х — хп^х,

хп—1 хп—1

п = 2т, т = 1,...,М, М = ^ — 1)/2.

При подстановке (14) в (15) получена система из N — 1 алгебраических линейных уравнений

'хп + 1 \ хп+1

f К(x)dx Ц = f д?( х^х,

хп—1 / хп— 1

(16)

'хп + 1 \ хп+1

f К(х)(х — хп^х Ц = f дг(х)(х — хп^х,

ХХп— 1 / хп— 1

п = 2т, т = 1,...,М, М = N — 1)/2.

Вектор параметров сплайна Ц имеет N + 4 компонент, для определения которых необходимо иметь N + 4 уравнений.

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

?г( х1) = дz (хО. (17)

С учетом (14) уравнение (17) принимает вид

К(х)Ц = дг (хО.

Недостающие еще 4 уравнения должны формироваться на основе учета краевых условий.

В итоге получается система из N + 4 алгебраических уравнений, однозначно определяющая вектор 3 и, следовательно, саму сплайн-функцию WЪ,l{x).

4. Модельная задача, имеющая точное решение

При изучении точности метода сплайнов рассматривалась частная по отношению к дифференциальному уравнению (1) модельная задача об изгибе нагруженного распреде-ленной нагрузкой qz, защемленного на концах стержня длиной € = 2 м. Балка считалась стальной и имеющей постоянный по ее объему модуль упругости E = 2 ■ 1011 Па= const. Поперечное сечение имело прямоугольную форму высотой h = 0,02 м и шириной b = 0,03 м. При этих данных, согласно [1], рассчитан осевой момент инерции поперечного сечения

bh3 о

1у = -j-^- =2-10 м = const.

При постоянных модуле упругости E = const и осевом моменте инерции сечения Iy = const уравнение (1) приняло вид

d4 ш

Ely Ihd = ^

Для защемленной по концам балки краевые условия имеют вид [1]

со = ^ = о при х = 0 и х = £. (19)

dx

Рассмотрен случай изменения поперечной нагрузки по степенному закону д2 = Лхк,

к = 0,1,..., где величина коэффициента Л определялась по формуле Л = 6 ■

• 102^ .

€к

При совместном рассмотрении (18) и (19) записана система расчетных уравнений

й^ш и

Е1У1?=Л^,к=й.\....

йш

при со = —— = 0 ДЛЯ X = 0 И X = -с,

йх

имеющая точное аналитическое решение следующего вида:

хк+4 — (к + 2)£к+1 х3 + (к + 1)£к+2 х2 , Л ,

со =---------------------------------, к = 0,1,... . (20)

(к + 1 )(к + 2 )(к + 3 )(к + 4)

На основе (20) определяются производные от ш по x

i/co _ A (к + А)хк+Ъ - 3{к + 2)€к+1х2 + 2(к + \)Єк+2х dx Ely (к + 1 ){к + 2 )(к + 3 ){к + 4)

d2ш A {к + 3)(k + 4)x^2 - 6(k + 2)€k+1 x + 2(k + \)fk+2

dx2 EIy (k + 1)(k + 2)(k + 3)(k + 4)

d3№ A (k + 2)(k + 3)(k + 4)xk+1 - 6(k + 2Xk+1

dx3 Ely (k + 1)(k + 2)(k + 3)(k + 4)

Определение главных векторов сил и главных моментов

Главный вектор сил Fn от распределенной нагрузки qz = Axk, действу-

ющей на отрезке [хп—1, хп+1], определяется выражением

хп+1 1 / ч

Е„ = / Ах^х = А^—^ ((хп+1)к+1 -(х„_і/+1)

Хп— 1

при п = 2т, т = 1,...,М, М = N - 1)/2.

Главный момент сил Мп от распределенной нагрузки д2 = Ахк, действующей на отрезке [хп-і, хп+і], при приведении к центру с координатой хп определяется выражением

хп+1

Мп = f Ахк(х — хп)^х =

хп— 1

при n = 2m, m = 1,... ,M, M = (N - l)/2.

Для уравнения (18) с постоянными коэффициентами E = const и Iy = = const согласно (13), получается

K(x) = EIyR4(x).

Согласно (9), определяются компоненты вектора R4(x)

(21)

(22)

x є [x„ x/+1], i = 1,...,N - 1, j = 1,...,N + 4,

или

R4(x) = 4!Сг(4 + 5!C(,5y)(x - xt),

x є [xi, xi+1], i = 1,...,N - 1, j = 1,...,N + 4.

В соответствии с (16) формулируются векторы-строки и ¿М)

' 4Р = ], ] = + 4),

- ¿М) = (]’п), ] = 1,...,И + 4) (24)

при п = 2т, т = 1,...,М, М = (И - 1)/2,

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

хп+1

¿Р = / К( х)йх,

X п-1

хп+1 (25)

¿М) = / К(х)(х - хп)йх

X п-1

при п = 2т, т = 1,...,М, М = (И - 1)/2.

В соответствии с (21) - (25) определяются выражения для вычисления значений компо-нентов векторов и ¿М)

хп

¿рп) = Е1у / (4!Сп4.)1,. + 5!С п5)1,у(х - хп-1)) йх+

х п-1

хп+1

+Е1у / (4! С п4) + 5!^] (х - хп)) йх, ] = 1,...,И + 4,

х

хп

¿М'п) = Е1у / (4!Сп4_)1,] + 5!Сп-)1,](х - хп-1)) (х - хп)йх+

хп+1

х п-1

+Е1у / (4!Сп4]) + 5!С(^])(х - хп)) (х - хп)йх, ] = 1,...,И + 4 х

при п = 2т, т = 1,...,М, М = (И - 1)/2.

5. Методика оценки точности результатов численных расчетов

Численное решение уравнения (18) строилось на сетке Оц постоянного шага к с N узлами

Ок = {х = 1к, к = €/(И - 1), I = 1,2,... ,И}.

Расчеты выполнены при числе узлов сетки N = 5,11,21,31,51,101 для значений коэффициентов к = 0,1,2,5,10,20,50,100.

Точность оценивалась по максимальным относительным ошибкам 65, 5 = = 0,1,2,3 значений функции ш = й0ш/йх0 и ее производных й5ш/йх5, 5 = = 1,2,3 в узловых точках х,, I = 1,...,И, вычисляемых по формуле

6( 5) = тах

г=1,...,И

й5 ш, й5 шТ

йх

йх

й5шТ

тах

5 = 0,1,2,3,

где

точное значение производной порядка 5 от функции ш в рас-

й5ш,

сматриваемои точке с координатой хг-,

йх

расчетное значение производ-

ной по рядка 5 от функции ш в рассматриваемой точке с координатой х й5 шТ й5шТ1

йх

тах

= тах

г=1,...,М

йх

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

6. Обсуждение результатов

На рис. 1-4 в координатах ^ ^65) ~ - 1), 5 = 0,1,2,3 представлены

зависимости ошибки расчетов от размерности сетки узлов N и показателя степени к для функции перемещения ш = ш(х) и производных от нее.

Из рис. 1 видно, что при к = 0 и к = 1 для сетки с числом узлов N = 5 относительная погрешность расчетов для функции перемещений ш = ш(х) характеризуется величиной 6(0) < 1 ■ 10-13, то есть в данном случае метод дает практически абсолютно точные результаты. В дальнейшем при росте N для случаев к = 0 и к = 1 погрешность расчетов возрастает, и это связано с накоплением арифметической ошибки расчетов. Чтобы эта ошибка увеличивалась менее интенсивно, необходимо более тщательно прорабатывать алгоритм арифметических вычислений.

В случае к > 1 погрешность расчетов с увеличением числа узлов сетки N монотонно уменьшается, что свидетельствует о сходимости метода. Из рис.1 видно, что при определении функции перемещения ш = ш(х) наблюдается шестой порядок сходимости. В случае N > 51 для к = 2,5,10 имеет место уменьшение точности результатов расчетов. Это связано с влиянием ошибки арифметических расчетов. Таким образом, для повышения точности расчетов при высоких значениях N необходимо повышать точность вычислительной процедуры, что можно сделать, например, вычислениями с большим числом значащих цифр.

На рис. 2 представлены результаты расчетов для первой производной йш/йх от функции перемещения ш = ш(х). Наблюдается тот же шестой порядок сходимости метода.

На рис. 3 представлены результаты расчетов для второй производной й2 ш/йх2 от функции перемещения ш = ш(х). В данном случае метод имеет четвертый порядок сходимости.

На рис. 4 представлены результаты расчетов для третьей производной й3 ш/йх3 от функции перемещения ш = ш(х). Для третьей производной метод обеспечил четвертый порядок сходимости.

Таким образом, результаты численного эксперимента показали, что предлагаемый в работе метод численного решения дифференциального уравнения четвертого порядка, описывающего изгиб балки, обеспечивает при вычислении самой функции ш = ш(х) и первой производной от нее d^/dx шестой порядок сходимости, а при расчете второй й^ш/dx2 и третьей d3ш/dx3 производных — четвертый порядок сходимости.

Литература

[1] Работнов, Ю.Н. Механика деформируемого твердого тела / Ю.Н. Ра-ботнов. - М., Наука, 1979.

[2] Павлов, В.П. Метод сплайнов и другие численные методы решения одномерных задач механики деформируемых твердых тел / В.П. Павлов. -Уфа, Уфимск. гос. авиац. техн. ун-т, 2003.

[3] Яблонский, А.А. Курс теоретической механики. Ч.1. Статика. Кинематика / А.А. Яблонский , В.М. Никифорова. - М., Высшая школа, 1977.

Поступила в редакцию 13/XT7/2006; в окончательном варианте — 13/X///2006.

A SPLINE METHOD VARIANT FOR BEAMS BENDING3

© 2007 A.A. Abdrakhmanova4

A numerical solution method for differential equation of the fourth order is considered, describing a bending of a beam. A spline-method in the integrated form, basing on spline-functions of the fifth degree and providing sixth order of convergence, is discussed. For estimation of accuracy of the considered method a modelling problem having the exact analytical solution is used.

Paper received 13/XTT/2006. Paper accepted 13/X///2006.

3Communicated by Dr. Sci. (Phys. & Math.) Prof. L.S.Pulkina.

4Abdrakhmanova Aliya Aidarovna, Dept. of Mathematics, Ufa State Aviation Technical University (USATU), Ufa, 450000, Russia.