高压油管的matlab仿真代码,基于状态模拟方法分析高压油管压力控制问题---源代码（matlab）...

1.m1.m

clear;clc;

load data1

y=E;

x=[ones(size(y)) p.^2 p];

[b,bint,r,rint,stats]=regress(y,x)

P=b(1)+b(2)*(p.^2)+b(3)*p;

figure;

plot(p,E,'r--',p,P,'b');

title("tan xing mo liang he ya qiang");

xlabel("ya qiang");

ylabel("tan xing mo liang");

legend('shi ji zhi','hui gui qu xian');

2.m2.m

clear;clc;

rho=0.85:0.005:0.90;

P=-53.23+227.0*tan(1.659+6.558*log(rho));

figure;

plot(rho,P);

title('???????');

xlabel('??(mg/mm^3)');

ylabel('??(MPa)');

3.m3.m

clear;clc;

load data1;

x=100:0.1:150;

y=1./(b(2).*x.^2+b(3).*x+b(1));

rho0=0.85;

rho=exp(trapz(x,y)+log(rho0))

x=100:0.1:160;

y=1./(b(2).*x.^2+b(3).*x+b(1));

rho0=0.85;

rho=exp(trapz(x,y)+log(rho0))

4.m4.m

clear;clc;

options=struct();

options.yfun=@yfun_1;

options.target=@target_1;

options.isplot=0;

options.interval=0.1;

options.duration=2000;

[fm,mx]=stepsolve(@targetfun,0.26,0.31,0.0005,0.02,5,options);

options.isplot=1;

targetfun(mx,options);

ylabel('??(MPa)');

xlabel('??(ms)');

title('???????????');

5.m5.m

clear;clc;

options=struct();

options.yfun=@yfun_1;

options.target=@target_2;

options.isplot=0;

options.interval=1;

options.duration=10000;

[fm,mx]=stepsolve(@targetfun,0.6,1,0.0005,0.05,5,options);

options.isplot=1;

targetfun(mx,options);

ylabel('??(MPa)');

xlabel('??(ms)');

title('???????????');

6.m6.m

clear;clc;

options=struct();

options.yfun=@yfun_2;

options.target=@target_1;

options.isplot=0;

options.interval=1;

options.duration=10000;

[fm,mx]=stepsolve(@targetfun,0.019,0.029,0.000001,0.001,3,options);

options.isplot=1;

targetfun(mx,options);

ylabel('??(MPa)');

xlabel('??(ms)');

title('???????????');

7.m7.m

clear;clc;

options=struct();

options.yfun=@yfun_3;

options.target=@target_1;

options.isplot=0;

options.interval=1;

options.duration=10000;

[fm,mx]=stepsolve(@targetfun,0.039,0.059,0.00001,0.002,3,options);

options.isplot=1;

targetfun(mx,options);

ylabel('高压油管压力(MPa)');

xlabel('运行时间(ms)');

title('两个喷油口工作时高压油管内部压力');

8.m8.m

clear;clc;

options=struct();

options.yfun=@yfun_4;

options.target=@target_1;

options.isplot=0;

options.interval=1;

options.duration=10000;

[fm,mx]=stepsolve(@targetfun,97,103,0.1,1,3,options);

options.isplot=1;

targetfun(mx,options);

ylabel('高压油管压力(MPa)');

xlabel('运行时间(ms)');

title('增加减压阀后高压油管内部压力');

9.Pfun.m

function P=Pfun(x)

if (x<0.7)

P=0;

else

P=(-53.23+227.0*tan(1.659+6.558*log(x)));

end

end

10.Q_in_1.m

function Q=Q_in_1(tA,t,P_in)

%tA????

%t????

%P_in

T=tA+10;

if mod(t,T)

Q=0.85*(pi/4)*(1.4)^2*sqrt(2*(160-P_in)/0.871);

else

Q=0;

end

end

11.Q_in_2.m

function Q=Q_in_2(w,t,P_in)

%tA????

%t????

%P_in

T=tA+10;

if mod(t,T)

Q=0.85*(pi/4)*(1.4)^2*sqrt(2*(160-P_in)/0.871);

else

Q=0;

end

end

12.Q_out_1.m

function Q=Q_out_1(t,P)

r=mod(t,100);

if r<=0.2

Q=100*r;

elseif r<=2.2

Q=20;

elseif r<=2.4

Q=-100*r+240;

else

Q=0;

end

end

13.Q_out_2.m

function Q=Q_out_2(t,p)

%t:时间

%p:密度

r=mod(t,100);

b1=[16.3667705292548,-2.27403710032532,0.0443842537666548];

b2=[16.4024620318978,-78.0705445458018,92.8628884740531];

if r<=0.44

h=polyval(b1,r);

elseif r<=2

h=2;

elseif r<=2.45

h=polyval(b2,r);

else

h=0;

end

dk=1.4;

dz=2.5;

alpha=9/180*pi;

A=min([pi*dk^2/4,pi/4*(4*dz*h*tan(alpha)+2*h^2*tan(alpha)^2)]);

Q=0.85*A*(2*Pfun(p)/p)^0.5;

end

14.Q_out_3.m

function Q= Q_out_3(p,rho)

Q=0.85*0.7^2*pi*sqrt(p/rho);

end

15.stepsolve.m

function [mf,M]=stepsolve(fun,L,R,tol,h0,d,options)

if nargin<=6

options=struct();

end

optfun=@(x)fun(x,options);

M=inf;

fprintf("Initiating...\n");

mf=optfun(M);

while R-L>tol

fprintf("Calculating.. From %.5f to %.5f step %.5f\n",L,R,h0);

t_A=L:h0:R;

ff=arrayfun(optfun,t_A);

mf=min(ff);

M=t_A(mf==ff);

L=M-(M-L)/d;

R=M+(R-M)/d;

h0=h0/d;

fprintf("Current: Min: %.5f At %.5f\n",mf,M);

end

end

16.target_1.m

function [f]=target_1(TT,PP)

f=norm(PP-100,2)/(max(TT)-min(TT));

end

17.target_2.m

function [f]=target_2(TT,PP)

if max(PP)>150

f=sum(diff(TT).*PP(1:end-1));

else

f=1e12;

end

end

18.targetfun.m

function f=targetfun(w,options)

if nargin==1

options=struct();

end

interval=initfun(options,'interval',1);

duration=initfun(options,'duration',20000);

isplot=initfun(options,'isplot',1);

target=options.target;

yfunopt=options.yfun;

m=0;

global flag;

flag=0;

yfun=@(t,y,m)yfunopt(t,y,m,w,interval);

y=0.850;

steps=0:interval:duration;

TT=zeros(size(steps));

YY=zeros(size(steps));

count=0;

for t=steps

count=count+1;

[y,m]=yfun(t,y,m);

TT(count)=t;

YY(count)=y;

end

PP=Pfun(YY);

f=target(TT,PP);

if isplot

figure;

plot(TT,PP,'r');

xlim([min(TT),max(TT)]);

ylim([min(PP),max(PP)]);

end

end

19.V_in_2.m

function V=V_in_2(w,t)

A=2.5^2*pi;

V0=20+2.413*A*2;

rho=-2.413*cos(w*t)+2.413;

V=V0-rho*A;

end

20.yfun_1.m

function [y,m]=yfun_1(t,y,m,w,dt)

Q_in=@Q_in_1;

Q_out=@Q_out_1;

V_0=500*(pi/4)*(10)^2;

y=y+((1/V_0)*(0.87*Q_in(w,t,Pfun(y))-y*Q_out(t,Pfun(y))))*dt;

end

21.yfun_2.m

function [y,m]=yfun_2(t,y,m,w,dt)

if (mod(t,2*pi/w)

m=V_in_2(w,0)*0.850;

end

V_0=500*(pi/4)*(10)^2;

V_A=V_in_2(w,t);

rho_A=m/V_A;

m_out=Q_out_2(t,y)*y;

P_A=Pfun(rho_A);

P_N=Pfun(y);

if (P_A>P_N)

Q_A=0.85*0.7*0.7*3.1415*sqrt(2*(P_A-P_N)/rho_A);

y=y+(1/V_0)*(rho_A*Q_A-m_out)*dt;

m=m-rho_A*Q_A*dt;

if (m<0)

m=0;

end

else

if (rho_A==0)

pause;

end

y=y+(1/V_0)*(-m_out)*dt;

end

end

22.yfun_3.m

function [y,m]=yfun_3(t,y,m,w,dt)

if (mod(t,2*pi/w)

m=V_in_2(w,0)*0.850;

end

V_0=500*(pi/4)*(10)^2;

V_A=V_in_2(w,t);

rho_A=m/V_A;

m_out=Q_out_2(t,y)*y*2;

P_A=Pfun(rho_A);

P_N=Pfun(y);

if (P_A>P_N)

Q_A=0.85*0.7*0.7*3.1415*sqrt(2*(P_A-P_N)/rho_A);

y=y+(1/V_0)*(rho_A*Q_A-m_out)*dt;

m=m-rho_A*Q_A*dt;

if (m<0)

m=0;

end

else

if (rho_A==0)

pause;

end

y=y+(1/V_0)*(-m_out)*dt;

end

end

23.yfun_4.m

function [y,m]=yfun_4(t,y,m,P_L,dt)

global flag;

w=0.05045;

if (mod(t,2*pi/w)

m=V_in_2(w,0)*0.850;

end

V_0=500*(pi/4)*(10)^2;%高压油管容积

V_A=V_in_2(w,t);%高压油泵内体积

rho_A=m/V_A;%高压油泵内密度

m_out_1=0;%喷油嘴堵塞

m_out_2=Q_out_3(Pfun(y),y)*y;%减压阀出油质量

m_out_3=Q_out_2(t,y)*y*2;%喷油嘴修复后出油质量

P_A=Pfun(rho_A);%高压油泵内压力

P_N=Pfun(y);%高压油管内压力

P_H=180;

if(P_N>P_H)

flag=1;

end

if (flag==0)%喷油嘴堵塞

if (P_A>P_N)

Q_A=0.85*0.7*0.7*pi*sqrt(2*(P_A-P_N)/rho_A);

P_A=Pfun(rho_A);

y=y+(1/V_0)*(rho_A*Q_A-m_out_1)*dt;

P_N=Pfun(y);

m=m-rho_A*Q_A*dt;

if (m<0)

m=0;

end

else

y=y+(1/V_0)*(-m_out_1)*dt;

end

end

%喷油嘴和减压阀同时打开

if (flag==1&&P_N>=P_L)

if (P_A>P_N)

Q_A=0.85*0.7*0.7*pi*sqrt(2*(P_A-P_N)/rho_A);

y=y+(1/V_0)*(rho_A*Q_A-m_out_2-m_out_3)*dt;

m=m-rho_A*Q_A*dt;

if (m<0)

m=0;

end

else

y=y+(1/V_0)*(-m_out_2-m_out_3)*dt;

end

end

%减压阀关闭

if (P_N

if (P_A>P_N)

Q_A=0.85*0.7*0.7*pi*sqrt(2*(P_A-P_N)/rho_A);

y=y+(1/V_0)*(rho_A*Q_A-m_out_3)*dt;

m=m-rho_A*Q_A*dt;

if (m<0)

m=0;

end

else

y=y+(1/V_0)*(-m_out_3)*dt;

end

end

end