3.1.3. Опис програмного забезпечення
Моделювання поводження двоступеневого гірокомпаса (ДГК) здійснюється за допомогою програми (для середовища Matlab), яка здійснює чисельне інтегрування рівняння (3.13), в якому до моментів сил навколо вимірювальної осі ДГК включені моменти сил в’язкого і сухого тертя.
Пакет програм і процедур мовою системи MatLAB, складається з наступних програмних одиниць:
– керуючої програми DGK_upr. m, що здійснює присвоювання значень усім параметрам рівняння, організує цикл чисельного інтегрування за допомогою процедури rko4 методу Рунге-Кутта четвертого порядку і виведення результатів інтегрування у графічне;
– процедури DGK, що обчислює значення правих частин диференційних рівнянь двоступеневого гірокомпаса, приведених до форми Коши. На її спирається процедура rko4 методу чисельного інтегрування.;
– процедури RKO4, що виконує чисельне інтегрування ЗДР, праві частини яких задані процедурою "ypfun", методом Рунге-Кутта 4-го порядку.
С
хема
взаємодії програм наведена нижче.
Тексти зазначених файлів наведено нижче.
Програма DGK_upr. m
% DGK_upr1
% Керуюча програма інтегрування диф. рівнянь двоступеневого гірокомпаса
% Воропаєва Є. О.
% Почато 09-01-2009 Останні зміни 19-01-2009
clear all, clc
global H J1 J2 omZ_fig fiN fiE Fiv om e al N F
sname='Воропаєва Є. О.';
% Задание значений исходных данных
omZ=7.29e-5; fig= 60*pi/180; omZ_fig=omZ*cos(fig);
% 1. Кінетичний момент
H = 4000;
% 2. Моменти інерції
J1=10; J2=3;
% 3. Коефіцієнт тертя
F=1.5;
% Момент сил сухого тертя
N=0.0;
% 4. Початкові умови
bt0=0.*pi/180; bt0T=0;
% 3. Параметри руху основи
fiN=0.5*pi/180/60; fiE=0.5*pi/180/60; Fiv=0*pi/180;
om=1; e=0*pi/180; al=0*pi/180;
% 4. Параметри інтегрування
% Кількість кроків інтегрування у найменшому періоді коливань
ni=100;
% Кількість найбільших періодів коливань
ntk=1.9;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Підготовка даних для інтегрування
OM0=sqrt(H*omZ_fig/J1); hi=2*pi/om/ni; TK=2*pi/OM0*ntk;
t=0; k=1;
y=[bt0, bt0T]; Beta(1)=y(1); BetaT(1)=y(2); tt(1)=t;
% Цикл інтегрування
while t<=TK
k=k+1; [tout,yout] = rko4('DGK11',hi,t,y);
y=yout; t=tout; Beta(k)=y(1); BetaT(k)=y(2); tt(k)=tout;
end
% Виведення графіка
subplot(3,1,[2:3])
plot(tt,Beta*180/pi,'k.'),grid
set(gca, 'fontname','MS Sans Serif')
title('Відхилення від площини меридіану','Fontsize',12,'FontWeight','bold')
xlabel('Час (с)','Fontsize',12), ylabel('\beta (gradus)','fontname','Helvetica','Fontsize',12)
% Текстове оформлення графічного вікна
subplot(3,1,1), axis off
h=text(0.2,1.17,' Рух Двоступеневого ГіроКомпаса',...
'Fontsize',14,'FontWeight','bold','fontname','MS Sans Serif');
h=text(0.7,1.17,sprintf('(\\phi_g = %g^o) ',fig*180/pi),'Fontsize',12);
h=text(0.0,0.9,'J_1(\beta" + u_y'') + (H - J_2 u_x_1)u_z_1 = M_y_1','Fontsize',14);
h=text(0.0,0.7,'u_x_1= (\omega_3cos\phi_g + u_x)cos\beta -u_zsin\beta;','Fontsize',12);
h=text(0.0,0.5,'u_z_1=(\omega_3cos\phi_g + u_x)sin\beta + u_zcos\beta','Fontsize',12);
h=text(0.0,0.3,'u_x= \phi_N_m\omega*sin(\omegat+e)','Fontsize',12);
h=text(0.0,0.1,'u_y= \phi_v_m\omega*cos(\omegat+a)','Fontsize',12);
h=text(0.0,-0.1,'u_z= \phi_E_m\omega*sin(\omegat)','Fontsize',12);
h=text(0.5,0.9,sprintf('H = %g',H),'Fontsize',12);
h=text(0.5,0.7,sprintf('J_1 = %g',J1),'Fontsize',12);
h=text(0.5,0.5,sprintf('J_2 = %g',J2),'Fontsize',12);
h=text(0.7,0.9,sprintf('\\phi_N_m = %g ''',fiN*180/pi*60),'Fontsize',12);
h=text(0.7,0.7,sprintf('\\phi_v_m = %g ''',Fiv*180/pi*60),'Fontsize',12);
h=text(0.7,0.5,sprintf('\\phi_E_m = %g ''',fiE*180/pi*60),'Fontsize',12);
h=text(0.9,0.9,sprintf('\\omega = %g',om),'Fontsize',12);
h=text(0.9,0.7,sprintf('e = %g^o',e*180/pi),'Fontsize',12);
h=text(0.9,0.5,sprintf('a = %g^o',al*180/pi),'Fontsize',12);
h=text(0.5,0.3,sprintf('Тертя f = %g M_t_r = %g',[F,N]),'Fontsize',12,'fontname','MS Sans
Serif');
h=text(0.5,0.2,sprintf('Початкові умови \\beta_0 = %g \beta_0''=
%g',[bt0*180/pi,bt0T]),'Fontsize',12,'fontname','MS Sans Serif');
h=text(0.25,0.2,sprintf('ПеріодT_0=%g с',round(2*pi*sqrt(J1/H/omZ_fig))),
'Fontsize',12,'fontname','MS Sans Serif');
h=text(0.4,0.1,'==================================================');
h1=text(0.4,0.03,'Програма DGK-upr.m; Автор - Воропаєва Є. О., каф. ПСОН',
'fontname','MS Sans Serif');
h1=text(0.4,-0.04,['Виконав ',sname],'fontname','MS Sans Serif');
tm=fix(clock); Tv=tm(4:5); h1=text(0.7,-0.04,[sprintf(' %g:',Tv),' ',date]);
h=text(0.4,-0.1,'==================================================');
Програма DGK
function z=DGK1(t,y)
% Процедура обчислювання правих частин диференційних рівнянь
% руху Двоступеневого ГіроКомпаса
% z(1)=y(2);
% z(2)= -ucycT+(My1-(H-J2*ucx1)*ucz1)/J1;
% Вхідні параметри:
% t - поточний час
% y - масив з двох змінних стану
% y(1) = beta - кут відхилення від меридіану
% y(2) = betaT - кутова швидкість
% Воропаєва Є. О. 09-01-2009
global H J1 J2 omZ_fig fiN fiE Fiv om e al N F
sB=sin(y(1)); cB=cos(y(1));
Uksi=fiN*om*sin(om*t+e); Udz=fiE*om*sin(om*t); Ucksi=omZ_fig+Uksi;
ucx1=Ucksi*cB-Udz*sB; ucz1=Ucksi*sB+Udz*cB;
fivT=Fiv*om*cos(om*t+al); ucycT=-Fiv*om^2*sin(om*t+al);
alT=y(2)-fivT; My1=- F*alT-N*sign(alT);
z(1)=y(2);
z(2)= -ucycT+(My1-(H-J2*ucx1)*ucz1)/J1;
Програма RKO4
function [tout,yout] = rko4(ypfun,h,t,y)
% процедура RKO4 чисельного інтегрування ЗДР, праві частини яких задани
% процедурою ypfun, методом Рунге-Кутта 4-го порядку.
% Вхідні змінні:
% ypfun - символьний рядок, що містить імя п/п правих частин ЗДР
% Виклик: z = fun(t,y), де ypfun = 'fun'. t - поточний час;
% y - вектор поточних значень змінних стану,
% z - розраховані значень похідних; z(i) = dy(i)/dt.
% h - крок інтегрування; t - поточний час;
% y - значеня вектора змінних стану на попередньому кроці
% Вихідні змінні:
% tout - наступний момент часу;
% yout - розраховане значення вектора змінних стану у цей момент
% Y. F. Lazarev, 236-67-31.
% Розрахунок проміжних значень похідних
s1 = feval(ypfun, t, y);
s2 = feval(ypfun, t+h/2, y+h/2*s1);
s3 = feval(ypfun, t+h/2, y+h/2*s2);
s4 = feval(ypfun, t+h, y+h*s3);
% Розрахунок значення вектора змінних стану у наступний момент часу
tout = t + h;
yout = y + h*(s1 + 2*s2 + 2*s3 + s4)/6;