 |
Моделирование свободных колебаний цепочки связанных гармонических осцилляторов
|
|
|
|
Лабораторная работа №5
Модели, представляющие собой линейные цепочки (рисунке 1), состоящие из конечного или бесконечного числа связанных осцилляторов, оказались весьма эффективными и в настоящее время используются в различных областях физики: физике твёрдого тела, физике сплошных сред, химической физике, радиофизике и др. Используя модели линейных цепочек, оказывается возможным естественным образом осуществить переход к волновым процессам и ввести такие понятия, как длина волны, групповая скорость, фазовая скорость, дисперсия и др.
Отмеченные обстоятельства определяют целесообразность рассмотрения данных моделей в соответствующих курсах физики и компьютерного моделирования. Однако необходимо отметить два важных обстоятельства. Во первых, аналитические решения уравнений движения длинных линейных цепочек (
) могут быть получены только для относительно небольшого числа случаев.
Запишем уравнения движения для каждой массы колебательной системы, представленной на рисунке 1:
Рисунок 1
(1)
Для удобства дальнейшего решения запишем уравнение (1), введя обозначение 
, в следующем виде:
(2)
Следуя общему подходу к решению рассматриваемой задачи, ищем решение системы дифференциальных уравнений в виде:
(3)
Подставив (3) в систему (2), сгруппировав члены, пропорциональные 
и записав систему в матричном виде, получим:
(4)
где
– трёхдиагональная матрица, элементы которой вычисляются по следующим правилам:
(5)
Необходимым и достаточным условием существования решения системы уравнений (4) является равенство нулю определителя матрицы 
(6)
Уравнение (6), называемое характеристическим уравнением, является уравнением степени 
относительно 
. Оно имеет в общем случае различных вещественных положительных корней 
Каждому собственному числу 
соответствует собственный вектор 
, являющийся решением уравнения
|
|
|
(7)
где 
– трёхдиагональная матрица элементы которой вычисляются по следующим правилам:
(8)
Частоту 
называют частотой нормальных колебаний, а вектор 
– вектором нормального колебания, отвечающего 
й частоте. Вектор нормального колебания 
меняется во времени по закону
(9)
Общее решение системы дифференциальных уравнений (2) 
, есть суперпозиция всех векторов нормальных колебаний 
:
(10)
где 
– произвольные постоянные, определяемые из начальных условий.
Скорость движения масс можно определить, продифференцировав (10) по времени:
(11)
Для решения задачи Коши системы дифференциальных уравнений (2) необходимо задать значения координат 
и скоростей 
каждого тела системы в начальный момент времени и решить систему уравнений
(12)
относительно неизвестных .
Запишем (12) в матричном виде
(13)
где 
– нулевая матрица, размерности 
(14)
(15)
(16)
(17)
Система уравнений (13) оказывается нелинейной, однако блочная структура матрицы позволяет найти решение данной системы, не прибегая к численным методам. Для этого, сначала, решим две линейные системы уравнений
(18)
(19)
найдём векторы 
, затем координаты вектора 
(20)
и далее значения фаз каждого нормального колебания:
(21)
Отметим, что функция 
на интервале 
является двузначной (рисунок 2), поэтому для выбора правильных значений данной функции необходимо контролировать знаки числителя и знаменателя дроби в выражении (21). Как видно из рисунка 2, правильное значение угла выбирается по следующим правилам:
Рисунок 2
. (22)
Предваряя описание решения задачи о колебаниях цепочки связанных осцилляторов, приведём алгоритм её решения:
1. Задать число тел, образующих цепочку 
|
|
|
2. Задать массы тел 
3. Задать значения коэффициентов жёсткости пружин 
4. Вычислить элементы матрицы 
в соответствии с (8).
5. Найти собственные числа 
матрицы .
6. Найти собственные векторы 
, соответствующие набору собственных частот .
7. Задать начальные условия 
8. Решить систему линейных уравнений (18), (19) относительно векторов 
и 
соответственно.
9. Вычислить координаты вектора 
в соответствии с (20).
10. Вычислить значения фаз нормальных колебаний 
в соответствии с (21).
11. Определить законы движения тела, образующих колебательную систему в соответствии с (10) и (11).
12. Провести анализ полученных законов движения.
Данный алгоритм в пакете FreeMat реализуется следующей последовательностью действий:
>> clear all
>> N=3; % число тел колебательной системы
>> m=[1 2 1]; % массы тел колебательной системы
>> k=[1 1 1 1]; % жёсткости пружин колебательной системы
>> R0=[-0.2 0 -0.3]; % смещение тел в момент времени 
>> v0=[1 -3 0]; % скорости тел в момент времени
% вычисление элементов матрицы 
>> for alpha=1:N+1
for beta=1:N
omega(alpha,beta)=k(alpha)/m(beta);
end;
end;
>> i=1:N;
>> j=1:N;
% вычисление элементов матрицы OMEGA в
>> OMEGA(i,j)=0;
>> for i=1:N
if i==1
OMEGA(i,i)=omega(1,1)+omega(2,1);
OMEGA(1,2)=-omega(2,1);
end;
if i>1
if i<N
OMEGA(i,i-1)=-omega(i,i);
OMEGA(i,i)=omega(i,i)+omega(i+1,i);
OMEGA(i,i+1)=-omega(i+1,i);
Else
OMEGA(i,i-1)=-omega(i,i);
OMEGA(i,i)=omega(i,i)+omega(i+1,i);
end;
end;
end;
>> [Sigma,Teta]=eig(OMEGA); %вычисление собственных значений и собственных векторов матрицы OMEGA
>> Teta= Teta^0.5; %вычисление собственных частот
>> for i=1:N
for j=1:N
SigmaV(j,i)=-Teta(i,i)*Sigma(j,i);
end;
end;
>> C1=Sigma^-1*R0’; % решение системы уравнений (18)
>> C2=SigmaV^-1*v0’; % решение системы уравнений (19)
>> C=(C1.^2+ C2.^2).^0.5; % вычисление координат вектора С
>> clear alpha
% вычисление фазы нормальных колебаний в соответствии с (21), (22 )
>> for i=1:N
if C(i)==0
alpha(i)=0;
Else
alpha(i)=atan(C2(i)./C1(i));
if C1(i)<0
alpha(i)=pi+alpha(i);
end;
if C1(i)>0
if C2(i)<0
alpha(i)=2*pi+alpha(i);
end;
end;
end;
end;
>> N=length(OMEGA);
>> N1=2^13; % число узлов временной сетки
>> j=1:N1;
>> Tmax=80; % правая граница временного интервала
>> t(j)=(j-1)/(N1-1)*Tmax; % координаты узлов временной сетки
% вычисление значений координат тел в узлах временной сетки
>> for j=1:N1
s=zeros(N,1);
for i=1:N
s=s+C(i)*Sigma(:,i).*cos(Teta(i,i)*t(j)+alpha(i));
end;
X(:,j)=s;
end;
% вычисление значений скоростей тел в узлах временной сетки
|
|
|
>> for j=1:N1
s=zeros(N,1);
for i=1:N
s=s+C(i)*Sigma(:,i).*sin(Teta(i,i)*t(j)+alpha(i));
end;
Xv(:,j)=-s;
end;
% визуализация зависимостей мгновенных значений смещений скорости от времени
>> figure(1); plot(t,X(1,:),’-k’,t,X(2,:),’--k’,t,X(3,:),’:k’)
>> figure(2); plot(t,Xv(1,:),’-k’,t,Xv(2,:),’--k’,t,Xv(3,:),’:k’)
% построение траектории движения тел на фазовой плоскости
>> figure(3); plot(t,X(1,:),Xv(1,:))
>> figure(4); plot(t,X(2,:),Xv(2,:))
>> figure(5); plot(t,X(3,:),Xv(3,:))
% вычисление спектров зависимостей смещений тел системы от времени
>> с1=fft(X(1,:));
>> с2=fft(X(2,:));
>> с3=fft(X(3,:));
>> j=2:N1/2;
% вычисление спектральной плотности смещений тел
>> Cm1(j-1)=abs(c1(j-1))/(N1/2);
>> Cm2(j-1)=abs(c2(j-1))/(N1/2);
>> Cm3(j-1)=abs(c3(j-1))/(N1/2);
>> Freq(j-1)=(j-1)/Tmax; % вычисление частот спектральных гармоник
% визуализация спектральных плотностей смещений тел
>> figure(3); semilogy(Freq,Cm1,’-k’,Freq,10*Cm2,’--k’, Freq,500*Cm3,’:k’)
>>axes([0 2.5 10^-3 2000])
|
|
|
