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

Свободные колебания балки при наличии вязкого трения и с учетом инерции вращения описываются однородным дифференциальным уравнением в частных производных гиперболического типа [7–9], которое для балки переменного сечения из однородного материала примет вид1

(Bu″)″ – Pu″ – rü″ + mü + εmu̇ = 0,

(x, t) ∈ Q ≡ [ (x, t): x ∈ L ≡ (0, l), t > 0 ], (1)

где

B(x) = EJ(х), r(x) = ρJ(х), m(x) = ρF(х),

P – осевая растягивающая сила;

E, ρ – модуль упругости и плотность материала;

F, J – площадь и осевой момент инерции поперечного сечения балки;

ε – коэффициент удельного линейно-вязкого трения.

Точки соответствуют дифференцированию по времени t, штрихи – по пространственной координате x.

В уравнении (1) в порядке следования слагаемых учтены силы упругости, осевая сила, силы инерции вращения, инерционная сила линейных перемещений, силы линейно-вязкого трения [10–12].

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

u(0,t) = u″(0,t) = u(1,t) = u″(1,t). (2)

Решение задачи (1), (2) отыскивается с помощью метода разделения переменных как произведение

u(x,t) = X(x)eλt. (3)

Здесь характеристический показатель

λ = – μ + jω, (4)

где j – мнимая единица;

μ и ω – подлежащие определению коэффициент затухания и частота свободных колебаний.

Подстановка (3) в (1), (2) дает

(BX″) – PX″ – rλ2X″ +(mλ2 + εmλ)X, (5) X(0) = X″(0) = X(1) = X″(1) = 0. (6) Введем обозначения g(x) = P + rλ2, s(x) = mλ2 + εmλ

и перепишем уравнение (5)

(BX″)″ – gX″ + sX = 0. (7)

Далее воспользуемся методом конечных разностей. С этой целью вместо области L + Г(Г – граничные точки) введем дискретную область lh в виде множества узлов равномерной сетки с шагом h

lh ≡ [ xi: xi = (i – 1)h, i = 1, 2, …, n, N; N = n + 1], h = l/n,

где N – количество узлов сетки.

Значения функций и производных в уравнении (7) заменим приближенными или точными (если возможно)

X(xi) ≈ yi,

X″(xi) ≈,

g(xi) = gi, s(xi) = si

и получим на первом этапе равенство

. (8)

Повторное применение процедуры замены второй производной конечноразностными соотношениями и приведение подобных членов дает однородную систему алгебраических уравнений

ai yi–2 – bi yi–1 + ci yi – di yi+1+ ei yi+2 = 0, i = 4, 5, …, n – 2, (9)

где ai = Bi–1, bi = 2Bi–1 + 2Bi + gih2, ci = Bi–1 + 4Bi + Bi+1 + 2gih2 + sih4, di = 2Bi + 2Bi+1 + gih2, ei = Bi+1.

Аналогичные замены проведем в граничных условиях (6), воспользуемся уравнением (9) для точек i = 2, 3, n – 1, n и получим четыре уравнения

q2 y2 – d2 y3 + e2 y4 = 0,

– b3 y2 + c3 y3 – d3 y4 + e3 y5 = 0,

an–1 yn–3 – bn–1 yn–2 + cn–1 yn–1 – dn–1 yn = 0,

an yn–2 – bn yn–1 + pn yn = 0.

Здесь обозначено

q2 = c2 – a2,

pn = cn – en.

Появляющиеся при этом значения функции во внесеточных узлах y0, yn+2 выражаются через значения функции в узлах сетки. Так, например, использование граничного условия на левом конце

X″(0) = 0

при i = 1 приводит к соотношению

из которого при y1 = 0 легко получается, что

y0 = – y2.

Аналогично граничные условия на правом конце дают

yn+2 = – yn.

Исключим с помощью данных соотношений y0 и yn+2 из числа компонентов вектора y, а также равные нулю y1 и yN и получим вектор yT = {y2, y3, …, yn–1, yn}, компонентами которого являются отклонения балки. Тогда систему уравнений можно переписать в матрично-векторной форме

G(λ) y = 0, (10)

где G(λ) – квадратная матрица порядка n – 1. Можно легко показать, что матрица коэффициентов является ленточной и пятидиагональной вида

Здесь нулевые элементы не выписаны.

Элементы матрицы G являются функциями характеристического показателя λ и через него – коэффициента затухания μ и частоты колебаний ω в соответствии с (4).

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

det G(λ) = 0, (11)

из которого определяются собственные значения, т. е. спектр {λ1, λ2, …}.

Левая часть уравнения (11) является полиномом порядка 2(n – 1), само же уравнение тогда оказывается алгебраическим вида f(λ) = 0. При больших значениях n его написание в развернутой форме, хотя и возможно, представляет громоздкую процедуру. Кроме того, получить его аналитические решения удается лишь в простейших случаях. Выход из таких затруднений состоит в применении численных методов и ЭВМ. При этом будем ориентироваться на алгоритмические языки и программные системы, позволяющие непосредственно пользоваться функциями комплексной переменной и проводить алгебраические и другие действия над ними (например, MathCad, C++, MATLAB [13–15] и т. д.).

Уравнение (11) с учетом того, что его левая часть представляет комплексную функцию, можно переписать в виде

f1(μ, ω) + j f2(μ, ω) = 0.

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

f1(μ, ω) = 0, f2(μ, ω) = 0. (12)

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

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

Решение системы уравнений найдем с помощью метода покоординатного спуска [16–18]. С этой целью образуем вспомогательную неотрицательную функцию

Φ(μ, ω) = | det (μ, ω) | = [f 12(μ, ω) + f 22(μ, ω)]1/2, обращающуюся в нуль лишь в том случае, если f1(μ, ω) = f2(μ, ω) = 0. Таким образом, решение исходной системы (12) будет одновременно точкой нулевого минимума скалярной функции векторного аргумента Ф(μ, ω), так что теперь необходимо решать более простую задачу Ф(z) = 0, z = {μ, ω} ∈ R2 ∩ (μ > 0, ω > 0).

Здесь вектор z соответствует искомым величинам, Ф(z) представляет мультимодальную функцию, имеющую множество локальных минимумов.

Суть применяемого метода состоит в построении последовательности точек (приближений к решению) zk, k = 0, 1, ..., сходящейся к точке локального минимума z*. При этом в процессе вычислений необходимо добиваться, чтобы значения вспомогательной функции были монотонно убывающими и ограниченными снизу

Ф(z0) ≥ Ф(z1) ≥.... ≥ Ф(zk) ≥.... ≥ Ф(z*). (13)

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

zk+1 = zk + Δk,

где Δk – вектор приращений координат, обеспечивающий выполнение условий (13).

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

Ф(zk+1) < δ, (14)

причем δ – априорно задаваемое малое положительное число (точность вычислений).

Многократное повторение такой процедуры вычислений дает спектры коэффициентов затухания и собственных частот

{(μ1, ω1), (μ2, ω2),…}.

При реализации предлагаемой численной схемы наибольшая сложность состоит в правильном выборе приращений координат Δ = {Δμ, Δω} на каждом шаге по величине и по знаку. В этом месте обычные подходы к составлению компьютерных программ вычислений требуют высокого программистского искусства, а сами программы получаются сложными и громоздкими. Трудности существенно упрощаются, если используемая система программирования позволяет оперативно визуализировать результаты вычислений. Поясним сказанное на примере данной задачи при пользовании программной системой MATLAB. Один из аргументов функции Φ(μ, ω), например μ, фиксируется, вычисляются ординаты функции Φ(ω | μ), на экран компьютера выводится ее график (рис. 1а), вычисления приостанавливаются. По кривой, являющейся сечением поверхности Φ(μ, ω), отчетливо видны условные минимумы функции Φ(ω | μ). Теперь фиксируем аргумент ω на значении, соответствующем условному минимуму и строим новый график Φ(μ | ω) (рис. 1б) и т. д. Такая процедура продолжается до выполнения условия (14).

Пример 1. В качестве тестовой проверки рассмотрим балку постоянного сечения из двутавра № 14 с численными данными

l = 4 м, n = 400, ρ = 7800 кг/м3, E = 200 ГПа, F = 17,4 см2, J = 572 см4, ε = 0,01 с-1, P = 0.

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

В табл. 2 показаны значения пяти элементов спектра собственных частот для балки с такими же численными данными при Р = –0,25Рэ (Рэ – критическая эйлерова сила).

Примечательно, что для получения результатов требуется небольшое количество шагов в методе покоординатного спуска (4–5 для каждого k).

Тот же спектр (при Р = 0), полученный точными аналитическими методами, приведен в табл. 3.

Ниже в табл. 4 показана разница между собственными частотами колебаний балки, полученными точными аналитическими методами и методом конечных разностей.

Видно, что разница небольшая, не имеющая практической значимости.

Результаты данного примера убедительно показывают, что метод конечных разностей дает надежные результаты, позволяющие универсальным и простым способом определять динамические характеристики свободных колебаний балок переменного сечения [19–20].

Спектр собственных форм может быть установлен с помощью системы уравнений (10). Ввиду того, что определитель матрицы G равен нулю, собственные формы могут быть найдены лишь с точностью до сомножителя. Тогда в (10) можно принять, например, что y2 = 1, и остальные неизвестные найти из системы уравнений, из которой исключена последняя строка (см. на стр. 16).

Заметим, что правая часть системы уравнений образована из элементов первого столбца матрицы G.

Пример 2. По данным примера 1 вычислим ординаты трех первых собственных форм φk(х) = yk(x).

Результаты счета в виде графиков показаны на рис. 2. Номера кривых при этом совпадают с номерами собственных частот; проведена нормировка форм, удовлетворяющая условию

| yi | = 1.

Данные кривые почти неразличимо совпадают с собственными формами, найденными с помощью точных аналитических методов.

Пример 3. Рассмотрена балка из сварного двутавра длиной 4 м, состоящая из трех участков, поперечное сечение которых изменяется по длине (рис. 3), с численными данными ρ = 7800 кг/м3, E = 200 ГПа, n = 500, ε = 0,01 c-1, l1 = l3 = 1 м, F1 = F3 = 61,9 см2,J1 = J3 = 13 380 см4, l2 = 2 м, F2 = 84,7 см2, J2 = 27 696 см4.

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

В табл. 6 представлены те же пять элементов спектров коэффициентов затухания и собственных частот для классической модели.

В табл. 7 показана разница между собственными частотами балки по двум моделям.

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

1. В данной модели колебаний деформации сдвига не учитываются.