Matlab  |  Mathcad  |  Maple  |  Mathematica  |  Statistica  |  Другие пакеты Поиск по сайту
Internet-класс  |  Примеры  |  Методики  |  Банк задач  |  Консультации & Форум  |  Download  |  Ссылки  |  Конкурсы
Научно-практический журнал "Exponenta Pro. Математика в приложениях". Вышел 1/2004 номер журнала
Карта сайта | На первую страницу | Поиск |О проекте |Сотрудничество |
Exponenta Pro | Matlab.ru

Наши баннеры


Copyright © 2000-2003. Компания SoftLine. Все права защищены.

Дата последнего обновления информации на сайте: 11.05.04
Сайт начал работу 1.09.00

Программное обеспечение Microsoft, Macromedia, VERITAS, Novell, Borland, Symantec, Oracle и др.

 

 
Условия и точки стабильности порядка 1/2 механической системы
выполнил: Чан Туан Чунг,
Московский энергетический институт(ТУ), 2 курс
2004

archive.gif (75 bytes) Архив разработки (115 Кб, Maple)


Решение задачи с одной обобщённой координатой

Данные значения

>    restart; Digits:=5:

>    MDz:=Mo-k*omega3:
Mnz:=-mu*omega4:
Fnx:=-nu*v5x:

Начальные условия

>    Yslo1:=m2=16,m3=34,m4=26,nu=8000,mu=13,k=15,Mo=10,R1=0.38,I1=9,r1=0.27,R3=0.28,
R4=0.2,r4=0.12,i4=0.15:
Yslo2:=m2=16,m3=34,m4=26,nu=8000,mu=13,k=k,Mo=10,R1=0.38,I1=9,r1=0.27,R3=0.28,
R4=0.2,r4=0.12,i4=0.15:
phi1o:=1.3:omega1o:=0.5:
NaYs:=phi1(0)=phi1o,D(phi1)(0)=omega1o:

Имеем

>    omega1:=diff(phi1(t),t):
v2x:=-r1*omega1*sin(phi1(t)):
omega3:=omega1*R1/R3:
omega4:=v2x/(2*R4):
v4x:=omega4*R4:
v5x:=omega4*(R4-r4):

Кинетическая энергия системы

>    T1:=0.5*I1*omega1^2:
T2:=0.5*m2*v2x^2:
T3:=0.25*(m3*R3^2)*omega3^2:
T4:=0.5*m4*v4x^2+0.5*m4*i4^2*omega4^2:
T:=T1+T2+T3+T4;

T := .5*I1*diff(phi1(t),t)^2+.5*m2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2+.25*m3*diff(phi1(t),t)^2*R1^2+.12500*m4*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2+.12500*m4*i4^2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))...
T := .5*I1*diff(phi1(t),t)^2+.5*m2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2+.25*m3*diff(phi1(t),t)^2*R1^2+.12500*m4*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2+.12500*m4*i4^2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))...

Обобщенные силы

>    Q:=(MDz*omega3+Mnz*omega4+Fnx*v5x)/omega1;

Q := ((Mo-k*diff(phi1(t),t)*R1/R3)*diff(phi1(t),t)*R1/R3-1/4*mu*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2/R4^2-1/4*nu*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2/R4^2*(R4-r4)^2)/diff(phi1(t),t)
Q := ((Mo-k*diff(phi1(t),t)*R1/R3)*diff(phi1(t),t)*R1/R3-1/4*mu*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2/R4^2-1/4*nu*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2/R4^2*(R4-r4)^2)/diff(phi1(t),t)

Уравнение Лагранжа 2-го рода

>    eps:= (Q-0.5*B*omega1^2*sin(2*phi1(t)))/(A+B*sin(phi1(t))^2);

eps := (((Mo-k*diff(phi1(t),t)*R1/R3)*diff(phi1(t),t)*R1/R3-1/4*mu*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2/R4^2-1/4*nu*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2/R4^2*(R4-r4)^2)/diff(phi1(t),t)-.5*B*diff(phi1...
eps := (((Mo-k*diff(phi1(t),t)*R1/R3)*diff(phi1(t),t)*R1/R3-1/4*mu*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2/R4^2-1/4*nu*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2/R4^2*(R4-r4)^2)/diff(phi1(t),t)-.5*B*diff(phi1...
eps := (((Mo-k*diff(phi1(t),t)*R1/R3)*diff(phi1(t),t)*R1/R3-1/4*mu*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2/R4^2-1/4*nu*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2/R4^2*(R4-r4)^2)/diff(phi1(t),t)-.5*B*diff(phi1...

  Коэффициенты :

>    A:=evalf(coeff(algsubs(omega1^2*sin(phi1(t))^2=x,T*2),omega1^2));
B:= evalf(coeff(algsubs(omega1^2*sin(phi1(t))^2=x,T*2),x));
Qo   := evalf(subs(omega1=omega1o,phi1(t)=phi1o,Q));
Epso := evalf(subs(omega1=omega1o,phi1(t)=phi1o,eps));

A := 1.0*I1+.50*m3*R1^2

B := r1^2*(.25000*m4*i4^2+1.0*m2*R4^2+.25000*m4*R4^2)/R4^2

Qo := 1.0000*(Mo-.5*k*R1/R3)*R1/R3-.11606*mu*r1^2/R4^2-.11606*nu*r1^2/R4^2*(R4-1.*r4)^2

Epso := (1.0000*(Mo-.5*k*R1/R3)*R1/R3-.11606*mu*r1^2/R4^2-.11606*nu*r1^2/R4^2*(R4-1.*r4)^2-.64438e-1*r1^2*(.25000*m4*i4^2+1.0*m2*R4^2+.25000*m4*R4^2)/R4^2)/(1.0*I1+.50*m3*R1^2+.92845*r1^2*(.25000*m4*i4...
Epso := (1.0000*(Mo-.5*k*R1/R3)*R1/R3-.11606*mu*r1^2/R4^2-.11606*nu*r1^2/R4^2*(R4-1.*r4)^2-.64438e-1*r1^2*(.25000*m4*i4^2+1.0*m2*R4^2+.25000*m4*R4^2)/R4^2)/(1.0*I1+.50*m3*R1^2+.92845*r1^2*(.25000*m4*i4...
Epso := (1.0000*(Mo-.5*k*R1/R3)*R1/R3-.11606*mu*r1^2/R4^2-.11606*nu*r1^2/R4^2*(R4-1.*r4)^2-.64438e-1*r1^2*(.25000*m4*i4^2+1.0*m2*R4^2+.25000*m4*R4^2)/R4^2)/(1.0*I1+.50*m3*R1^2+.92845*r1^2*(.25000*m4*i4...

   Дифференциальное уравнение

>    Urav:=diff(omega1,t) - eps;

Urav := diff(phi1(t),`$`(t,2))-(((Mo-k*diff(phi1(t),t)*R1/R3)*diff(phi1(t),t)*R1/R3-1/4*mu*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2/R4^2-1/4*nu*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2/R4^2*(R4-r4)^2)/diff(p...
Urav := diff(phi1(t),`$`(t,2))-(((Mo-k*diff(phi1(t),t)*R1/R3)*diff(phi1(t),t)*R1/R3-1/4*mu*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2/R4^2-1/4*nu*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2/R4^2*(R4-r4)^2)/diff(p...
Urav := diff(phi1(t),`$`(t,2))-(((Mo-k*diff(phi1(t),t)*R1/R3)*diff(phi1(t),t)*R1/R3-1/4*mu*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2/R4^2-1/4*nu*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2/R4^2*(R4-r4)^2)/diff(p...
Urav := diff(phi1(t),`$`(t,2))-(((Mo-k*diff(phi1(t),t)*R1/R3)*diff(phi1(t),t)*R1/R3-1/4*mu*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2/R4^2-1/4*nu*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2/R4^2*(R4-r4)^2)/diff(p...

  Коэффициент а1 вариации omega1

>    a1:= factor(coeff(D(Urav),D(omega1)));

a1 := (R1^2*R4^2*k+.25000*mu*r1^2*sin(phi1(t))^2*R3^2+.25000*nu*r1^2*sin(phi1(t))^2*R3^2*R4^2-.50000*nu*r1^2*sin(phi1(t))^2*R3^2*R4*r4+.25000*nu*r1^2*sin(phi1(t))^2*R3^2*r4^2+.25000*r1^2*diff(phi1(t),t...
a1 := (R1^2*R4^2*k+.25000*mu*r1^2*sin(phi1(t))^2*R3^2+.25000*nu*r1^2*sin(phi1(t))^2*R3^2*R4^2-.50000*nu*r1^2*sin(phi1(t))^2*R3^2*R4*r4+.25000*nu*r1^2*sin(phi1(t))^2*R3^2*r4^2+.25000*r1^2*diff(phi1(t),t...
a1 := (R1^2*R4^2*k+.25000*mu*r1^2*sin(phi1(t))^2*R3^2+.25000*nu*r1^2*sin(phi1(t))^2*R3^2*R4^2-.50000*nu*r1^2*sin(phi1(t))^2*R3^2*R4*r4+.25000*nu*r1^2*sin(phi1(t))^2*R3^2*r4^2+.25000*r1^2*diff(phi1(t),t...
a1 := (R1^2*R4^2*k+.25000*mu*r1^2*sin(phi1(t))^2*R3^2+.25000*nu*r1^2*sin(phi1(t))^2*R3^2*R4^2-.50000*nu*r1^2*sin(phi1(t))^2*R3^2*R4*r4+.25000*nu*r1^2*sin(phi1(t))^2*R3^2*r4^2+.25000*r1^2*diff(phi1(t),t...
a1 := (R1^2*R4^2*k+.25000*mu*r1^2*sin(phi1(t))^2*R3^2+.25000*nu*r1^2*sin(phi1(t))^2*R3^2*R4^2-.50000*nu*r1^2*sin(phi1(t))^2*R3^2*R4*r4+.25000*nu*r1^2*sin(phi1(t))^2*R3^2*r4^2+.25000*r1^2*diff(phi1(t),t...

  Найдём коэффициенты и нарисуем графики зависимости угла и коэффициента a1(t) от времени

 

  Коэффициенты

>    A:=subs(Yslo1,A);
B:=subs(Yslo1,B);
Qo:=subs(Yslo1,Qo);
Epso:=subs(Yslo1,Epso);

A := 11.455

B := 1.9067

Qo := -13.821

Epso := -1.0544

Нарисуем графики зависимости угла и   коэффициента a1(t) от времени

>    Urav1:=subs(Yslo1,Urav);
F1:= dsolve({Urav1,NaYs},phi1(t),type = numeric);
a11:=subs(Yslo1,a1);
plots[odeplot](F1,[[t,phi1(t)],[t,a11]],t=0..2*Pi,title="График",legend=["Phi(t)","a1(t)"],thickness=2);

Urav1 := diff(phi1(t),`$`(t,2))-((1.3571*(10-20.357*diff(phi1(t),t))*diff(phi1(t),t)-29.251*diff(phi1(t),t)^2*sin(phi1(t))^2)/diff(phi1(t),t)-.95335*diff(phi1(t),t)^2*sin(2*phi1(t)))/(11.455+1.9067*sin...
Urav1 := diff(phi1(t),`$`(t,2))-((1.3571*(10-20.357*diff(phi1(t),t))*diff(phi1(t),t)-29.251*diff(phi1(t),t)^2*sin(phi1(t))^2)/diff(phi1(t),t)-.95335*diff(phi1(t),t)^2*sin(2*phi1(t)))/(11.455+1.9067*sin...

F1 := proc (x_rkf45) local res, solnproc, outpoint, ndsol, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 14; if _EnvInFsolve =...

a11 := 12.755*(.86640e-1+.9177e-1*sin(phi1(t))^2+.59795e-2*diff(phi1(t),t)*sin(2.*phi1(t)))/(.45820+.76270e-1*sin(phi1(t))^2)

[Maple Plot]

 Найдём условие нестабильности при значении t

>    ak:=subs(Yslo2,a1);

>    k:=solve(ak,k);

>    Urav2:=subs(Yslo2,Urav);

>    F2:=dsolve({Urav2,NaYs},phi1(t),type = numeric,output=listprocedure);

>    tk:=5.7;

>    phi1:=subs(F2,phi1(t)): phi1tk:=phi1(tk);

>    Dphi1:=subs(F2,diff(phi1(t),t)): Dphi1tk:=Dphi1(tk);

>    k:=evalf(subs(diff(phi1(t),t)=Dphi1tk,phi1(t)=phi1tk,k));

>    atk=evalf(subs(diff(phi1(t),t)=Dphi1tk,phi1(t)=phi1tk,ak));

>    plots[odeplot](F2,[t,ak],0..2*Pi,title="коэффициента a(t), при касании cо осью при t",color=blue,thickness=2,legend="Коэффициент a1");

>    plots[odeplot](F2,[t,ak],5.5..2*Pi,title="коэффициента a(t), при касании cо осью при t",color=blue,thickness=2,legend="Коэффициент a1");

ak := 12.755*(.5776e-2*k+.9177e-1*sin(phi1(t))^2+.59795e-2*diff(phi1(t),t)*sin(2.*phi1(t)))/(.45820+.76270e-1*sin(phi1(t))^2)

k := -15.888*sin(phi1(t))^2-1.0352*diff(phi1(t),t)*sin(2.*phi1(t))

Urav2 := diff(phi1(t),`$`(t,2))-((1.3571*(10-1.3571*(-15.888*sin(phi1(t))^2-1.0352*diff(phi1(t),t)*sin(2.*phi1(t)))*diff(phi1(t),t))*diff(phi1(t),t)-29.251*diff(phi1(t),t)^2*sin(phi1(t))^2)/diff(phi1(t...
Urav2 := diff(phi1(t),`$`(t,2))-((1.3571*(10-1.3571*(-15.888*sin(phi1(t))^2-1.0352*diff(phi1(t),t)*sin(2.*phi1(t)))*diff(phi1(t),t))*diff(phi1(t),t)-29.251*diff(phi1(t),t)^2*sin(phi1(t))^2)/diff(phi1(t...
Urav2 := diff(phi1(t),`$`(t,2))-((1.3571*(10-1.3571*(-15.888*sin(phi1(t))^2-1.0352*diff(phi1(t),t)*sin(2.*phi1(t)))*diff(phi1(t),t))*diff(phi1(t),t)-29.251*diff(phi1(t),t)^2*sin(phi1(t))^2)/diff(phi1(t...
Urav2 := diff(phi1(t),`$`(t,2))-((1.3571*(10-1.3571*(-15.888*sin(phi1(t))^2-1.0352*diff(phi1(t),t)*sin(2.*phi1(t)))*diff(phi1(t),t))*diff(phi1(t),t)-29.251*diff(phi1(t),t)^2*sin(phi1(t))^2)/diff(phi1(t...

F2 := [t = proc (t) option `Copyright (c) 1993 by the University of Waterloo. All rights reserved.`; t end proc, phi1(t) = proc (t) local res, solnproc, outpoint, ndsol; option `Copyright (c) 2000 by W...

tk := 5.7

phi1tk := 21.6138735378912444

Dphi1tk := 6.51303094351156098

k := 2.4623

atk = 0.

[Maple Plot]

[Maple Plot]

Выыод : При t=5.7[c] система не стабильна, если k=2.4623

Решение задачи с двумя обобщёнными координатами

Данные значения

>    restart; Digits:=5:

>    MDz:=Mo-k*omega3:
Mnz:=-mu*omega4:
Fnx:=-nu*v5x:

Начальные условия

>    Yslo2:=m2=16,m3=34,m4=26,nu=8000,mu=13,k=15,Mo=10,R1=0.38,I1=9,r1=0.27,R3=0.28,R4=0.2,
r4=0.12,i4=0.15,c=100:
phi1o:=1.3:omega1o:=0.5:
phi3o:=1.3:omega3o:=0.5:
NaYs:=phi1(0)=phi1o,D(phi1)(0)=omega1o,phi3(0)=phi3o,D(phi3)(0)=omega3o:

Имеем

>    omega1:=diff(phi1(t),t):
v2x:=-r1*omega1*sin(phi1(t)):
omega3:=diff(phi3(t),t):
omega4:=v2x/(2*R4):
v4x:=omega4*R4:
v5x:=omega4*(R4-r4):

Кинетическая энергия системы

>    T1:=0.5*I1*omega1^2:
T2:=0.5*m2*v2x^2:
T3:=0.25*(m3*R3^2)*omega3^2:
T4:=0.5*m4*v4x^2+0.5*m4*i4^2*omega4^2:
T:=T1+T2+T3+T4:

  Потенциальная энергия системы:  

>    П := 0.5*c*(R1*phi1(t)-R3*phi3(t))^2:

  Энергия системы

>    L:= T-П;

L := .5*I1*diff(phi1(t),t)^2+.5*m2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2+.25*m3*R3^2*diff(phi3(t),t)^2+.12500*m4*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2+.12500*m4*i4^2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))...
L := .5*I1*diff(phi1(t),t)^2+.5*m2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2+.25*m3*R3^2*diff(phi3(t),t)^2+.12500*m4*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2+.12500*m4*i4^2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))...
L := .5*I1*diff(phi1(t),t)^2+.5*m2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2+.25*m3*R3^2*diff(phi3(t),t)^2+.12500*m4*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2+.12500*m4*i4^2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))...

  Система уравнений

>    eq1:= diff(subs(w1=omega1,diff(subs(omega1=w1,L),w1)),t) - subs(f1=phi1(t),diff(subs(phi1(t)=f1,L),f1));
eq2:= diff(subs(w3=omega3,diff(subs(omega3=w3,L),w3)),t) - subs(f3=phi3(t),diff(subs(phi3(t)=f3,L),f3));

eq1 := 1.0*I1*diff(phi1(t),`$`(t,2))+1.0*m2*r1^2*diff(phi1(t),`$`(t,2))*sin(phi1(t))^2+1.0*m2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))*cos(phi1(t))+.25000*m4*r1^2*diff(phi1(t),`$`(t,2))*sin(phi1(t))^2+.2500...
eq1 := 1.0*I1*diff(phi1(t),`$`(t,2))+1.0*m2*r1^2*diff(phi1(t),`$`(t,2))*sin(phi1(t))^2+1.0*m2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))*cos(phi1(t))+.25000*m4*r1^2*diff(phi1(t),`$`(t,2))*sin(phi1(t))^2+.2500...
eq1 := 1.0*I1*diff(phi1(t),`$`(t,2))+1.0*m2*r1^2*diff(phi1(t),`$`(t,2))*sin(phi1(t))^2+1.0*m2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))*cos(phi1(t))+.25000*m4*r1^2*diff(phi1(t),`$`(t,2))*sin(phi1(t))^2+.2500...
eq1 := 1.0*I1*diff(phi1(t),`$`(t,2))+1.0*m2*r1^2*diff(phi1(t),`$`(t,2))*sin(phi1(t))^2+1.0*m2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))*cos(phi1(t))+.25000*m4*r1^2*diff(phi1(t),`$`(t,2))*sin(phi1(t))^2+.2500...

eq2 := .50*m3*R3^2*diff(phi3(t),`$`(t,2))-1.0*c*(R1*phi1(t)-R3*phi3(t))*R3

   Вариация:

>    eq11:= subs(sin(phi1(t))=sin@phi1(t),cos(phi1(t))=cos@phi1(t),eq1):
eq21:= subs(sin(phi1(t))=sin@phi1(t),cos(phi1(t))=cos@phi1(t),eq2):
eq11a:= D(eq11):
eq21a:= D(eq21):
eq11:=subs(sin@phi1(t)=sin(phi1(t)),(-sin)@phi1(t)=-sin(phi1(t)),cos@phi1(t)=cos(phi1(t)),eq11a):
eq21:=subs(sin@phi1(t)=sin(phi1(t)),(-sin)@phi1(t)=-sin(phi1(t)),cos@phi1(t)=cos(phi1(t)),eq21a):

   Найдём точки нестабильности

>    a1:= coeff(eq11,D(phi3(t)));
a2:= coeff(eq11,D(phi1(t)));
a3:= coeff(eq21,D(phi1(t)));
a4:= coeff(eq21,D(phi3(t)));

a1 := -1.0*c*R3*R1

a2 := 1.0*c*R1^2-1.0*m2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2-.25000*m4*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2-.25000*m4*i4^2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2/R4^2+1.0*m2*r1^2*diff(phi1(t),t)^2*cos...
a2 := 1.0*c*R1^2-1.0*m2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2-.25000*m4*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2-.25000*m4*i4^2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2/R4^2+1.0*m2*r1^2*diff(phi1(t),t)^2*cos...
a2 := 1.0*c*R1^2-1.0*m2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2-.25000*m4*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2-.25000*m4*i4^2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2/R4^2+1.0*m2*r1^2*diff(phi1(t),t)^2*cos...
a2 := 1.0*c*R1^2-1.0*m2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2-.25000*m4*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2-.25000*m4*i4^2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2/R4^2+1.0*m2*r1^2*diff(phi1(t),t)^2*cos...
a2 := 1.0*c*R1^2-1.0*m2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2-.25000*m4*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2-.25000*m4*i4^2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2/R4^2+1.0*m2*r1^2*diff(phi1(t),t)^2*cos...

a3 := -1.0*c*R3*R1

a4 := 1.0*c*R3^2

   Oпределитель

>    det:= a1*a3-a2*a4;

det := 1.00*c^2*R3^2*R1^2-1.0*(1.0*c*R1^2-1.0*m2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2-.25000*m4*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2-.25000*m4*i4^2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2/R4^2+1.0*m2*r...
det := 1.00*c^2*R3^2*R1^2-1.0*(1.0*c*R1^2-1.0*m2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2-.25000*m4*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2-.25000*m4*i4^2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2/R4^2+1.0*m2*r...
det := 1.00*c^2*R3^2*R1^2-1.0*(1.0*c*R1^2-1.0*m2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2-.25000*m4*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2-.25000*m4*i4^2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2/R4^2+1.0*m2*r...
det := 1.00*c^2*R3^2*R1^2-1.0*(1.0*c*R1^2-1.0*m2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2-.25000*m4*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2-.25000*m4*i4^2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2/R4^2+1.0*m2*r...
det := 1.00*c^2*R3^2*R1^2-1.0*(1.0*c*R1^2-1.0*m2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2-.25000*m4*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2-.25000*m4*i4^2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2/R4^2+1.0*m2*r...
det := 1.00*c^2*R3^2*R1^2-1.0*(1.0*c*R1^2-1.0*m2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2-.25000*m4*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2-.25000*m4*i4^2*r1^2*diff(phi1(t),t)^2*sin(phi1(t))^2/R4^2+1.0*m2*r...

   Численно

>    eq1:=subs(Yslo2,eq1);
eq2:=subs(Yslo2,eq2);
det:=subs(Yslo2,det);

eq1 := 9.0*diff(phi1(t),`$`(t,2))+1.9067*diff(phi1(t),`$`(t,2))*sin(phi1(t))^2+1.9067*diff(phi1(t),t)^2*sin(phi1(t))*cos(phi1(t))+14.440*phi1(t)-10.640*phi3(t)
eq1 := 9.0*diff(phi1(t),`$`(t,2))+1.9067*diff(phi1(t),`$`(t,2))*sin(phi1(t))^2+1.9067*diff(phi1(t),t)^2*sin(phi1(t))*cos(phi1(t))+14.440*phi1(t)-10.640*phi3(t)

eq2 := 1.3328*diff(phi3(t),`$`(t,2))-10.640*phi1(t)+7.8400*phi3(t)

det := 14.949*diff(phi1(t),t)^2*sin(phi1(t))^2-14.949*diff(phi1(t),t)^2*cos(phi1(t))^2-29.899*diff(phi1(t),`$`(t,2))*cos(phi1(t))*sin(phi1(t))
det := 14.949*diff(phi1(t),t)^2*sin(phi1(t))^2-14.949*diff(phi1(t),t)^2*cos(phi1(t))^2-29.899*diff(phi1(t),`$`(t,2))*cos(phi1(t))*sin(phi1(t))

    Pешение системы дифференциальных уравнений:

>    F:=dsolve({eq1,eq2,NaYs},{phi1(t),phi3(t)},type=numeric,output=listprocedure):

>    plots[odeplot](F,[[t,phi1(t)],[t,phi3(t)],[t,det]],0..2*Pi,title="График",legend=["Phi1(t)","phi3(t)","Det12(t)"],
thickness=2);

[Maple Plot]

Вывод : На промежутке времени [0,2Pi]  система имеет 2 точки нестабильности порядка [1/2]

Моделирование механизма в движении

>    restart;
pi:= evalf(Pi):
with(plots): with(plottools):

Warning, the name changecoords has been redefined
Warning, the name arrow has been redefined
 

 Процедура изображения окружности с центром в точке i

>    Circle:= proc(i,R) PLOT(circle([x[i],y[i]],R,color=blue)): end:

  Процедура изображения  точек на вращающемся колесе радиуса R c центром  в точке  i

>    Rot:= proc(i,R,f) local j:
 PLOT(seq(POINTS([x[i]+R*cos(f+pi/3*j),y[i]+R*sin(f+pi/3*j)]),j=0..5)):
end:

  Процедура изображения линии от точки i   к   j

>    Line:= proc(i,j,c) PLOT(CURVES([[x[i],y[i]],[x[j],y[j]]]),COLOR(HUE,c/10.)): end:

   Процедура изображения неподвижной опоры в точке    i

>    Rest:= proc(i,R) local x0,x1,y0,h,N:
x0:=x[i]-R*0.8: x1:= x[i]+R*0.8:
y0:=y[i]-R*0.8: h:= 3*R: N:= 2:
display(PLOT(circle([x[i],y[i]],R,color=blue)),PLOT(CURVES([[x0,y0],[x[i]-h,y[i]-h]],[[x1,y0],
[x[i]+h,y[i]-h]],[[x0-h,y[i]-h],[x1+h,y[i]-h]],seq([[x0-h*(-j/N+1),y[i]-h],[x0-h*(-j/N+1)-h/4,y[i]-h*1.3]],j=0..2*N+1)))):
end:

   Процедура изображения горизонтальной опорной поверхности (вверху h<0 , внизу h>0)

>    Surf:= proc(x,y,L,h) local N:
 N:= round(abs(L/h))-2:  display(PLOT(CURVES([[x,y],[x+L,y]],seq([[x+L*j/N/2,y],[x+L*j/N/2-h,y-h]],j=0..2*N)))):
end:

   Процедура изображения прямоугольника a x b c центром  в точке  i

>    Box:= proc(i,A,B) PLOT(POLYGONS([[x[i]-A,y[i]+B],[x[i]+A,y[i]+B],[x[i]+A,y[i]-B],[x[i]-A,y[i]-B]],COLOR(HUE,0.2))): end:

   Начальные значения

>    nam:=[O3,O1,A,E]; Nnam:= [1,6,7,10];
R1:= 38: r1:=27: R3:= 28: R4:= 20: r4:=12: l:=50:

nam := [O3, O1, A, E]

Nnam := [1, 6, 7, 10]

>    phio:= 1.3: omegao:= 0.5:

   Координаты основных точек

>    x[1]:= 0:                       y[1]:= 0:
x[2]:= 0:                       y[2]:= R3:
x[3]:= -R1-R3-l:                y[3]:= R1:
x[4]:= 0:                       y[4]:= -R3:
x[5]:= x[3]:                    y[5]:= -R1:
x[6]:= x[5]:                    y[6]:=0:
                                y[8]:= 0:
                                y[9]:=0:
                                y[10]:= y[9]+R4:
                                y[11]:= y[9]+R4+r4:
                                y[12]:= y[11]:

   Количество кадров

>    N:= 40:

  Создаем все кадры

>    for i from 0 to N do
 t:= 2*pi*i/N;
 phi:= phio + t:
 x[7]:= x[6] + r1*cos(phi): y[7]:= y[6] + r1*sin(phi):
 x[8]:= x[7]:
 x[9]:= x[8] - 5*l:
 x[10]:= x[9]/2-3*l:
 x[11]:= x[10]*(R4-r4)/R4-6*l:
 x[12]:= x[11] + 5*l:
 F:= arrow([x[11]-10,y[11]],evalm([30,0]),2,6,0.3):
 M:=arc([x[10],y[10]],R4-5,Pi/2..Pi,color=red):
 Mm:=arrow([x[10],y[10]+R4-5],evalm([7,-1]),1,5,0.5,color=red):
 M1:=arc([x[1],y[1]],R4-5,Pi/2..Pi,color=red):
 M1m:=arrow([x[1]-R4+5,y[1]],evalm([-1,-7]),1,5,0.5,color=red):
 P[i]:= display(F,Box(8,2,r1),Circle(1,R3), Rot(1,R3,phi*R1/R3),
                              Circle(6,R1), Rot(6,R1,phi),
                              Circle(10,R4),Rot(10,R4,x[10]/R4),
                              Circle(10,r4),Rot(10,r4,x[10]/R4),
                              Circle(10,1),
                              Mm,M,M1m,M1,  
                  seq(Line(2*k,2*k+1,1),k=1..2),
                      Line(6,7,0),
                      Line(8,9,0),
                      Line(11,12,0),
                  seq(TEXT([x[Nnam[j]]+8,y[Nnam[j]]],nam[j]),j=1..4));
od:

   Изображение механизма в движении

>    PP:=display(seq(P[i],i=0..N),insequence=true,thickness=2,axes=none):

>    display(                       PP,Rest(1,2),Rest(6,2),
                               Surf(-8.4*l,y[10]+R4,4*l,-6),
                               Surf(-3.5*l,y[8]+2,10,-2),
                               Surf(-3.5*l,y[8]-2,10,2),
                               Surf(-5*l,y[11]-2,10,2),
                               Surf(-5*l,y[11]+2,10,-2),
                               Surf(-7.8*l,y[11]-2,10,2),
                               Surf(-7.8*l,y[11]+2,10,-2));

[Maple Plot]

Наверх