К расчету частот и форм колебаний балки с произвольным числом упруго закрепленных тел

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

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

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

1    Получение уравнения на собственные частоты и соответствующих им собственных форм колебаний

Рассмотрим механическую систему, представляющую собой однородную упругую балку длины l, плотности р, модулем упругости E, моментом инерции J поперечного сечения балки относительно нейтральной оси сечения, перпендикулярной плоскости колебаний, с закрепленными на нем nтелами массами mi,i = 1,...,n , посредством пружин с жесткостя ми ci, i = 1,...,n . Концы балки жестко закреплены. Массы mi могут перемещаться вертикально в направлении осей zi. Колебания масс характери зуются функциями zi (t), перемещения точек оси балки описываются функцией u (x, t) [13].

Уравнения движения механической системы имеют вид: d ² z

~г+рЛz. -u(a,t)) =0, d2u , ,54u                        _      .

tt + b — = J e , ⁽ zi ^- и ⁽ x , t >⁽ x ^- a ), d t        d x      i = i

• 2    c к EJ где Pi = — , b = — , e = m i     p F

  c_i

  p F

  
  

  .

  
  

Краевые условия на концах балки:

u (0, t ) = u ( l , t ) = 0, d u _ d u _ — (0, t ) = —( l , t ).

d x       d x

Представив zi(t), u(x, t) в виде zi (t) = At sin( ^t + ai), u(x, t) = V(x) sin(mt + в), подставив в (1), после преобразований получаем

• -    _Ю ² A , + p , ²( A , - V ( a , )) = 0

• -    го2 V(x) + b d V4x) = J e,(A, - V(x))S(x - a,) dxi

с краевыми условиями на концах стержня

V (0) = V (l) = 0, dVdV

—(0) = -- ( l ) = 0.

dxdx

Ввиду линейности второго уравнения в системе (3) его обобщенное решение имеет вид n

V (x) = J V( x - a,) e ,(A,- V (a,)),(5)

i = 1

где функции V , -( x ) являются решениями краевых задач

• - to ² V i ( x ) + b d V¹ 4 ⁽ x ) = A ( x ) , dx ⁴

  
  

V (-a = V(l - ai) = 0, dV      dV

^{( - a}i ) =      ⁽ l ^{- a}i ) = ⁰ .

dx        dx

Подставляя в (5) последовательно x = a i , получаем систему линейных алгебраических уравнений относительно V ( a_i ) :

(1 + V (0) e )V ( a ) + V ( a - a ₂) e ₂ V ( a 2 ) +

+ ... V ⁽ a ^{- a} n ) e n ^{V ( a} n ) =

=ZV(a- a) eA, i=1

V ( a 2 - a)eV ( a) + (1 + V(0)e₂)V ( a 2 ) +

+ ... V ^{( a} 2 ^{- a} n ) e n V ^{( a} n ) =

= £V(a2 -ai)eA, i=1

V ( a n - a ) eV ( a ) + V ( a n - a 2 )e-V ( a 2 ) +

+ ...(1 + V (0) e n )V ( a n ) =

=j^ V ^{( a} n ^- a) e i A i .

i = 1

Определитель матрицы системы (7) имеет вид: A = det|| bj II, где bj = §ij + V(ai - aj)ej, i, j = 1,...n , 5ij = <

' ¹, i = j

0, i * j

.

Обозначим

A k = detPkA.+ (1 - Aj) bj, где

n ai =E V(ai - aj )emAm .

m = 1

Тогда по формуле Крамера

V ( a * )   Л ■

Имея в виду зависимость A k от A m согласно (7), получаем для V ( a k )

выражение вида:

n

V ( a k ) = Z Pha ,                    (8)

l = 1

det || j(a,-am)e + (1 - 5y)be || где Ры =----------------;----------------.

A

Подставляя (8) в первое уравнение системы (3), получаем [ to ² - p 12 (1 - в 11 )] A 1 + P 2 в 12 A 2 +...............+ P 12 в n A n = 0

P 22 P 21 A + [ to ² - p 2 (1 - P 22 )] A 2 . +..............+ p 2 в n A n = 0

P 2 P n1 A 1 + P 2 P n 2 A 2 ++ [ to ² - P 2 (1 - P nn )] A n = 0

Из условия нетривиальности решения системы (9) относительно амплитуд Al получаем частотное уравнение det|| (^2 - p2) + p2 M= 0, (10) решая которое можно найти собственные частоты to исследуемой механической системы.

В частности, при n = 3 подставляя в (5) последовательно x = a 1 , x = a 2 и x = a 3 , получаем систему линейных алгебраических уравнений относительно V ( a 1 ), V ( a 2 ) и V ( a 3 )

( 1 + V ( 0 ) e 1 ) V ( a 1 ) + V ₂ ( a 1 - a ₂ ) e ₂ V ( a ₂ ) + V 3 ( a 1 - a 3 ) e 3 V ( a 3 ) =

- £ V ( a - a , ) e , A , ,

I = 1

V ( a ₂ - a 1 ) e 1 V ( a 1 ) + ( 1 + V 2 ( 0 ) e ₂ ) V ( a ₂ ) + V 3 ( a 1 - a 3 ) e 3 V ( a 3 ) =

3                                                                         (11)

= E V ( ^a 2 ^- a ) ^e A,

I = 1

V ( a ₂ - a 1 ) e 1 V ( a 1 ) + V ₂ ( a 3 - a ₂ ) e₂V ( a ₂ ) + ( 1 + V 3 ( 0 ) e 3 ) V ( a 3 ) =

= ]Т V , ( a 3 - a , ) e , A,

I = 1

Решая систему (11), получаем

V ^{( a}1 ) = ^в 11 A1 ^{+ в} 12 A 2 ^{+ в} 13 A3,

V ⁽ a 2 ^{) = в} 21 A 1 ^{+ в} 22 A 2 ^{+ в} 23 A 3,           ⁽¹²⁾

^{+ в} 33 A 3,

V ^{( а} 3 ) ^в 31 A 1 ^{+ в} 32 A 2

где

в

¹¹ = А

_ V1 Фк

V 1 ( a ₂ - a 1 ) e 1

V 1 ( a ₃ - a 1 ) e 1

в

¹² = А

V 2 ( a 1 - a ₂ ) e ₂

_ V ( 0 ) e 2

V 2 ( a ₃ - a ₂ ) e ₂

в

¹³ = А

V ₃ ( a 1 - a ₃ ) e ₃

V ₃ ( a ₂ - a ₃ ) e ₃

V ( 0 ) e 3

в

²¹ = А

J + V 1 ( 0 ) e 1

V 1 ( a ₂ - a 1 ) e 1

V 1 ( a ₃ - a 1 ) e 1

в

²² = А

J + V 1 (Q ^

V 1 ( a ₂ - a 1 ) e 1

V 1 ( a ₃ - a 1 ) e 1

в

²³ = А

J + V 1 (Q^

V 1 ( a ₂ - a 1 ) e 1

V 1 ( a ₃ - a 1 ) e 1

V ₂ ( a 1 - a ₂ ) e ₂ 1 + V₂ ( 0 ) e 2

V 2 ( a ₃ - a ₂ ) e ₂

V 2 ( a 1 - a ₂ ) e ₂

1 + V ( 0 ) e 2

V ₂ ( a ₃ - a ₂ ) e ₂

V 2 ( a 1 - a ₂ ) e ₂ 1 + V ( 0 ) e 2

V 2 ( a ₃ - a ₂ ) e ₂

_ V 1 (0> 1

V 1 ( a ₂ - a 1 ) e 1

V 1 ( a ₃ - a 1 ) e 1

V 2 ( a 1 - a ₂ ) e ₂

_ ( 0> 2

V 2 ( a ₃ - a ₂ ) e ₂

V ₃ ( a 1 - a ₃ ) e ₃

V ₃ ( a ₂ - a ₃ ) e ₃

V ( 0 ) e 3

V ₃ ( a 1 - a ₃ ) e ₃

V ₃ ( a ₂ - a ₃ ) e ₃ 1 + V , ( 0 ) e ₃

V ₃ ( a 1 - a ₃ ) e ₃

V ₃ ( a ₂ - a ₃ ) e ₃

1 + V₃ ( 0 ) e ₃

V ₃ ( a 1 - a ₃ ) e ₃

V ₃ ( a ₂ - a ₃ ) e ₃

1 + V₃ ( 0 ) e ₃

V ₃ ( a 1 - a ₃ ) e ₃

V ₃ ( a ₂ - a ₃ ) e ₃ 1 + V₃ ( 0 ) e ₃

V ₃ ( a 1 - a ₃ ) e ₃

V ₃ ( a ₂ - a ₃ ) e ₃ 1 + V₃ ( 0 ) e ₃

V ₃ ( a 1 - a ₃ ) e ₃

V ₃ ( a ₂ - a ₃ ) e ₃ 1 + V₃ ( 0 ) e ₃

^в 4

J + V i (0> 1

V 1 ( a 2 - a 1 ) e 1

V 1 ( a ₃ - a 1 ) e 1

V₂ ( a 1 - a 2 ) e 2       V 1 ( 0 ) e 1

1 + V 2 ( 0 ) e 2     V 1 ( a 2 - a 1 ) e 1

V₂ ( a ₃ - a 2 ) e 2 V 1 ( a ₃ - a 1 ) e 1

^ 32 4

J + V 1 (0> 1

V 1 ( a ₂ - a 1 ) e 1

V 1 ( a ₃ - a 1 ) e 1

V 2 ( a 1 - a ₂ ) e 2

1 + V 2 ( 0 ) e 2

V 2 ( a ₃ - a ₂ ) e 2

V 2 ( a 1 - a ₂ ) e 2

_ V 2 ( 0 ) e 2

V 2 ( a ₃ - a ₂ ) e 2

в

²² = A

J + V 1 (0> 1

V 1 ( a ₂ - a 1 ) e 1

V 1 ( a ₃ - a 1 ) e 1

V 2 ( a 1 - a ₂ ) e 2 V₃ ( a 1 - a ₃ ) e ₃ 1 + V 2 ( 0 ) e 2 V ₃ ( a ₂ - a ₃ ) e ₃

V 2 ( a ₃ - a 2 ) e 2       V ₃ ( 0 ) e ₃

A =

J + V 1 ( 0 ) 6 , V 1 ( a 2 - a 1 ) e 1 V 1 ( a ₃ - a 1 ) e 1

V 2 ( a 1 - a 2 ) e 2 V ₃ ( a 1 - a ₃ ) e ₃

1 + V 2 ( 0 ) e 2     V ₃ ( a 2 - a ₃ ) e ₃

V 2 ( a ₃ - a 2 ) e 2 1 + V ₃ ( 0 ) e ₃

Подставляя (12) в уравнения системы (3), получаем следующую систему:

222  22

^{( p}1 ^- » ^{- p} 1 ^e _n) A1 ^{- p} 1 ^в 12 A 2 ^{- p} 1 ^в 13 A 3 = ⁰,

^{( p} 2 ®    ^p 2 ^в 22 ) A 2 ^p 2 ^в 23 A 3 ⁰,

^- p 2 ^в 21 A 1 ⁺

2   2    222

^p 3 ^в 31 A 1 p 3 ^в 32 A 2 ^{+( p} 3 ®    p 3 ^в 33 ) A 3 ^0.

Отсюда получим уравнение для собственных частот:

P 2 ^- ю ^{2 -} p 1^{2 в} 11

^- p 2 ^в 21

^- p 3 ^в 31

^p 1 ^в 12

p 2 - ю ² - p 2 Д

^{- p} 3 ^в 32

^p 1 ^в 13

^- p 2 ^в 23

^p 3 - ^ю - ^p 3 ^в 33

= 0 .

По найденным собственным частотам ю определяем формы собственных колебаний по формуле (6):

V ( x ) = V ( x - a 1 )e 1 (A ] - V ( a j )) + V , ( x - a ₂) e ₂( A ₂ - V ( a ₂)) + + V 3 ( x - a 3 ) e 3 ( A 3 - V ( a 3 )).

2 Сравнительный анализ

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

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

С целью проведения сравнительного анализа произведен расчет собственных частот балки с тремя телами, то есть при n = 3, опираясь на данные и результаты расчетов из зарубежной статьи [2].

Параметры балки с упруго закрепленными тремя телами в системе СИ: к1 = 3kb,m1 = 0,2mb, k2 = 4,5kb,m2 = 0,5mb, k3 = 6kb,m3 = 1 mb, kb = 63476,1, mb = 15,3875, a1 = 0,1,a2 = 0,4,a3 = 0,8, l = 1,EJ = 63476,1,pF = 15,3875 .

Сравнение частот, полученных в зарубежной статье и в данной работе

	« х	^ 2	Ю 3	® 4	^ 5
Зарубежная статья [2]	156,6703	190,6994	248,6622	1454,2932	3968,4732
Данная работа	156,6705	190,6995	248,6625	1454,3015	3968,4815
Погрешность	1,2740-6	5,2440-7	1,21-10 " 6	5,71 •Ю-6	2,0940-6

Заключение

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

Программа const maxElement = 3;

type

MyFloat = Extended;

TArrayMyFloat = array[0..maxElement, 0..maxElement] of MyFloat;

var

• E, J, ro, A : MyFloat;

• b, beta : MyFloat;

• w, dw, w1, w2 : MyFloat;

M, M0 : TArrayMyFloat;

• L, a1, a2, a3 : MyFloat;

e1, e2, e3 : MyFloat;

k1, k2, k3 : MyFloat;

m1, m2, m3 : MyFloat;

p1, p2, p3 : MyFloat;

c_ij, beta_ij : TArrayMyFloat;

//функции крылова function S1(y : MyFloat) : MyFloat;

begin

Result := (exp(y) + exp(-y) + 2 * cos(y)) / 4; end;

function S2(y : MyFloat) : MyFloat;

begin

Result := (exp(y) - exp(-y) + 2 * sin(y)) / 4; end;

function S3(y : MyFloat) : MyFloat;

begin

Result := (exp(y) + exp(-y) - 2 * cos(y)) / 4; end;

function S4(y : MyFloat) : MyFloat;

begin

Result := (exp(y) - exp(-y) - 2 * sin(y)) / 4; end;

//вычисляем (-1)^i function od(i : Integer) : Integer; begin if ((i + 1) mod 2) = 0 then Result := 1 else Result := -1;

end;

//удаляем строку и столбец из матрицы function Del(A : TArrayMyFloat; l : Integer; m : Integer) : TArrayMyFloat;

var

• j, g : Integer;

begin for g := l to m - 1 do for j := 0 to m do

A[g, j] := A[g + 1, j];

for g := 0 to m - 1 do for j := 0 to m - 1 do

A[g, j] := A[g, j + 1];

Result := A;

end;

//рекурсивно вычисляем детерминант function Det(A : TArrayMyFloat; m : Integer) : My-Float;

Var i : Integer;

buf : MyFloat;

begin buf := 0;

if m = 0 then buf := A[0, 0] else for i := 0 to m do buf := buf + od(i + 1) * A[i,

• 0] * Det(Del(A, i, m), m - 1);

Result:=buf;

end;

//очистка временной таблицы procedure ClearM0;

var

• i,    j : Integer;

begin for i := 0 to maxElement do for j := 0 to maxElement do

M0[i, j] := 0;

end;

//подготовка временной матрицы procedure PrepareM0(a : MyFloat);

begin

ClearM0;

M0[0, 0] := S1(-beta * a); M0[0, 1] := S2(-beta * a); M0[0, 2] := S3(-beta * a); M0[0, 3] := S4(-beta * a); M0[1, 0] := S1(beta * (L - a)); M0[1, 1] := S2(beta * (L - a)); M0[1, 2] := S3(beta * (L - a)); M0[1, 3] := S4(beta * (L - a)); M0[2, 0] := S4(-beta * a); M0[2, 1] := S1(-beta * a); M0[2, 2] := S2(-beta * a); M0[2, 3] := S3(-beta * a); M0[3, 0] := S4(beta * (L - a)); M0[3, 1] := S1(beta * (L - a)); M0[3, 2] := S2(beta * (L - a)); M0[3, 3] := S3(beta * (L - a)); end; function get_ai(i : Integer) : MyFloat; begin if i = 0 then Result :=a1;

if i = 1 then Result :=a2;

if i = 2 then Result :=a3;

end;

//вычисление коэффициентов разложения Ci procedure Calc_c;

var

• a, delta, delta_col : MyFloat;

elem1, elem3 : MyFloat;

row, col : Integer;

begin for row := 0 to 2 do begin a := get_ai(row);

PrepareM0(a);

(L - a)) / (b * beta * beta

(L - a)) / (b * beta * beta

elem1 := -S4(beta *

• *    beta);

elem3 := -S3(beta *

• *    beta);

delta := Det(M0, 3);

for col := 0 to 3 do begin

M := M0;

M[0, col] := 0;

M[1, col] := elem1;

M[2, col] := 0;

M[3, col] := elem3;

delta_col := Det(M, 3);

c_ij[row, col] := delta_col / delta; end;

end; end;

//функция хэвисайда function hevi(x : MyFloat) : MyFloat;

begin if x < 0 then Result := 0 else Result := 1; end;

//вычисление Vi(x)

function Vi(row : Integer; x : MyFloat) : MyFloat; begin

Result := c_ij[row, 0] * S1(beta * x) + c_ij[row, 1]  * S2(beta * x) + c_ij[row, 2]  * S3(beta * x) + c_ij[row, 3] * S4(beta * x) + S4(beta * x) * hevi(x) / (b * beta * beta * beta);

end;

//подготовка временной матрицы procedure PrepareM0_again;

begin

ClearM0;

M0[0, 0] = 1 + Vi(0, 0) * e1; M0[0, 1] = Vi(1, a1 - a2) * e2; M0[0, 2] = Vi(2, a1 - a3) * e3; M0[1, 0] = Vi(0, a2 - a1) * e1; M0[1, 1] = 1 + Vi(1, 0) * e2; M0[1, 2] = Vi(2, a2 - a3) * e3; M0[2, 0] = Vi(0, a3 - a1) * e1; M0[2, 1] = Vi(1, a3 - a2) * e2; M0[2, 2] end; = 1 + Vi(2, 0) * e3; function get_ei(i : Integer) : MyFloat; begin if i = 0 then Result := e1;

if i = 1 then Result := e2;

if i = 2 then Result := e3;

end;

//расчет коэффициентов beta_ij procedure Calc_beta;

var i, row, col : Integer;

delta, delta_colrow : MyFloat;

begin

PrepareM0_again;

delta := Det(M0, 2);

for row := 0 to 2 do for col := 0 to 2 do begin

M := M0;

if row = col then M[row, col] := Vi(row, 0)

* get_ei(row) else for i := 0 to 2 do M[i, row] := Vi(col, get_ai(i) - get_ai(col)) * get_ei(col);

delta_colrow := Det(M, 2);

beta_ij[row, col] := delta_colrow / delta; end;

end;

//основная процедура с циклом по частоте w procedure TMainF.StartBtn1Click(Sender: TObject); var

• f, f0 : MyFloat;

Step : Integer;

begin

//инициализация параметров

E := 2.069E+11;

J := 3.06796E-7;

ro := 7.8367E+3;

A := Pi * sqr(0.05) / 4;

a1 := 0.1;

a2 := 0.4;

a3 := 0.8;

• L := 1;

k1 := 3 * 6.34761E+4;

k2 := 4.5 * 6.34761E+4;

k3 := 6 * 6.34761E+4;

m1 := 0.2 * 15.3875;

m2 := 0.5 * 15.3875;

m3 := 1.0 * 15.3875;

//параметры итератора dw :=0.01;

w1 :=100;

w2 :=4000;

//расчет параметров p1 := sqrt(k1 /m1);

p2 := sqrt(k2 /m2);

p3 := sqrt(k3 /m3);

b := E * J / (ro* A);

e1 := k1 /  (ro *A);

e2 := k2 /  (ro *A);

e3 := k3 /  (ro *A);

//цикл по w с шагом dw w := w1;

f0 := 1E-9;

Step := 0;

repeat beta := sqrt(w) / sqrt(sqrt(b));

//вычисляем с1, с2, с3, с4 для a1, a2, a3

Calc_c;

//вычисляем beta11, beta12 .. beta33

Calc_beta;

//вычисляем f(w)=

M[0,  0]   := sqr(p1) *  (1  - beta_ij[0, 0])  - sqr(w);

M[0, 1] := - sqr(p1) * beta_ij[0, 1];

M[0, 2] := - sqr(p1) * beta_ij[0, 2];

M[1, 0] := -sqr(p2) * beta_ij[1, 0];

M[1,  1]   := sqr(p2) *  (1  - beta_ij[1, 1])  - sqr(w);

M[1, 2] := -sqr(p2) * beta_ij[1, 2];

M[2, 0] := - sqr(p3) * beta_ij[2, 0];

M[2, 1] := - sqr(p3) * beta_ij[2, 1];

M[2,  2]   := sqr(p3) *  (1  - beta_ij[2, 2])  - sqr(w);

f := Det(M, 2);

f0 := f;

• w := w + dw;

Inc(Step);

until w > w2;

end;