clc; clear; clear all;
Ixx=1e3; Iyy=3e3; Ixy=2e3; // осевые и центробежные моменты инерции шасси;
ri=1e-5; // радиус инерции колеса относительно оси вращения
T=0.1; r=0.2; // вынос и радиус колеса
a=1e1; b=1e3; // поперечная и продольная жёсткость колеса;
k=9e-1; // коэффициент бокового увода колеса
h=1e2; // коэффициент вязкого трения демпфера шимми
l=0.5; // длина стойки
V=200; // скорость движения самолёта
Cpsi=1e-5;
Ctheta=1e-5;
A=[1 0 0 0 0;
0 1 0 0 0;
0 0 1 0 0;
0 0 0 Ixx Ixy;
0 0 0 Ixy,Iyy];
global dlambda;
function dq=difur(t,q0)
global dlambda;
dlambda=-(l*q0(4)+t*q0(5)+V*a/k*q0(3));
dq=inv(A)*...
[
q0(4);
q0(5);
dlambda;
-(Cpsi*q0(1)+V*ri/r*q0(5)-a*l*q0(3));
(q0(4)*ri*V)/r-(q0(5)*k*T^2)/V-(q0(4)*k*l*T)/V-(dlambda*k*T)/V-...
(b*q0(5)*T)/V-...
(b*q0(4)*l)/V-...
(b*dlambda)/V-...
Ctheta*q0(2)+...
q0(5)*h;
]
endfunction
//Начальные значения dpsi, dtheta, psi, theta:
dpsi0=0;
dtheta0=0;
//dlambda0=0;
dtheta10=0;
psi0=0;
theta0=1e-1;
lambda0=0;
theta1=0;
U0=[psi0 theta0 lambda0 dpsi0 dtheta0]';
//время начальное, итоговое, шаг
t0=0;
tk=60;
dt=1e-2
Result=zeros((tk-t0)/dt+1,8);// сохраняем результаты
for i=1🙁tk-t0)/dt do
U0=ode(U0,(i-1)*dt,i*dt,difur);
Result(i,: )=[U0(1:3)', -(h/Ctheta*U0(5)-U0(2)), U0(4:5)', dlambda, 0];
end
i=0;
i=i+1; subplot(2,4,i); plot(Result🙂,i)); xtitle("psi")
i=i+1; subplot(2,4,i); plot(Result🙂,i)); xtitle("theta")
i=i+1; subplot(2,4,i); plot(Result🙂,i)); xtitle("lambda")
i=i+1; subplot(2,4,i); plot(Result🙂,i)); xtitle("theta1")
i=i+1; subplot(2,4,i); plot(Result🙂,i)); xtitle("dpsi")
i=i+1; subplot(2,4,i); plot(Result🙂,i)); xtitle("dtheta")
i=i+1; subplot(2,4,i); plot(Result🙂,i)); xtitle("dlambda")
