logo
[ЛР]ЭиИСУЛА[Илюхин А

Программа для расчета динамических характеристик исполнительного двигателя системы рулевого привода.

Описание программы.

Программа предназначена для математического моделирования пневматического исполнительного двигателя с клапанным газораспределителем (время движения затвора – пренебрежимо мало) при отработке ступенчатого входного сигнала. Математическое моделирование осуществляется путем численного решения системы дифференциальных уравнений, описывающих работу привода (п.1.2.).

Структура программы.

Программа состоит из следующих подпрограмм.

Главная (вызывающая) – выполняет ввод исходных данных, обращение к процедуре EU, вывод результатов.

Процедура EUРеализация метода Ейлера, обращение к процедуре RP.

Процедура RPвычисление функций правых частей системы уравнений, обращение к функции G.

Функция G – вычисление секундного массового расхода газа.

Программные обозначения.

Параметры, которые задаются в процедуре RP.

Ds1, Ds2 – диаметры впускных отверстий в первую и вторую полости соответственно (м);

Dpotr1, Dpotr2 - диаметры выпускных отверстий из первой и второй полости соответственно (М).

Параметры, которые задаются в вызывающей программе.

Jm – порядок истемы уравнений;

t0 начальное значение времени (с);

tk :=0.1- время интегрирования (с);

h:= 0.0000001 – шаг интегрирования (с);

xp := 0.0 – «время» печати (с);

hp := 0.0010 – «шаг» печати (с);

mu:=0.9 – коэффициент расхода;

Pv:=20.0*100000 – давление газа в источнике питания (первый операнд, бар);

Pa:=1.0*100000 – давление газа в отсеке (первый операнд, бар);

Tv:= 293.16 - температура газа в источнике питания (град. К);

Ta:= 293.16 - температура газа в отсеке (град. К);

m:=0.05 – масса подвижных частей (кг);

Xmax:=0.007 – максимальный ход поршня (м);

per:=0.02 – длительность периода (с);

gam:=0.5 - скважность;

nu:=0.0 – коэффициент нагрузки (н/м);

xs:=0.0 – начальные условия по координате;

v:=0.0 – начальные условия по скорости;

dp:=0.007 – диаметр поршня (м);

Sp:=(pi*Dp*Dp)/4;

k:=1.4 – коэффициент абиабаты;

Bk:=0.528 – критическое отношение давлений;

k0:=exp(0.5*ln(k)+((k+1)*(ln(2/(k+1)))/(2*(k-1))));

R:=287.0 – газовая постоянная (дж/кг*К);

Текст программы.

Program PnevmoPrivod;

Label 1;

type

mas = array[1..6] of real;

var

t0, h, t, tk :real;

yin, yout, fout :mas;

alfa,x,xp, hp :real;

F :text;

i,jm :integer;

per,gam,tf,nnn,Fdv,Fro,Fsum,mu,sv,Pv,Pnas,Dpotr1,Dpotr2,Pa,Tv:real;

rn,v,tn,pn,T1,Ta,Dp,Dp1,Dp2,Ds,Ds1,Ds2,Sp,Ss,Dpotr,Spotr,nu,htr,k,Bk:real;

xs,Gv,Gn,Xmax,k0,ss1,ss2,spotr1,spotr2,W,T2,R,M,W1,Fp,pp2,pn1:real;

function g(p1,p2,r,t,mu,s:real):real;

var p,y:real;

begin

p:=p2/p1;

if p>=1.0 then Y:=0.0

else if p>bk then

Y:=(1/k0)*(sqrt((2*k)/(k-1)))*(sqrt((exp(ln(p)*2/k))-(exp(ln(p)*(k+1)/k))))

else Y:=1;

G:=(mu*k0*s*p1*Y)/(Sqrt(r*t));

end;

procedure RP(var t0 :real;var yin,fout :mas);

var ppv,pp1,x,y,z:real;

begin

tf:=t0-(per*nnn);

if tf>=(per*gam) then begin

Ds1:=0.002;

Ds2:=0.0;

Dpotr1:=0.0;

Dpotr2:=0.004;

Ss1:=(pi*Ds1*Ds1)/4;

Ss2:=(pi*Ds2*Ds2)/4;

Spotr1:=(pi*Dpotr1*Dpotr1)/4;

Spotr2:=(pi*Dpotr2*Dpotr2)/4;

end else begin

Ds1:=0.0;

Ds2:=0.002;

Dpotr1:=0.004;

Dpotr2:=0.0;

Ss1:=(pi*Ds1*Ds1)/4;

Ss2:=(pi*Ds2*Ds2)/4;

Spotr1:=(pi*Dpotr1*Dpotr1)/4;

Spotr2:=(pi*Dpotr2*Dpotr2)/4;

end;

if tf>=per then nnn:=nnn+1;

T1:=yin[1]/(yin[2]*r);

ppv:=(k*r*Tv)/(k-1);

pp1:=(k*r*T1)/(k-1);

Gv:=g(pv,yin[1],r,Tv,mu,ss1);

Gn:=g(yin[1],pa,r,T1,mu,spotr1);

fout[1]:=(k-1.0)*(ppv*Gv-pp1*Gn-(k*yin[1]*Sp*yin[5])/(k-1.0))/(W+Sp*yin[6]);

fout[2]:=(Gv-Gn-yin[2]*Sp*yin[5])/(W+Sp*yin[6]);

T2:=yin[3]/(yin[4]*r);

ppv:=(k*r*Tv)/(k-1);

pp1:=(k*r*T2)/(k-1);

Gv:=g(pv,yin[3],r,Tv,mu,ss2);

Gn:=g(yin[3],pa,r,T2,mu,spotr2);

fout[3]:=(k-1.0)*(ppv*Gv-pp1*Gn+(k*yin[3]*Sp*yin[5])/(k-1.0))/(W-Sp*yin[6]);

fout[4]:=(Gv-Gn+yin[4]*Sp*yin[5])/(W-Sp*yin[6]);

Fdv:=Sp*(yin[1]-yin[3])-h*yin[5]-nu*yin[6];

{1} if (yin[6]>-xmax) and (yin[6]<xmax) then

Fro:=0.0

else

{2} if (Fdv>0.0) and (yin[5]>0.0) and (yin[6]>=xmax)

then begin

Fro:=Fdv;

yin[6]:=xmax;

yin[5]:=0.0

end

else

{3} if (yin[6]=xmax) and (Fdv>0.0) and (yin[5]=0.0)

then begin

Fro:=Fdv;

yin[6]:=xmax;

yin[5]:=0.0

end

else

{4} if (yin[6]>=xmax) and (Fdv<0.0) and (yin[5]>0.0)

then begin

Fro:=0.0;

yin[6]:=xmax;

yin[5]:=0.0

end

else

{5} if (yin[6]=xmax) and (Fdv<=0.0) and (yin[5]<=0.0)

then

Fro:=0.0

else

{6} if (yin[6]<=-xmax) and (Fdv<0.0) and (yin[5]<0.0)

then begin

Fro:=Fdv;

yin[6]:=-xmax;

yin[5]:=0.0

end

else

{7} if (yin[6]=-xmax) and (Fdv<0.0) and (yin[5]=0.0)

then begin

Fro:=Fdv;

yin[6]:=-xmax;

yin[5]:=0.0

end

else

{8} if (yin[6]<=-xmax) and (Fdv>0.0) and (yin[5]<0.0)

then begin

Fro:=0.0;

yin[6]:=-xmax;

yin[5]:=0.0

end

else

{9} if (yin[6]=-xmax) and (Fdv>=0.0) and (yin[5]>=0.0)

then

Fro:=0.0;

fout[5]:=1/m*(Fdv-Fro);

fout[6]:=yin[5];

{ end; }

end;

procedure eu(var t0, t, h, xp :real; varjm :integer;var yin, yout :mas);

var i :integer;

{ fout :mas;}

begin

while t <= xp do

begin

rp(t0, yin, fout);

for i := 1 to jm do

yout[i] := yin[i] + h * fout[i];

t:= t0 + h;

for i := 1 to jm do

yin[i] := yout[i];

t0 :=t;

end;

end;

BEGIN

jm:=6;

t0 := 0.0;

tk :=0.1;

h:= 0.0000001;

xp := 0.0;

hp := 0.001;

mu:=0.9;

Pv:=2.0*100000*10;

{Pnas:=1.0*100000*20;}

Pa:=1.0*100000*1;

Tv:= 293.16;

Ta:= 293.16;

T1:= 293.16;

t2:=t1;

Ds1:=0.002;

Ds2:=0.0;

Dpotr1:=0.0;

Dpotr2:=0.002;

Ss1:=(pi*Ds1*Ds1)/4;

Ss2:=(pi*Ds2*Ds2)/4;

m:=0.05;

Spotr1:=(pi*Dpotr1*Dpotr1)/4;

Spotr2:=(pi*Dpotr2*Dpotr2)/4;

Xmax:=0.007;

per:=0.02;

gam:=0.5;

nu:=0;

nnn:=0;

xs:=0;

v:=0;

dp:=0.007;

Sp:=(pi*Dp*Dp)/4;

k:=1.4; {}

Bk:=0.528;

k0:=exp(0.5*ln(k)+((k+1)*(ln(2/(k+1)))/(2*(k-1)))); {}

writeln(k0);

R:=287; {}

W:=Sp*Xmax*1.2; {0.2 litra}

yin[1]:=Pa;

yin[2]:=yin[1]/(R*t1);

yin[3]:=Pa;

yin[4]:=yin[3]/(R*t2);

yin[5]:=xs;

yin[6]:=v;

Writeln(yin[1],' ',yin[2]);

Writeln(yin[3],' ',yin[4]);

writeln(yin[5],' ',yin[6]);

writeln;

{ Writeln(f,' P 1 ',' P 2 ',' X ',' V ',' t '); }

Writeln(' P 1 ',' P 2 ',' X ',' V ',' t ');

while t<=tk do

begin

eu(t0, t, h, xp, jm, yin, yout);

pn:=yout[1]/100000.0; {Pa v atmosfer}

pn1:=yout[3]/100000.0;

tn:=yout[1]/(yout[2]*r);

for i := 1 to jm do

yin[i] := yout[i];

{If yin[5]>Xmax then goto 1;}

t0:= t;

xp := xp + hp;

writeln(pn:5:2,' ',pn1:5:2,' ',yin[6]:5,' ',yin[5]:5,' ',t:3{,' ',Fro,' ',Fdv},tf:3);

{ writeln(f,pn:5:2,' ',pn1:5:2,' ',yin[6]:5,' ',yin[5]:5,' ',t:3);

writeln(f,Ds1,' ',Ds2,' ',Dpotr1,' ',Dpotr2,' ',t:3); }

{readkey;}

end;

{Close(f);}

END.