数模-高压油管
题目
问题1：
首先分析密度和压力传播的变化特点，在高压环境下，将压力的非均匀变化问题简化为均匀变化问题，然后进行模型的建立与求解：第一步是用最小二乘法对附件3中弹性模量与压力的对应关系进行拟合，得出二次多项式关系；第二步是根据注1中给出的密度变化量和压力变化量的正比例关系计算出附件2中压力值对应的密度值，再次拟合得到密度与压力的关系式；第三步是通过进出高压油管的流量、单向阀和喷油嘴的工作时间等变量建立高压油管燃油量的变化模型；第四步是分析题目中要求高压油管内压力稳定在规定值的条件，由质量守恒定律得出密度随时间变化的模型，再结合压力与密度的关系模型建立高压油管内压力波动模型；第五步是确定第一小问模型，以总时间内压力波动最小为目标函数，以密度随时间变化的模型为约束条件，通过多次循环遍历的方法求得单向阀开启时长为0.2875ms时，高压油管内的压力能够维持在100MPa；第六步是在第一小问的模型的基础上加以修改，将总时间分别固定为2s、5s、10s时，同样通过多次循环遍历的方法得到当单向阀开启时长为0.7369ms时，高压油管内的压力可以从100MPa变为150MPa，并最终维持在150MPa左右；
附件3的内容：
导入数据：
%引入数据
data_3=xlsread('附件3-弹性模量与压力.xlsx',1,'A2:B402');
Pression=data_3(:,1);%压力
elastic_Mo=data_3(:,2);%弹性模量
然后进行压力和密度的拟合
clear,clc,close
%程序初始化
drou=0.00001;%密度变化步长
rou1(1)=0.85;%压强(0-100)对应的密度
rou2(1)=0.85;%压强(100-200)对应的密度
p1(1)=100;%压强(0-100)
p2(1)=100;%压强(100-200)
n=1;%计数器
%%算法主体
while(p1(end)>=0)
n=n+1;
rou1(n)=rou1(n-1)-drou;
E=fun(p1(n-1));
dp=E/((rou1(n-1)+rou1(n))/2)*drou;
p1(n)=p1(n-1)-dp;
end
n=1;
while(p2(end)<=220)
n=n+1;
rou2(n)=rou2(n-1)+drou;
E=fun(p2(n-1));
dp=E/((rou2(n-1)+rou2(n))/2)*drou;
p2(n)=p2(n-1)+dp;
end
%%数据拼接 这里相当于把0-100个的压力和密度倒序输出来
for i=1:size(p1,2)
Result(1,i)=p1(size(p1,2)+1-i);
Result(2,i)=rou1(size(p1,2)+1-i);
end
%这里是因为都有压强为100时候的压力和密度数据
p2(:,1)=[];
rou2(:,1)=[];
Result=[Result [p2;rou2]];
%去除首尾无效数据
flag=[];
for k=1:length(Result)
if Result(1,k)<0 | Result(1,k)>200
flag=[flag,k];
end
end
Result(:,flag)=[];
clear dp drou E i n p1 p2 rou1 rou2;
A=Result(1,:);
B=Result(2,:);
plot(A,B)
%分段多项拟合的函数
function[E]=fun(P)
E=0.02893*P^2+3.007*P+1572;
end
极径与极角关系式的确定
由附件1中给定的凸轮边缘曲线中极径和极角的数据进行线性拟合，绘制极径与极角的变化关系图，得到拟合曲线如下图所示：
data_1=xlsread('附件1-凸轮边缘曲线.xlsx',1,'A2:B629');
x=data_1(:,1);%极角
y=data_1(:,2);%极径
第一题的第一小问代码：
clear;
tic;
%%数据载入
load P_rou.mat %P与Rou的对应关系
P_ref=Result(1,:);
Rou_ref=Result(2,:);
dtt=0.00001;%对应关系里的时间精度
%%数据初始化
%可调参数
T=0.288;%单向阀开启时间
P(1)=100;%初始压强
rou=0.85;%初始密度(调节压强的时候请同时条件密度)
dt=0.001;%时间步长（精度)
Tall=30000;%总模拟时间
t0=0;%喷油开始时刻
%其他参数
V=12500*pi;%油管体积
m=rou*V;%初始质量
P_I=160;%入口处的压强
rou_I=0.87143;%160MPa时的密度
C=0.85;%流量系数
A=pi*0.7*0.7;%小孔面积
Itimes=1;%进油次数
Etimes=1;%出油次数
I_prev=0;%前一时刻的进油量
E_prev=0;%前一时刻的出油量
%进油第一个周期的时长
Tin_start=0;%进油开始时刻
Tin_end=T;%进油结束时刻
%出油第一个周期的时长
Tout_start=t0;%第一次喷油开始的时刻
Tout_1=Tout_start+0.2;%出油速率第一段
Tout_2=Tout_1+2;%出油速率第二段
Tout_end=Tout_2+0.2;%出油速率第三段，完成一个周期内的出油工作
%%主程序
for i=2:1:Tall/dt %0-Tal1ms
%时间步进
t=(i-1)*dt;
%计算t时刻的燃油质量m(t)
%判断是否同时有进油和出油
if t>Tin_start && t<=Tin_end && t>Tout_start && t<=Tout_end
%同时有进油和出油
I_now=C*A*sqrt(2*(P_I-P(i-1))/rou_I);%当前进油流量
I=(I_now+I_prev)/2; %总的进油
I_prev=I_now;
%计算出油流量速率
if t<Tout_1
E_now=100*(t-t0-100*(Etimes-1));
E=(E_now+E_prev)/2; %总的出油
E_prev=E_now;
elseif t<Tout_2
E_now=20;
E=(E_now+E_prev)/2;
E_prev=E_now;
else
E_now=240-100*(t-t0-100*(Etimes-1));
E=(E_now+E_prev)/2;
E_prev=E_now;
end
elseif t>Tin_start && t<=Tin_end %只有进油
I_now=C*A*sqrt(2*(P_I-P(i-1))/rou_I);%当前进油流量
I=(I_now+I_prev)/2;
I_prev=I_now; %总的进油
E_now=0;
E=(E_now+E_prev)/2;
E_prev=E_now; %总的出油
%出油是否进入下一个周期
if abs(t-(Tout_start+100))<0.001
%按照时间序列逐步进行计算，所以当刚好达到要转变的那个时间点的时候，abs(t-(Tout_start+100))=0，此时判断条件为<0.001，然后进行新一轮周期的运算
Etimes=Etimes+1;%新一次出油周期即将开始
Tout_start=t;
Tout_1=Tout_start+0.2;%出油速率第一段
Tout_2=Tout_1+2;%出油速率第二段
Tout_end=Tout_2+0.2;%出油速率第三段,完成一个周期内的出油工作
end
elseif t>Tout_start && t<=Tout_end %只有出油
I_now=0;
I=(I_now+I_prev)/2;
I_prev=I_now;
if t<Tout_1
E_now=100*(t-t0-100*(Etimes-1));
elseif t<Tout_2
E_now=20;
else
E_now=240-100*(t-t0-100*(Etimes-1));
end
E=(E_now+E_prev)/2;
E_prev=E_now;
%进油是否进入下一个周期
if abs(t-(Tin_end+10))<0.001
Itimes=Itimes+1;%新一次进油周期即将开始
Tin_start=t;
Tin_end=t+T; %进油结束时刻
end
else%没有进油和出油
I_now=0;
I=(I_now+I_prev)/2;
I_prev=I_now;
E_now=0;
E=(E_now+E_prev)/2;
E_prev=E_now;
if abs(t-(Tin_end+10))<0.001
Itimes=Itimes+1;%新一次进油周期即将开始
Tin_start=t;
Tin_end=t+T;%进油结束时刻
end
if abs(t-(Tout_start+100))<0.001
Etimes=Etimes+1;%新一次出油周期即将开始
Tout_start=t;
Tout_1=Tout_start+0.2;%出油速率第一段
Tout_2=Tout_1+2; %出油速率第二段
Tout_end=Tout_2+0.2;%出油速率第三段,完成一个周期内的出油工作
end
end
m=m+dt*(rou_I*I-rou*E) ;
rou=m/V;
%计算t时刻的密度
P(i)=P_ref(int16(floor(rou/dtt)-Rou_ref(1)/dtt+1));%计算t时刻的压强 直接利用下标索引进行导出
end
clear A C dt dtt E E_now E_prev Etimes i I I_now I_prev ...
Itimes k m P_I P_ref P_rou Pe qw rou rou_I Rou_ref t T t0...
Tall Tin_end Tin_start Tout_1 Tout_2 Tout_end Tout_start V W Z
%%图形显示
plot(P);
toc
这里需要用到压强和密度的关系，也就是前面探讨密度和压强关系的代码
对于单向阀的0.288的求解可以通过遍历的方法
再进一步调整精度
代码为：
clc,clear,close
kkk=1;
for T=0.28:0.001:0.31
%%数据载入
load P_rou.mat %P与Rou的对应关系
P_ref=Result(1,:);
Rou_ref=Result(2,:);
dtt=0.00001;%对应关系里的时间精度
%%数据初始化
%可调参数
%T=0.288;%单向阀开启时间
P(1)=100;%初始压强
rou=0.85;%初始密度(调节压强的时候请同时条件密度)
dt=0.001;%时间步长（精度)
Tall=30000;%总模拟时间
t0=0;%喷油开始时刻
%其他参数
V=12500*pi;%油管体积
m=rou*V;%初始质量
P_I=160;%入口处的压强
rou_I=0.87143;%160MPa时的密度
C=0.85;%流量系数
A=pi*0.7*0.7;%小孔面积
Itimes=1;%进油次数
Etimes=1;%出油次数
I_prev=0;%前一时刻的进油量
E_prev=0;%前一时刻的出油量
%进油第一个周期的时长
Tin_start=0;%进油开始时刻
Tin_end=T;%进油结束时刻
%出油第一个周期的时长
Tout_start=t0;%第一次喷油开始的时刻
Tout_1=Tout_start+0.2;%出油速率第一段
Tout_2=Tout_1+2;%出油速率第二段
Tout_end=Tout_2+0.2;%出油速率第三段，完成一个周期内的出油工作
%%主程序
for i=2:1:Tall/dt %0-Tal1ms
%时间步进
t=(i-1)*dt;
%计算t时刻的燃油质量m(t)
%判断是否同时有进油和出油
if t>Tin_start && t<=Tin_end && t>Tout_start && t<=Tout_end
%同时有进油和出油
I_now=C*A*sqrt(2*(P_I-P(i-1))/rou_I);%当前进油流量
I=(I_now+I_prev)/2; %总的进油
I_prev=I_now;
%计算出油流量速率
if t<Tout_1
E_now=100*(t-t0-100*(Etimes-1));
E=(E_now+E_prev)/2; %总的出油
E_prev=E_now;
elseif t<Tout_2
E_now=20;
E=(E_now+E_prev)/2;
E_prev=E_now;
else
E_now=240-100*(t-t0-100*(Etimes-1));
E=(E_now+E_prev)/2;
E_prev=E_now;
end
elseif t>Tin_start && t<=Tin_end %只有进油
I_now=C*A*sqrt(2*(P_I-P(i-1))/rou_I);%当前进油流量
I=(I_now+I_prev)/2;
I_prev=I_now; %总的进油
E_now=0;
E=(E_now+E_prev)/2;
E_prev=E_now; %总的出油
%出油是否进入下一个周期
if abs(t-(Tout_start+100))<0.001
%按照时间序列逐步进行计算，所以当刚好达到要转变的那个时间点的时候，abs(t-(Tout_start+100))=0，此时判断条件为<0.001，然后进行新一轮周期的运算
Etimes=Etimes+1;%新一次出油周期即将开始
Tout_start=t;
Tout_1=Tout_start+0.2;%出油速率第一段
Tout_2=Tout_1+2;%出油速率第二段
Tout_end=Tout_2+0.2;%出油速率第三段,完成一个周期内的出油工作
end
elseif t>Tout_start && t<=Tout_end %只有出油
I_now=0;
I=(I_now+I_prev)/2;
I_prev=I_now;
if t<Tout_1
E_now=100*(t-t0-100*(Etimes-1));
elseif t<Tout_2
E_now=20;
else
E_now=240-100*(t-t0-100*(Etimes-1));
end
E=(E_now+E_prev)/2;
E_prev=E_now;
%进油是否进入下一个周期
if abs(t-(Tin_end+10))<0.001
Itimes=Itimes+1;%新一次进油周期即将开始
Tin_start=t;
Tin_end=t+T; %进油结束时刻
end
else%没有进油和出油
I_now=0;
I=(I_now+I_prev)/2;
I_prev=I_now;
E_now=0;
E=(E_now+E_prev)/2;
E_prev=E_now;
if abs(t-(Tin_end+10))<0.001
Itimes=Itimes+1;%新一次进油周期即将开始
Tin_start=t;
Tin_end=t+T;%进油结束时刻
end
if abs(t-(Tout_start+100))<0.001
Etimes=Etimes+1;%新一次出油周期即将开始
Tout_start=t;
Tout_1=Tout_start+0.2;%出油速率第一段
Tout_2=Tout_1+2; %出油速率第二段
Tout_end=Tout_2+0.2;%出油速率第三段,完成一个周期内的出油工作
end
end
m=m+dt*(rou_I*I-rou*E) ;
rou=m/V;
%计算t时刻的密度
P(i)=P_ref(int16(floor(rou/dtt)-Rou_ref(1)/dtt+1));%计算t时刻的压强 直接利用下标索引进行导出
end
clear A C dt dtt E E_now E_prev Etimes i I I_now I_prev ...
Itimes k m P_I P_ref P_rou Pe qw rou rou_I Rou_ref t T t0...
Tall Tin_end Tin_start Tout_1 Tout_2 Tout_end Tout_start V W Z
%%图形显示
error=(P-100).^2;
sum_error(kkk)=sum(error)/10000000;
kkk=kkk+1;
end
x=0.28:0.001:0.31;
plot(x,sum_error)
在上述的情况中。设置初始情况为0.1的时候，会报错，因为
当开启时长只有0.1的时候，因为时间太短了，所以喷油嘴工作的更多，会出现压力为负数的情况，而上面这个图能够把负数的压力求出来，是因为他们把压力和密度的对应关系求出来了，而我们只是导入了一些数据，而不是关系式。
第一小题的第二问
代码（初始压强为130MPa）：
clear;
tic;
%%数据载入
load P_rou.mat %P与Rou的对应关系
P_ref=Result(1,:);
Rou_ref=Result(2,:);
dtt=0.00001;%对应关系里的时间精度
%%数据初始化
%可调参数
T=0.288;%单向阀开启时间
P(1)=130;%初始压强
rou=0.8612;%初始密度(调节压强的时候请同时条件密度)
dt=0.001;%时间步长（精度)
Tall=30000;%总模拟时间
t0=0;%喷油开始时刻
%其他参数
V=12500*pi;%油管体积
m=rou*V;%初始质量
P_I=160;%入口处的压强
rou_I=0.87143;%160MPa时的密度
C=0.85;%流量系数
A=pi*0.7*0.7;%小孔面积
Itimes=1;%进油次数
Etimes=1;%出油次数
I_prev=0;%前一时刻的进油量
E_prev=0;%前一时刻的出油量
%进油第一个周期的时长
Tin_start=0;%进油开始时刻
Tin_end=T;%进油结束时刻
%出油第一个周期的时长
Tout_start=t0;%第一次喷油开始的时刻
Tout_1=Tout_start+0.2;%出油速率第一段
Tout_2=Tout_1+2;%出油速率第二段
Tout_end=Tout_2+0.2;%出油速率第三段，完成一个周期内的出油工作
%%主程序
for i=2:1:Tall/dt %0-Tal1ms
%时间步进
t=(i-1)*dt;
%计算t时刻的燃油质量m(t)
%判断是否同时有进油和出油
if t>Tin_start && t<=Tin_end && t>Tout_start && t<=Tout_end
%同时有进油和出油
I_now=C*A*sqrt(2*(P_I-P(i-1))/rou_I);%当前进油流量
I=(I_now+I_prev)/2; %总的进油
I_prev=I_now;
%计算出油流量速率
if t<Tout_1
E_now=100*(t-t0-100*(Etimes-1));
E=(E_now+E_prev)/2; %总的出油
E_prev=E_now;
elseif t<Tout_2
E_now=20;
E=(E_now+E_prev)/2;
E_prev=E_now;
else
E_now=240-100*(t-t0-100*(Etimes-1));
E=(E_now+E_prev)/2;
E_prev=E_now;
end
elseif t>Tin_start && t<=Tin_end %只有进油
I_now=C*A*sqrt(2*(P_I-P(i-1))/rou_I);%当前进油流量
I=(I_now+I_prev)/2;
I_prev=I_now; %总的进油
E_now=0;
E=(E_now+E_prev)/2;
E_prev=E_now; %总的出油
%出油是否进入下一个周期
if abs(t-(Tout_start+100))<0.001
%按照时间序列逐步进行计算，所以当刚好达到要转变的那个时间点的时候，abs(t-(Tout_start+100))=0，此时判断条件为<0.001，然后进行新一轮周期的运算
Etimes=Etimes+1;%新一次出油周期即将开始
Tout_start=t;
Tout_1=Tout_start+0.2;%出油速率第一段
Tout_2=Tout_1+2;%出油速率第二段
Tout_end=Tout_2+0.2;%出油速率第三段,完成一个周期内的出油工作
end
elseif t>Tout_start && t<=Tout_end %只有出油
I_now=0;
I=(I_now+I_prev)/2;
I_prev=I_now;
if t<Tout_1
E_now=100*(t-t0-100*(Etimes-1));
elseif t<Tout_2
E_now=20;
else
E_now=240-100*(t-t0-100*(Etimes-1));
end
E=(E_now+E_prev)/2;
E_prev=E_now;
%进油是否进入下一个周期
if abs(t-(Tin_end+10))<0.001
Itimes=Itimes+1;%新一次进油周期即将开始
Tin_start=t;
Tin_end=t+T; %进油结束时刻
end
else%没有进油和出油
I_now=0;
I=(I_now+I_prev)/2;
I_prev=I_now;
E_now=0;
E=(E_now+E_prev)/2;
E_prev=E_now;
if abs(t-(Tin_end+10))<0.001
Itimes=Itimes+1;%新一次进油周期即将开始
Tin_start=t;
Tin_end=t+T;%进油结束时刻
end
if abs(t-(Tout_start+100))<0.001
Etimes=Etimes+1;%新一次出油周期即将开始
Tout_start=t;
Tout_1=Tout_start+0.2;%出油速率第一段
Tout_2=Tout_1+2; %出油速率第二段
Tout_end=Tout_2+0.2;%出油速率第三段,完成一个周期内的出油工作
end
end
m=m+dt*(rou_I*I-rou*E) ;
rou=m/V;
%计算t时刻的密度
P(i)=P_ref(int16(floor(rou/dtt)-Rou_ref(1)/dtt+1));%计算t时刻的压强 直接利用下标索引进行导出
end
clear A C dt dtt E E_now E_prev Etimes i I I_now I_prev ...
Itimes k m P_I P_ref P_rou Pe qw rou rou_I Rou_ref t T t0...
Tall Tin_end Tin_start Tout_1 Tout_2 Tout_end Tout_start V W Z
%%图形显示
plot(P);
toc
先判断达到150MPa的时候，最好的阀门开启时间
代码：
clc,clear,close
kkk=1;
for T=0.74:0.001:0.78
%%数据载入
load P_rou.mat %P与Rou的对应关系
P_ref=Result(1,:);
Rou_ref=Result(2,:);
dtt=0.00001;%对应关系里的时间精度
%%数据初始化
%可调参数
%T=0.288;%单向阀开启时间
P(1)=100;%初始压强
rou=0.85;%初始密度(调节压强的时候请同时条件密度)
dt=0.001;%时间步长（精度)
Tall=30000;%总模拟时间
t0=0;%喷油开始时刻
%其他参数
V=12500*pi;%油管体积
m=rou*V;%初始质量
P_I=160;%入口处的压强
rou_I=0.87143;%160MPa时的密度
C=0.85;%流量系数
A=pi*0.7*0.7;%小孔面积
Itimes=1;%进油次数
Etimes=1;%出油次数
I_prev=0;%前一时刻的进油量
E_prev=0;%前一时刻的出油量
%进油第一个周期的时长
Tin_start=0;%进油开始时刻
Tin_end=T;%进油结束时刻
%出油第一个周期的时长
Tout_start=t0;%第一次喷油开始的时刻
Tout_1=Tout_start+0.2;%出油速率第一段
Tout_2=Tout_1+2;%出油速率第二段
Tout_end=Tout_2+0.2;%出油速率第三段，完成一个周期内的出油工作
%%主程序
for i=2:1:Tall/dt %0-Tal1ms
%时间步进
t=(i-1)*dt;
%计算t时刻的燃油质量m(t)
%判断是否同时有进油和出油
if t>Tin_start && t<=Tin_end && t>Tout_start && t<=Tout_end
%同时有进油和出油
I_now=C*A*sqrt(2*(P_I-P(i-1))/rou_I);%当前进油流量
I=(I_now+I_prev)/2; %总的进油
I_prev=I_now;
%计算出油流量速率
if t<Tout_1
E_now=100*(t-t0-100*(Etimes-1));
E=(E_now+E_prev)/2; %总的出油
E_prev=E_now;
elseif t<Tout_2
E_now=20;
E=(E_now+E_prev)/2;
E_prev=E_now;
else
E_now=240-100*(t-t0-100*(Etimes-1));
E=(E_now+E_prev)/2;
E_prev=E_now;
end
elseif t>Tin_start && t<=Tin_end %只有进油
I_now=C*A*sqrt(2*(P_I-P(i-1))/rou_I);%当前进油流量
I=(I_now+I_prev)/2;
I_prev=I_now; %总的进油
E_now=0;
E=(E_now+E_prev)/2;
E_prev=E_now; %总的出油
%出油是否进入下一个周期
if abs(t-(Tout_start+100))<0.001
%按照时间序列逐步进行计算，所以当刚好达到要转变的那个时间点的时候，abs(t-(Tout_start+100))=0，此时判断条件为<0.001，然后进行新一轮周期的运算
Etimes=Etimes+1;%新一次出油周期即将开始
Tout_start=t;
Tout_1=Tout_start+0.2;%出油速率第一段
Tout_2=Tout_1+2;%出油速率第二段
Tout_end=Tout_2+0.2;%出油速率第三段,完成一个周期内的出油工作
end
elseif t>Tout_start && t<=Tout_end %只有出油
I_now=0;
I=(I_now+I_prev)/2;
I_prev=I_now;
if t<Tout_1
E_now=100*(t-t0-100*(Etimes-1));
elseif t<Tout_2
E_now=20;
else
E_now=240-100*(t-t0-100*(Etimes-1));
end
E=(E_now+E_prev)/2;
E_prev=E_now;
%进油是否进入下一个周期
if abs(t-(Tin_end+10))<0.001
Itimes=Itimes+1;%新一次进油周期即将开始
Tin_start=t;
Tin_end=t+T; %进油结束时刻
end
else%没有进油和出油
I_now=0;
I=(I_now+I_prev)/2;
I_prev=I_now;
E_now=0;
E=(E_now+E_prev)/2;
E_prev=E_now;
if abs(t-(Tin_end+10))<0.001
Itimes=Itimes+1;%新一次进油周期即将开始
Tin_start=t;
Tin_end=t+T;%进油结束时刻
end
if abs(t-(Tout_start+100))<0.001
Etimes=Etimes+1;%新一次出油周期即将开始
Tout_start=t;
Tout_1=Tout_start+0.2;%出油速率第一段
Tout_2=Tout_1+2; %出油速率第二段
Tout_end=Tout_2+0.2;%出油速率第三段,完成一个周期内的出油工作
end
end
m=m+dt*(rou_I*I-rou*E) ;
rou=m/V;
%计算t时刻的密度
P(i)=P_ref(int16(floor(rou/dtt)-Rou_ref(1)/dtt+1));%计算t时刻的压强 直接利用下标索引进行导出
end
clear A C dt dtt E E_now E_prev Etimes i I I_now I_prev ...
Itimes k m P_I P_ref P_rou Pe qw rou rou_I Rou_ref t T t0...
Tall Tin_end Tin_start Tout_1 Tout_2 Tout_end Tout_start V W Z
%%图形显示
error=abs(P-150);
sum_error(kkk)=sum(error)/10000000;
kkk=kkk+1;
end
x=0.74:0.001:0.78;
plot(x,sum_error)
发现0.765的时候最好，所以将0.765带入
代码：
clear;
tic;
%%数据载入
load P_rou.mat %P与Rou的对应关系
P_ref=Result(1,:);
Rou_ref=Result(2,:);
dtt=0.00001;%对应关系里的时间精度
%%数据初始化
%可调参数
T=0.765;%单向阀开启时间
P(1)=100;%初始压强
rou=0.85;%初始密度(调节压强的时候请同时条件密度)
dt=0.001;%时间步长（精度)
Tall=30000;%总模拟时间
t0=0;%喷油开始时刻
%其他参数
V=12500*pi;%油管体积
m=rou*V;%初始质量
P_I=160;%入口处的压强
rou_I=0.87143;%160MPa时的密度
C=0.85;%流量系数
A=pi*0.7*0.7;%小孔面积
Itimes=1;%进油次数
Etimes=1;%出油次数
I_prev=0;%前一时刻的进油量
E_prev=0;%前一时刻的出油量
%进油第一个周期的时长
Tin_start=0;%进油开始时刻
Tin_end=T;%进油结束时刻
%出油第一个周期的时长
Tout_start=t0;%第一次喷油开始的时刻
Tout_1=Tout_start+0.2;%出油速率第一段
Tout_2=Tout_1+2;%出油速率第二段
Tout_end=Tout_2+0.2;%出油速率第三段，完成一个周期内的出油工作
%%主程序
for i=2:1:Tall/dt %0-Tal1ms
%时间步进
t=(i-1)*dt;
%计算t时刻的燃油质量m(t)
%判断是否同时有进油和出油
if t>Tin_start && t<=Tin_end && t>Tout_start && t<=Tout_end
%同时有进油和出油
I_now=C*A*sqrt(2*(P_I-P(i-1))/rou_I);%当前进油流量
I=(I_now+I_prev)/2; %总的进油
I_prev=I_now;
%计算出油流量速率
if t<Tout_1
E_now=100*(t-t0-100*(Etimes-1));
E=(E_now+E_prev)/2; %总的出油
E_prev=E_now;
elseif t<Tout_2
E_now=20;
E=(E_now+E_prev)/2;
E_prev=E_now;
else
E_now=240-100*(t-t0-100*(Etimes-1));
E=(E_now+E_prev)/2;
E_prev=E_now;
end
elseif t>Tin_start && t<=Tin_end %只有进油
I_now=C*A*sqrt(2*(P_I-P(i-1))/rou_I);%当前进油流量
I=(I_now+I_prev)/2;
I_prev=I_now; %总的进油
E_now=0;
E=(E_now+E_prev)/2;
E_prev=E_now; %总的出油
%出油是否进入下一个周期
if abs(t-(Tout_start+100))<0.001
%按照时间序列逐步进行计算，所以当刚好达到要转变的那个时间点的时候，abs(t-(Tout_start+100))=0，此时判断条件为<0.001，然后进行新一轮周期的运算
Etimes=Etimes+1;%新一次出油周期即将开始
Tout_start=t;
Tout_1=Tout_start+0.2;%出油速率第一段
Tout_2=Tout_1+2;%出油速率第二段
Tout_end=Tout_2+0.2;%出油速率第三段,完成一个周期内的出油工作
end
elseif t>Tout_start && t<=Tout_end %只有出油
I_now=0;
I=(I_now+I_prev)/2;
I_prev=I_now;
if t<Tout_1
E_now=100*(t-t0-100*(Etimes-1));
elseif t<Tout_2
E_now=20;
else
E_now=240-100*(t-t0-100*(Etimes-1));
end
E=(E_now+E_prev)/2;
E_prev=E_now;
%进油是否进入下一个周期
if abs(t-(Tin_end+10))<0.001
Itimes=Itimes+1;%新一次进油周期即将开始
Tin_start=t;
Tin_end=t+T; %进油结束时刻
end
else%没有进油和出油
I_now=0;
I=(I_now+I_prev)/2;
I_prev=I_now;
E_now=0;
E=(E_now+E_prev)/2;
E_prev=E_now;
if abs(t-(Tin_end+10))<0.001
Itimes=Itimes+1;%新一次进油周期即将开始
Tin_start=t;
Tin_end=t+T;%进油结束时刻
end
if abs(t-(Tout_start+100))<0.001
Etimes=Etimes+1;%新一次出油周期即将开始
Tout_start=t;
Tout_1=Tout_start+0.2;%出油速率第一段
Tout_2=Tout_1+2; %出油速率第二段
Tout_end=Tout_2+0.2;%出油速率第三段,完成一个周期内的出油工作
end
end
m=m+dt*(rou_I*I-rou*E) ;
rou=m/V;
%计算t时刻的密度
P(i)=P_ref(int16(floor(rou/dtt)-Rou_ref(1)/dtt+1));%计算t时刻的压强 直接利用下标索引进行导出
end
clear A C dt dtt E E_now E_prev Etimes i I I_now I_prev ...
Itimes k m P_I P_ref P_rou Pe qw rou rou_I Rou_ref t T t0...
Tall Tin_end Tin_start Tout_1 Tout_2 Tout_end Tout_start V W Z
%%图形显示
plot(P);
toc
后续不再介绍
转载请注明出处，欢迎讨论和交流!
