Продольные колебания стержня

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

Представим однородный стержень постоянного сечения A массы M в виде системы n материальных точек, соединенных невесомыми линейно-упругими элементами – пружинами. Масса материальных точек:

m_i = M/n

Жесткость пружин:

c_i = n \cdot \frac{EA}{L}

где E – модуль Юнга, L – длина стержня.

Уравнение движения массы i имеет вид (в проекции на горизонтальную ось x )

m_i \ddot{x}_i = -F_i + F_{i+1}, \quad i = 1, \ldots, n-1

Для последней массы (i=n)

m_n \ddot{x}_n = -F_n

Силы, действующие на массы вычисляются следующим образом

F_k = c_k (x_k - x_{k-1} - l_0), \quad k=2,\ldots,n

Для первой точки

F_1 = c_1 (x_1 - l_0)

где l_0 - свободная длина пружины

l_0 = \frac{L}{n}

Image

Программа (файл-скрипт)

% Модуль Юнга Н/м2
E  = 1000000;
% Площадь поперечного сечения стержня м2
A  = 0.001;
% Погонная масса кг/м
mu = 10.0;
% Длина стержня м
L  = 1.0;
% Масса стержня
M   = L*mu;
% Количество материальных точек, на которые разбивается стержень
N   = 30;
% Начальное расстояние между точками (l0)
p.L = L/N;
% Масса материальной точки
p.m = M/N;
% Жесткость пружин
p.c = N*E*A/L;
%
% Начальное положение точек
x0  = (1:N)'/N*(L+0.1);
% Начальная скорость
vx0 = zeros(N,1);
% Вектор состояния
q0  = [x0;vx0];

% Интегрирование
[t, q] = ode113(@(t,q) dqdt(t,q,p),0:0.002:2,q0,odeset('RelTol',1e-8));

% Вектор деформаций пружин
dx = [q(:,1)-p.L (q(:,2:N) - q(:,1:N-1))-p.L];
% Максимальная деформация
dx_max = max(max(dx));
% Минимальная деформация
dx_min = -dx_max;
% Открываем файл для записи видео
v = VideoWriter('rod_d.avi');
open(v);
% Фигура
figure('Position',[100 100 1920 1080]);
axis([0 L*1.3,-0.5, 0.5]);
set(gca,'FontSize',20);
grid on; hold on;
% Палитра
colormap('jet');
% Для каждого момента времени из таблицы результатов интегрирования
for i=1:size(t,1)
    % Очищаем рисунок
    cla;
    % Массив деформаций пружин
    dydx = ([q(i,1) diff(q(i,1:N))] - repmat(p.L,1,N));
    % Вершины прямоугольников, изображающих пружины
    % Каждый столбец содержит 4-ку координат для прямоугольника
    % Координаты x   
    vertex_x  = [0 q(i,1) q(i,1) 0; q(i,1:N-1)' q(i,2:N)' q(i,2:N)' q(i,1:N-1)'];
    % Координаты y
    vertex_y  = repmat([0.1 0.1 -0.1 -0.1],N,1);     
    % Рисуем все прямоугольники сразу
    % Третий аргумент - деформация каждой пружины для раскраски прямоугольника
    patch(vertex_x', vertex_y', dydx,'FaceColor','flat');
    % Приведение диапазона деформаций к диапазону цветов
    caxis([dx_min dx_max]);
    colorbar;
    text(0.1,0.45,sprintf('T=%5.3f c. Длина стержня L=%5.3f м',t(i),q(i,N)),'FontSize',20);   
    % График деформаций
    fplot(@(x) interp1([0 q(i,1:N)],[0 dydx]*0.1/dx_max,x,'next'),[0 q(i,N)],'k-','LineWidth',1);
    frame = getframe(gcf);
    writeVideo(v,frame);
end
close(v);

Файл-функция правых частей дифференциальных уравнений движения системы материальных точек

function dq = dqdt(t, q, p)

n  = size(q,1)/2;
x  = q(1:n);

% Расстояние между точками для определения деформаций пружин
% x1 (x2-x1) (x3-x2) (x4-x3) ... (xn - x_n-1)
xdif = [x(1);diff(x)];

% Из расстояний между точками вычитается свободная длина пружины
% определяется деформация
% x1-l0 (x2-x1)-l0 (x3-x2)-l0 (x4-x3)-l0 ... (xn - x_n-1)-l0
dx   = xdif - repmat(p.L,n,1);

% По деформации определяется сила
% F = [-(x1-l0)*c -(x2-x1-l0)*c ...]
F    = -dx*p.c;
% На все точки, кроме последей, действуют еще силы со стороны соседней пружины справа
% с противоположным знаком. На последнюю точку действует только одна пружина,
% поэтому к действующей на нее силе добавляется ноль
F    = F - [F(2:end); 0];

% Ускорения точек
a    = F/p.m;
% Результат работы функции -- скорости и ускорения точек
dq   = [q(n+1:end);a];

end

Стержень, как распределённая система

Бабаков И.М. Теория колебаний. 4-е изд., испр. - М.: Дрофа, 2004. - 591 с

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

\frac{\partial^2 y}{\partial t^2} - c^2  \frac{\partial^2 y}{\partial x^2} = 0

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

c^2 = \frac{E \cdot A}{\mu}

В последней формуле плотность материала стержня определяется отношением погонной массы \mu к площади поперечного сечения A.

Сложные малые колебания рассматриваемой системы с бесконечным числом степеней свободы представляются в виде суммы гармонических колебаний с частотами p_i и формами колебаний d_i(x), определяющих распределение деформаций по длине стержня

y_i(x,t) = d_i(x) \sin (p_i \cdot t + \alpha_i)

где d_i(x)i-я форма главного колебания – амлитудное значение деформации, как функция положения точки на стержне.

Подстановка y(x,t) = d(x) \sin (p \cdot t + \alpha) в уравнение движения приводит к уравнению собственных форм колебаний:

d''(x) + a^2 d(x) = 0

где

a^2 = \frac{p^2 \mu}{EA}

Решение полученного однородного дифференциального уравнения второго порядка имеет следующий вид:

d_i(x) = B \cos a_i x + C \sin a_i x

Постоянные B и C определяются из граничных условий – условий закрепления стержня. Например для закрпеленного с одной стороны стержня граничные условия будут иметь вид:

d(0) = 0
d'(L) = 0

Подставляя эти граничные условия в решение, получим

B = 0

и

C a_i \cos a_i L = 0

Нетривиальное решение этого уравнения имеет вид:

\cos a_i L  = 0 \quad a_i = \frac{1+2k}{2 L} \pi, \quad k=0,1,2,\ldots,\infty

Таким образом, ряд сообственных частот продольных колебаний стержня c закреплённым (x=0) будет определяться следующим образом:

p_k = \sqrt{\frac{EA}{\mu}} \frac{1+2k}{2 L} \pi, \quad k=0,1,2,\ldots,\infty

Общее решение имеет следующий вид

y(x, t) = \sum_{i=1}^{\infty} d_i(x) (M_i \cos p_i t + N_i \sin p_i t)

Постоянные M_i, N_i определяются из начальных условий – распределения начальной деформации по длиней стержня:

y(x,0) = f(x)

и распределением начальных скоростей:

\dot{y}(x,0) = \gamma(x)

В общем случае константы M_i, N_i определяются следующим образом:

M_i = \int_0^L d_i(x) f(x) dx
N_i = \frac{1}{p_i} \int_0^L d_i(x)\gamma(x) dx
% Модуль Юнга Н/м2
E  = 1000000;
% Площадь поперечного сечения м2
A  = 0.001;
% Погонная масса kg/м
mu = 10.0;
% Длина стержная м
L  = 1.0;
% Количество учитываемых собственных форм
k  = 1:10;
a  = (2*k-1)*pi*0.5/L;
% Массив частот
p  = a*sqrt(E*A/mu); % рад/с
% Начальная деформация
% К концу стержня приложена растягивающая сила F
F  = 100; % Н
y0 = @(xx) interp1(linspace(0,L,10),F*linspace(0,L,10)/(E*A),xx);
% Постоянные, определяемые по начальным условиям (для частоты i)
M  = @(i)  (2/L)*integral(@(xx) y0(xx).*sin(a(i).*xx),0, L);
M  = @(i)  (2*F/(L*E*A))*sin(a(i)*L)/(a(i)^2.0);
% Начальная скорость деформаций равна нулю, поэтому (для частоты i)
N  = @(i) 0;
% Слагаемые общего решения
yk = @(i,t,x) (M(i)*cos(p(i)*t)+N(i)*sin(p(i)*t))*sin((2*i-1)*pi*x*0.5/L);
% Слагаемые деформации
ek = @(i,t,x) (M(i)*cos(p(i)*t)+N(i)*sin(p(i)*t))*cos((2*i-1)*pi*x*0.5/L)*(2*i-1)*pi*0.5/L;
% Общее решение y(t,x)
y   = @(t,x) sum(arrayfun(@(i) yk(i,t,x),k));
% Деформация
eks = @(t,x) sum(arrayfun(@(i) ek(i,t,x),k));
% ------------------------------------------------------------------------
% Рисуем
% ------------------------------------------------------------------------
% Видео
v = VideoWriter('rod_c.avi');
open(v);

% Разделяем балку на 30 частей (для раскраски)
x = linspace(0,L,30);
figure('Position',[100 100 900 600]);
axis([0 1.3 -0.5 0.5]);
hold on;
% Цветовая палитра JET
cmap = colormap(jet(128));
% Максимальная деформация
maxd =  0.2;
mind = -0.2;
def2index = @(d) floor((d-mind)/((maxd-mind)/size(cmap,1)));
for t=0:0.002:2
    cla;
    yx = arrayfun(@(xx) y(t,xx),x);    
    pos = x + yx;    
    defx = arrayfun(@(xx) eks(t,xx),x);
    for j = 1:size(x,2)
        col = cmap( def2index(defx(j)),:);
        prev_point = 0;
        if j ~= 1
            prev_point = pos(j-1);
        end
        patch([prev_point pos(j) pos(j) prev_point],...
              [0.1 0.1 -0.1 -0.1],col,'EdgeColor','none');            
    end
    plot(pos, arrayfun(@(xx) eks(t,xx),x));
    xlabel('x, м');ylabel('y, м');box('on');

    text(0.1,0.45,sprintf('T=%5.3f c. Длина стержня L=%5.3f м',t,yx(end)+L));    
    text(0.1,0.35,sprintf('Частоты (Гц): '));    
    text(0.1,0.31,sprintf('%3.1f | ',(2*pi./p).^-1));    

    frame = getframe;    
    writeVideo(v,frame);    
end
close(v);

© 2018. All rights reserved.