% уравн. динамики двуcтороннего привода с нач.перепадом давлений
global k xi01 xi02 omega21 sigmaa N Pi21 chi razm sigmam perem iperem
global D0 D1 D2 D3 nuc nup
k=1.4;
nuc=0;nup=0;
chim=[0.1 0.2 0.3 0.4];
xi01m=[0.3 0.3 0.3 0.3];
xi02m=[0.3 0.3 0.3 0.3];
omega21m=[1 1 1 1];
Nm=[0.1 0.1 0.1 0.1];
sigmaam=[0.167 0.167 0.167 0.167];
sigmamm=[1 1 1 1];
Pi21m=[1 1 1 1];
nazv1='дв.прив.с нач.переп.давл.при ';
nazv2='перех.проц.в двустор.пн.приводе с нач.перепадом давл. при ';
%измен.параметр
perem=chim;
iperem='chi';
jk=4; %колич.расч.вар.
%задание знач.перем.
for j=1:jk;
xi01=xi01m(j); xi02=xi02m(j);
omega21=omega21m(j);
sigmaa=sigmaam(j); sigmam=sigmamm(j);
N=Nm(j); Pi21=Pi21m(j); chi=chim(j);
%тип привода
%с нач.перепадом давл.
D1=1;
%вых.в атмосферу
D2=sigmaa;
%давл.в выхл.п.конеч.
D3=sigmaa; D0=1;
%врем. до начала движ. поршня
% начальные усл.(давление в рабочей и выхлопной полоcти)
y30=[sigmaa sigmam];
% интервал интегрирования
int3=0:0.0005:15;
[T3,Y3]=ode45(@ndy3,int3,y30);
kol=size(T3,1);
i=1;
while (Y3(i,1)-Pi21*Y3(i,2)-chi)<=0
if Y3(i,1)>=sigmam Y3(i,1)=sigmam;end
i=i+1;end
k3=i-1;
T3(k3+1:kol)=[];
Y3(k3+1:kol,:)=[];
Y_3=zeros(k3,4);
Y_3(:,4)=Y3(:,2);
Y_3(:,3)=Y3(:,1);
t3=T3(k3);
%врем.движ.поршн
%начальные усл. (перемещение, скорость перемещения,...
%давление в рабочей полости, давление в выхлопной полости)
y0=[0 0 Y_3(k3,3) Y_3(k3,4)];
%интервал интегрирования
int=T3(k3):0.0005:T3(k3)+6;
[T,Y]=ode45(@ndy,int,y0);
kol1=size(T,1)
i=1;
while Y(i,1)<=1
i=i+1;end
k11=i-1;
T(k11+1:kol1)=[];
Y(k11+1:kol1,:)=[];
t=T(k11)-t3;
%врем. после оcтановки поршня
% начальные усл.(давление в рабочей и выхлопной полоcти)
ys0=[Y(k11,3) Y(k11,4)];
% интервал интегрирования
ints=T(k11):0.0005:T(k11)+2;
[Ts,Ys]=ode45(@NDY111,ints,ys0);
kols=size(Ts,1);
i=1;
while or(Ys(i,1)<=0.95*sigmam,Ys(i,2)>=1.05*D3)
if Ys(i,1)>=sigmam Ys(i,1)=sigmam;end
if Ys(i,2)<=D3 Ys(i,2)=D3;end
i=i+1;end
ks=i-1
Ts(ks+1:kols)=[];
Ys(ks+1:kols,:)=[];
Y111=zeros(ks,4);
Y111(:,1)=1;
Y111(:,3)=Ys(:,1);
Y111(:,4)=Ys(:,2);
ts=Ts(ks)-t3-t;
% объединение массив.
tdv=cat(1,T3,T,Ts);
ydv=cat(1,Y_3,Y,Y111);
kdv=size(tdv,1);
t3m(j)=t3;ttm(j)=t;tsm(j)=ts
%построен.граф.
figure(1);
subplot(2,jk/2,j)
plot(tdv(:),ydv(:,1),'-',tdv(:),ydv(:,2),':',...
tdv(:),ydv(:,3),'-.',tdv(:),ydv(:,4),'--')
title([nazv1 iperem '=' num2str(perem(j))]),xlabel('t')
axis([0 tdv(kdv) -0.05 1.05])
grid
switch j
case 1,kdv1=kdv;,tdv1=tdv;,ydv1=ydv;
case 2,kdv2=kdv;,tdv2=tdv;,ydv2=ydv;
case 3,kdv3=kdv;,tdv3=tdv;,ydv3=ydv;
case 4,kdv4=kdv;,tdv4=tdv;,ydv4=ydv;
end
end
figure(1);legend('x','v','p1','p2',0)
for i=1:4
figure(i+1);
subplot(3,1,[1 2]);
plot(tdv1,ydv1(:,i),'-',tdv2,ydv2(:,i),':',...
tdv3,ydv3(:,i),'-.',tdv4,ydv4(:,i),'--')
title([nazv2 'изменении ' iperem])
xlabel('t')
legend([iperem '=' num2str(perem(1))],[iperem '=' num2str(perem(2))],...
[iperem '=' num2str(perem(3))],[iperem '=' num2str(perem(4))],0)
grid
switch i
case 1,ylabel('перемещение,x')
case 2,ylabel('скорость,v')
case 3,ylabel('давл.раб.пол.,p1')
case 4,ylabel('давл.выхл.пол.,p2')
end
subplot(3,1,3);
axis('off')
h1=text(0,0.9,'вариант:','FontSize',9);
h1=text(0.45,0.9,' 1 2 3 4','FontSize',9);
h1=text(0,0.8,'нагрузка chi:','FontSize',9);
h1=text(0.45,0.8,sprintf(' %7.2f',chim),'FontSize',9);
h1=text(0,0.7,'нач.коорд.раб.п. xi01:','FontSize',9);
h1=text(0.45,0.7,sprintf(' %7.2f',xi01m),'FontSize',9);
h1=text(0,0.6,'нач.коорд.выхл.п. xi02:','FontSize',9);
h1=text(0.45,0.6,sprintf(' %7.2f',xi02m),'FontSize',9);
h1=text(0,0.5,'коэф.пропускн.спос-ти omega21:','FontSize',9);
h1=text(0.45,0.5,sprintf(' %8.1f',omega21m),'FontSize',9);
h1=text(0,0.4,'обобщ.констр.параметр N:','FontSize',9);
h1=text(0.45,0.4,sprintf(' %8.1f',Nm),'FontSize',9);
h1=text(0,0.3,'отн.атм.давл sigmaa:','FontSize',9);
h1=text(0.45,0.3,sprintf('%7.3f',sigmaam),'FontSize',9);
h1=text(0,0.2,'отн.магистр.давл. sigmam:','FontSize',9);
h1=text(0.45,0.2,sprintf(' %7.2f',sigmamm),'FontSize',9);
h1=text(0,0.1,'отн.F2/F1 Pi21:','FontSize',9);
h1=text(0.45,0.1,sprintf(' %7.2f',Pi21m),'FontSize',9);
h1=text(0,0.0,'врем.подгот.этапа:','FontSize',9);
h1=text(0.45,0.0,sprintf(' %7.2f',t3m),'FontSize',9);
h1=text(0,-0.1,'врем.движ.:','FontSize',9);
h1=text(0.45,-0.1,sprintf(' %7.2f',ttm),'FontSize',9);
h1=text(0,-0.2,'врем.заключ.этапа:','FontSize',9);
h1=text(0.45,-0.2,sprintf(' %7.2f',tsm),'FontSize',9);
h1=text(0,-0.3,'врем.перех.проц. tdv:','FontSize',9);
h1=text(0.45,-0.3,sprintf(' %7.2f',[tdv1(kdv1) tdv2(kdv2) tdv3(kdv3) tdv4(kdv4)]),'FontSize',9);
end
Подпрограмма NDY3 расчета подготовительного этапа прямого хода в безразмерных величинах.
% описание системы диф. уравнений
function dy=ndy3(t3,y3)
global k xi01 xi02 omega21 sigmaa N Pi21 chi D0 D1 D2 D3 nuc nup
dy=zeros(2,1);
dy=[k*ras(y3(1))/xi01;...
-k*D0*(omega21*y3(2)^((3*k-1)/(2*k))*ras(D2/y3(2))/(Pi21*D1))/(xi02+1)];
if (y3(1)-Pi21*y3(2)-chi)>0
return
end
Подпрограмма NDY расчета этапа движения поршня в безразмерных величинах.
% описание системы диф. уравнений
function dy=ndy(t,y)
global k xi01 xi02 omega21 sigmaa N Pi21 chi D0 D1 D2 D3 nuc nup
dy=zeros(4,1);
if y(1)>1.02
return
else
dy=[y(2);...
(y(3)-Pi21*y(4)-chi-nup*y(1)-nuc*y(2))/(N^2);...
k*(ras(y(3))-y(3)*y(2))/(xi01+y(1)); ...
-k*D0*(omega21*y(4)^((3*k-1)/(2*k))*ras(D2/y(4))/(Pi21*D1)-y(4)*y(2))/(xi02+1-y(1))];
end
Подпрограмма NDY111 расчета заключительного этапа прямого хода в безразмерных величинах.
% описание системы диф. уравнений
function dy=ndy111(ts,ys)
global k xi01 xi02 omega21 sigmaa sigmam N Pi21 chi razm D0 D1 D2 D3 nuc nup
dy=zeros(2,1);
dy=[k*ras(ys(1))/(xi01+1); ...
-k*D0*omega21*ys(2)^((3*k-1)/(2*k))*ras(D2/ys(2))/(Pi21*xi02*D1)];
Подпрограмма ras расчета расходной функции
function Rf=ras(x)
global k
if and(x>0,x<=0.528) Rf=0.2588;
else
if and(x>0.528,x<=1)
Rf=sqrt(x^(2/k)-x^((k+1)/k));
else Rf=0;
end
end
