Динамика ротора (Maxima) (Материалы Панкратьевой по ротору для курсовой работы)
Описание файла
Файл "Динамика ротора (Maxima)" внутри архива находится в папке "2". Документ из архива "Материалы Панкратьевой по ротору для курсовой работы", который расположен в категории "". Всё это находится в предмете "теоретическая механика" из 3 семестр, которые можно найти в файловом архиве НИУ «МЭИ» . Не смотря на прямую связь этого архива с НИУ «МЭИ» , его также можно найти и в других разделах. .
Онлайн просмотр документа "Динамика ротора (Maxima)"
Текст из документа "Динамика ротора (Maxima)"
Типовой расчёт. Динамические реакции в подшипниках
ротора.
Пример выполнения задания.
Введение параметров
--> b:0.267;
(%o1)
--> M0:1950;
k1:1.57;
k2:0.039;
m1:47;
m2:57;
m3:37;
R1:0.17;
R2:0.2;
R3:0.15;
O1C1:4·10^−4;
O3C3:3·10^−4;
alpha:0.019;
beta:0.68;
a:0.133;
zA:0;
zB:b;
c:0.093;
phi0:0;
omega_z0:600;
tau:0.48;
delta_t:0.02;
(%o2)
(%o3)
(%o4)
(%o5)
(%o6)
(%o7)
(%o8)
(%o9)
(%o10)
(%o11)
(%o12)
(%o13)
(%o14)
(%o15)
(%o16)
(%o17)
(%o18)
(%o19)
(%o20)
(%o21)
(%o22)
--> MDz:M0−k1·diff(phi(t),t);
MCz:−k2·diff(phi(t),t)^2;
(%o23)
(%o24)
--> Tfin:tau;
(%o25)
--> M:m1+m2+m3;
rC1:matrix([O1C1·sin(beta)],[−O1C1·cos(beta)],[a]);
rC2:matrix([0],[0],[b+c]);
rC3:matrix([0],[O3C3],[−c]);
xC:(m1·rC1[1,1]+m2·rC2[1,1]+m3·rC3[1,1])/M;
yC:(m1·rC1[2,1]+m2·rC2[2,1]+m3·rC3[2,1])/M;
(%o26)
(%o27)
(%o28)
(%o29)
(%o30)
(%o31)
--> I1C1x1y1z1: matrix([m1·R1^2/4,0,0],[0,m1·R1^2/4,0],[0,0,(m1·R1^2/4)·2]);
I2C2XiEtaZeta: matrix([m2·R2^2/4,0,0],[0,m2·R2^2/4,0],[0,0,(m2·R2^2/4)·2]);
I3C3x3y3z3: matrix([m3·R3^2/4,0,0],[0,m3·R3^2/4,0],[0,0,(m3·R3^2/4)·2]);
(%o32)
(%o33)
(%o34)
--> S: matrix([cos(alpha),0,sin(alpha)],[0,1,0],[−sin(alpha),0,cos(alpha)]);
(%o35)
--> ST:transpose(S)$
I2C2x2y2z2:ST.I2C2XiEtaZeta.S;
(%o37)
--> rC1T:transpose(rC1)$
rC2T:transpose(rC2)$
rC3T:transpose(rC3)$
/·Единичная матрица·/
E:matrix([1,0,0],[0,1,0],[0,0,1])$
I1Axyz:I1C1x1y1z1 + m1·((rC1T.rC1).E − rC1.rC1T);
I2Axyz:I2C2x2y2z2 + m2·((rC2T.rC2).E − rC2.rC2T);
I3Axyz:I3C3x3y3z3 + m3·((rC3T.rC3).E − rC3.rC3T);
(%o42)
(%o43)
(%o44)
--> IAxyz:I1Axyz + I2Axyz + I3Axyz;
(%o45)
Составление уравнений для определения динамических реакций
--> eqMz: IAxyz[3,3]·diff(phi(t),t,2) = MDz + MCz;
eqQx: −M·yC·diff(phi(t),t,2) − M·xC·diff(phi(t),t)^2 = XA + XB;
eqQy: M·xC·diff(phi(t),t,2) − M·yC·diff(phi(t),t)^2 = YA + YB;
eqMx: IAxyz[1,3]·diff(phi(t),t,2) − IAxyz[2,3]·diff(phi(t),t)^2 = −zA·YA − zB·YB;
eqMy: IAxyz[2,3]·diff(phi(t),t,2) + IAxyz[1,3]·diff(phi(t),t)^2 = zA·XA + zB·XB;
(%o46)
(%o47)
(%o48)
(%o49)
(%o50)
Решение уравнений для определения динамических реакций и выражение реакций
--> ratprint: false$
solR:solve([eqQx, eqQy, eqMx, eqMy],[XA, XB, YA, YB])$
expand(solR),numer;
(%o53)
Нахождение модулей RA,RB реакций в подшипниках А,B
--> XA:rhs(expand(solR[1][1]))$
XB:rhs(expand(solR[1][2]))$
YA:rhs(expand(solR[1][3]))$
YB:rhs(expand(solR[1][4]))$
--> RA(XA,YA):=sqrt(XA^2 + YA^2);
RB(XB,YB):=sqrt(XB^2 + YB^2);
YB,numer;
(%o58)
(%o59)
(%o60)
Численное решение дифференциального уравнения вращения ротора методом Рунге-Кутта.
Для обращения к rk,стандартной функции Maxima,реализующей этот метод, необходимо
записать уравнение в виде системы двух дифференциальных уравнений первого порядка в нормальной форме Коши.
Необходимо произвести замену переменных в виде
diff(phi(t),t) = omegaz
Тогда вторая производная будет равна
diff(phi(t),t,2) = diff(omegaz,t)
eqMz представим в виде
diff(omegaz,t) = (MDz + MCz)/IAxyz[3,3]
или
diff(omegaz,t) = (M0-k1*omegaz - k2*omegaz^2)/IAxyz[3,3]
Получаем систему из 2х уравнений, которую решаем методом Рунге-Кутта
diff(omegaz,t) = (M0-k1*omegaz - k2*omegaz^2)/IAxyz[3,3]
diff(phi(t),t) = omegaz
Начальные условия
phi(0) = 0
omegaz(0) = omega_z0
--> eq1:(M0−k1·omegaz − k2·omegaz^2)/IAxyz[3,3]$
sol:rk([eq1,omegaz],[omegaz,phi],[omega_z0,0],[t,0,Tfin, delta_t])$
Решение sol является списком значений
Каждое значение состоит из трех компонент
[t, omegaz(t), phi (t)]
Выводим их в отдельные списки для удобства
--> points:Tfin/delta_t+1$
t_list:makelist(sol[k][1],k,1,points)$
diff_phi_list: makelist(sol[k][2],k,1,points)$
phi_list: makelist(sol[k][3],k,1,points)$
Находим значение второй производной phi от t из уравнения
diff(omegaz,t) = (M0-k1*omegaz - k2*omegaz^2)/IAxyz[3,3]
--> diff2_phi_list: makelist(((M0−k1·diff_phi_list[k] − k2·diff_phi_list[k]^2)/IAxyz[3,3]),k,1,points)$
Находим значения реакций, подставляя полученные значения
первой и второй производной phi по t в соответствующие
уравнения
Производим замену производных на переменные для подстановки
--> XA(dphi1,dphi2):=subst([diff(phi(t),t)=dphi1,diff(phi(t),t,2)=dphi2],XA)$
XA_list:makelist(XA(diff_phi_list[k],diff2_phi_list[k]),k,1,points)$
XB(dphi1,dphi2):=subst([diff(phi(t),t)=dphi1,diff(phi(t),t,2)=dphi2],XB)$
XB_list:makelist(XB(diff_phi_list[k],diff2_phi_list[k]),k,1,points)$
YA(dphi1,dphi2):=subst([diff(phi(t),t)=dphi1,diff(phi(t),t,2)=dphi2],YA)$
YA_list:makelist(YA(diff_phi_list[k],diff2_phi_list[k]),k,1,points)$
YB(dphi1,dphi2):=subst([diff(phi(t),t)=dphi1,diff(phi(t),t,2)=dphi2],YB)$
YB_list:makelist(YB(diff_phi_list[k],diff2_phi_list[k]),k,1,points)$
RA_list:makelist(RA(XA_list[k],YA_list[k]),k,1,points),numer$
RB_list:makelist(RB(XB_list[k],YB_list[k]),k,1,points),numer$
Построение графиков функций по результатам интегрирования
--> wxplot2d([discrete,t_list,phi_list],[xlabel, t],[ylabel,phi])$
(%t78)
--> wxplot2d([discrete,t_list,diff_phi_list],[xlabel, t],[ylabel,omega_z])$
(%t79)
-->
wxplot2d([discrete,t_list,diff2_phi_list],[xlabel, t],[ylabel,epsilon_z])$
(%t80)
--> wxplot2d([discrete,t_list,RA_list],[xlabel, t],[ylabel,"RA(t)"])$
(%t81)
Построение таблицы значений функций по результатам интегрирования
--> N:["t"," phi"," omegaz"," epsilonz", " XA"," YA"," XB"," YB"," RA"," RB"]$
T1:matrix(
t_list,
phi_list,
diff_phi_list,
diff2_phi_list,
XA_list,
YA_list,
XB_list,
YB_list,
RA_list,
RB_list
)$
Table:transpose(T1)$
printf(true," ~10a~11a~9a~15a~11a~11a~11a~11a~11a~11a ~% ~{ ~{ ~9,2f ~} ~% ~}",
N[1],N[2],N[3],N[4],N[5],N[6],N[7],N[8],N[9],N[10],
args(Table))$
Контроль решения
--> /·Замена дифференциалов на переменные для подстановки ·/
eqMz(dphi1,dphi2):=subst([diff(phi(t),t)=dphi1,diff(phi(t),t,2)=dphi2],eqMz)$
eqMz(omega,0);
(%o87)
--> solve(eqMz(omega,0),omega),numer;
(%o88)