2019考研数学一，2019考研数学一真题及答案
高压油管的压力控制赛题解析及代码实现
问题一题解
分析
要使高压油管内的燃油压力稳定在给定值100MPa，我们的切入点应该是将压力转换成可直接计算的量。附录三中已经给出了已知燃油的压力与弹性模量的关系，我们又已知燃油的压力变化量与密度变化量满足一微分方程，根据这两个条件我们可以得到燃油压力与密度的关系。再根据进出口的流量函数，我们可以推出计算出油管内燃油密度的变化规律。到此，我们就可以将原问题转化为对燃油密度优化的问题了。
（关于求解微分方程这一步骤，我觉得直接近似求解和拟合都可以，但是拟合的话最好不要用线性，二次和指数的效果相对可能好一丢丢。）
问题一的第二小问，可以在第一小问的工作上求解。第一小问是优化因变量压力，第二小问则是优化自变量时间，我们需要系统达到稳定的时间精确控制在指定数值左右。
求解
求解微分方程
因为燃油的压力变化量与密度变化量成正比，易知微分方程：
由附录三知
根据已知条件当压力为100 MPa时，燃油的密度为 0.850mg/mm^3解得
求解管内燃油密度变化
根据流量计算公式
可以得到燃油进出速度
燃油进入时的速度(一个周期内)
燃油喷出时的速度(一个周期内)
油管增加的质量
油管减少的质量
油管内质量的变化
油管内密度的变化
所以给定单向阀开启时长的值，我们就能通过燃油进出流量的函数来求得油管内质量的变化，进而求得密度的变化。再根据密度与压力的关系函数，我们就能求得压力的变化。当我们取很小的时间间隔为步长时，压力可以取离散值，那么原问题可以看作优化
因为单向阀每次开启结束后需要关闭10ms，考虑到周期的数量级，我们先将最优解的查找范围限定在(0.10ms)之间，步长为0.5ms。取100s中后5s压力稳定值的平均值作为结果，画出图像如图1所示。
由图1可以看出压力稳定在100Mpa时时间区间应该是(0.27ms,0.29ms)，我们在此区间再次搜索最优解。步长为0.01ms。结果如图2所示。
由图2可以进一步缩小区间至(0.2875ms,0.2876ms)，我们在此区间再次搜索最优解。步长为0.00001ms。结果如图3所示。
得到最后结果0.28752ms。对应压力变化图如图4所示。
问题二题解
分析
问题二与问题一一样都是压力控制，不同的是我们需要从几何图形上入手。进油时，在柱塞的压油过程中，柱塞被凸轮驱动上下运动，从而改变柱塞腔内的压力，决定高压油泵内的燃油是否进入。那么，由凸轮边缘曲线与角度的关系可以求得柱塞的运动情况，从而得到腔内的体积变化情况，进而求得腔内燃油的密度变化，由燃油压力与密度关系式得到其压力变化。喷油时，由已知喷油器的喷嘴结构及其工作条件和附录可以得到针阀与密封座之间的空隙面积。将问题一中给定的喷油器工作次数、高压油管尺寸和初始压力等条件代入模型求解，搜索凸轮的角速度，使高压油管内的压力尽量稳定在100 MPa左右。
求解
燃油的进入
由附件1给出的凸轮轮廓线极径、极角的关系，用插值法处理数据，经过三次差值后得到拟合轮廓线。对于凸轮的运动规律，可由对拟合曲线进行反解得到。在柱塞的运动过程中，我们关注其升程的最大值和最小值。其对应于凸轮轮廓线变化过程中的最高点与最低点。通过观察凸轮轮廓线发现，两个最值对应于凸轮轮廓线上的两端圆弧。以水平和竖直方向建立直角坐标系，将柱塞与轮廓接触点标为M，则该点对应坐标由式(18)计算
图8为凸轮与柱塞连接驱动的几何示意图，M点为某时刻凸轮轮廓线上y坐标值最大的点。
燃油的喷出
对喷油器喷嘴的示意图分析，将针阀与密封座之间的空隙面积记为S1，密封座底部面积记为S2。图9、图10分别为喷油嘴的的纵向截面和横向截面的各部分面积示意图。
根据上图的几何关系，密封座截面半径为：
对应的针阀与密封座之间的面积S1为：
同理问题一，当我们取很小的时间间隔为步长时，压力可以取离散值，那么原问题可以看作优化
求解该问题得到凸轮的角速度。
模型的求解
以0.01为步长搜索凸轮运动角速度，反复尝试使管内压力波动尽可能小。求得凸轮角速度为0.027249rad/ms时，高压油管压力稳定在100MPa左右。压力时间曲线如图12所示。
问题三题解
分析
问题三在问题二的基础上增加了一个喷油嘴，除此之外均与问题二相同，增加喷油嘴可以更快地将供油排出，从而减少高压油管内压力变化峰值，使油管内压力更加稳定。因为两个喷油嘴的喷油规律相同，所以我们可以直接在问题二基础上进行下一步工作。在原有模型的基础上考虑单向减压阀的引入。由单向减压阀的工作原理，设置单向减压阀的开启条件，设置阈值压强,管内压强大于阈值压强为减压阀开启条件，即此时燃油从油管内回流进减压阀，由流量计算公式求解减压阀中回流的燃油流量，从而得到高压油管内的燃油密度，进而得其压力，得到压力波动最小时的最优解。
求解
若要使高压油管内的压力稳定在100MPa左右，需调整喷油和供油策略。通过改变凸轮转速可以调整燃油的进入，通过单向减压阀调整高压油管内的压力；通过控制两个喷油嘴的喷油时间和调整针阀下止点来改变喷油策略。
要寻求合适的供油策略，相当于寻找一个合适的凸轮转动角速度ω，在该部分的求解中，我们将ω的值与其他参量联合起来进行求解。
喷油策略调整的主要目标为调节两个喷油嘴的喷油时间。这里假设每个喷油嘴的喷油规律与问题2中的相同，即1s喷射10次。
由于两喷油嘴喷油存在时间间隔，故设两喷油嘴喷油时间间隔为Δ t ，喷油嘴1在t=0时刻开始喷射，喷油嘴2在t = Δt时刻开始喷射。
搜索两个喷油嘴的喷油间隙。当喷油嘴开始喷射时间间隔在0-50ms内变化时，取特征点12.5ms，25ms 37.5ms，50.0ms。其压力变化如图13-16所示。
由图可得，当喷油嘴开始喷射时间间隔在0-50ms内变化时，管内压力在稳定值附近的波动幅度逐渐减小借由前面的分析，可知道0-50ms时间间隔的大致趋势，我们画出0-100ms时压强差指标的图。如图17所示。
从图中可以看到，延迟时间在63ms附近压力差最小。
在上述调整策略的基础上引入单向减压阀，当高压油管内压力过高时降低压力。使用单向减压阀时，需设置一个压力阈值P ′
当管内压力大于阈值时开启减压阀，使得管内燃油回流至外部低压油路，从而减小管内压力。设回流的燃油流量为Q d
，单向减压阀单次开启时长为t d ，则
当油管内压力小于阈值P’时，单向减压阀关闭。因此需求凸轮旋转角速度ω 和压力阈值P ′ ，最小化压力变化，即
最终找到的最优解为角速度0.0545，Δ t = 63 m s ,对应压力变化如图18：
代码附录
Q1.1
%密度函数，压力函数，出油函数
P=@(x)(x-0.80853816)/(0.001859637768-0.00170*x);
density=@(x)0.80853816*(1+ (0.6*10^-3*x)/(1+0.0017*x));
O=@(x)100*x.*(x<=0.2)+20.*(x>0.2&x<=2.2)-(100*x-240).*(x>2.2&x<=2.4)+0.*(x>2.4);
%规定常量
delta_t=0.01;
C=0.85;
V=(5^2)*500*pi;
A=(1.4/2)^2*pi;
p0=100;
m0=density(100)*V;
ml=m0;
y_test=[]
x_test=[]
l=0;
for t=0.28752:0.28752
result=0;
t_out=0;t_in=0;flag=0;
x=[];y=[]; %画图坐标
l=l+1;
for i=0:delta_t:100000
flag=flag+1;
x(flag)=i;
y(flag)=p0;
vin=0;
if t_in>t+10
t_in=t_in-(t+10);
end
if t_in
vin=(C*A*(2*(160-p0)/0.8696)^(.5));
end
t_in=t_in+delta_t;
m_in=0.8696*vin*delta_t;
if t_out>100
t_out=t_out-100;
end
vout=O(t_out)*delta_t;
t_out=t_out+delta_t;
mout=vout*density(p0) ;
delta_m=m_in-mout;
rho_now=(m0+delta_m)/V;
p0=P(rho_now);
m0=m0+delta_m;
flag=flag+1;
x(flag)=i;
y(flag)=p0;
end
%搜索判定条件
y_ave=0;k=0;
for j=9950000:10000000
k=k+1;
y_ave=y_ave+y(j);
end
y_ave=y_ave/k;
y_test(l)=y_ave;
x_test(l)=t;
end
plot(x,y)
%plot(x_test,y_test);
Q1.2
%密度函数，压力函数，出油函数
P=@(x)(x-0.80853816)/(0.001859637768-0.00170*x);
density=@(x)0.80853816*(1+ (0.6*10^-3*x)/(1+0.0017*x));
O=@(x)100*x.*(x<=0.2)+20.*(x>0.2&x<=2.2)-(100*x-240).*(x>2.2&x<=2.4)+0.*(x>2.4);
%规定常量
delta_t=0.01;
C=0.85;
V=(5^2)*500*pi;
A=(1.4/2)^2*pi;
p0=100;
m0=density(100)*V;
ml=m0;
y_test=[]
x_test=[]
l=0;
t=0;
result=0;
t_out=0;t_in=0;flag=0;
x=[];y=[]; %画图坐标
l=l+1;
for i=0:delta_t:20000
flag=flag+1;
x(flag)=i;
y(flag)=p0;
%if i<5000
% t=0.0001*i+0.26485;
%end
%if i>=5000
% t=0.76485;
%end
%if i<2000
% t=-0.0001*i+0.96685;
%end
%if i>=2000
% t=0.76485;
%end
if i<10000
t=0.00004*i+0.36485;
end
if i>=10000
t=0.76485;
end
vin=0;
if t_in>t+10
t_in=t_in-(t+10);
end
if t_in
vin=(C*A*(2*(160-p0)/0.8696)^(.5));
end
t_in=t_in+delta_t;
m_in=0.8696*vin*delta_t;
if t_out>100
t_out=t_out-100;
end
vout=O(t_out)*delta_t;
t_out=t_out+delta_t;
mout=vout*density(p0) ;
delta_m=m_in-mout;
rho_now=(m0+delta_m)/V;
p0=P(rho_now);
m0=m0+delta_m;
flag=flag+1;
x(flag)=i;
y(flag)=p0;
end
%搜索判定条件
%y_ave=0;k=0;
%for j=9950000:10000000
% k=k+1;
% y_ave=y_ave+y(j);
%end
%y_ave=y_ave/k;
%y_test(l)=y_ave;
%x_test(l)=t;
plot(x,y)
%plot(x_test,y_test);
Q2
%密度函数，压力函数
density=@(x)0.80853816*(1+ (0.6*10^-3*x)/(1+1.7*10^-3*x));
P=@(x)-(90071992547409920000*x-72826643121816530000)/(153122387330596864*x-167501279180178019);
c=0.85;
V=(5^2)*500*pi;
A=(1.4/2)^2*pi;
A_oil=(5/2)^2*pi;
H=5.8446;
Vmax=H*A_oil;
M=Vmax*density(0.5);
m0=density(100)*V;
delta_t=0.01;
p_wai=0.1;
theta=xlsread(‘附件1-凸轮边缘曲线’,’sheet1′,’A2:A629′);
r=xlsread(‘附件1-凸轮边缘曲线’,’sheet1′,’B2:B629′);
L0=0.7/tan(9/180*pi)*2.5/1.4;
A_zhen=(2.5/2)^2*pi;
A_kong=0.7^2*pi;
k=[];x=[];y=[];
al=xlsread(‘附件2-针阀运动曲线’,’sheet1′,’B2:B46′);
bl=2*ones(156,1);
cl=xlsread(‘附件2-针阀运动曲线’,’sheet1′,’E2:E46′);
h_zhen=[al;bl;cl;zeros(9755,1)]’;
for i=1:628
x(i)=sin(theta(i))*r(i);
y(i)=cos(theta(i))*r(i);
end
s=[];u=1;
for i=0:0.01:6.27
xk=[];
yk=[];
for j=1:628
xk(j)=x(j)*cos(i)-y(j)*sin(i);
yk(j)=x(j)*sin(i)-y(j)*cos(i);
end
s(u)=max(yk)-2.4130;
u=u+1;
end
s=interp1([0:0.01:6.27],s,[0:0.000001:6.27],’spline’);
w=0.027249;
j=1;
u=1;
p0=100;
rho0=0.85;
m_oil=M;
fy=1;
for i=0:delta_t:100000
fy=fy+int32(w*delta_t/0.000001);
while fy>6270001
fy=fy-6270001;
m_oil=M;
end
h_now=H-s(fy);
V_now=A_oil*h_now;
rho_now=m_oil/V_now;
P_now=P(rho_now);
vin=0;
if P_now>p0
vin=c*A*(2*(P_now-p0)/rho_now)^(.5);
end
m_in=rho_now*vin*delta_t;
m_oil=m_oil-m_in;
rl=(h_zhen(u)+L0)*tan(9/180*pi);h_z=h_zhen(u);
A_out=pi*rl^2-A_zhen;
if pi*rl^2-A_zhen-A_kong>0
A_out=A_kong;
end
if h_zhen(u)==0
A_out=0;
end
u=u+1;
if u>10001
u=u-10001;
end
vout=c*A_out*(2*(p0-p_wai)/rho0)^(.5);
mout=vout*rho0*delta_t;
delta_m=m_in-mout;
rho0=(m0+delta_m)/V;
k(j)=p0;
p0=P(rho0);
j=j+1;
m0=m0+delta_m;
end
plot(k)
xlabel(‘时间(t/ms)’)
ylabel(‘压力(P/MPa)’)
Q3
%密度函数，压力函数
density=@(x)0.80853816*(1+ (0.6*10^-3*x)/(1+1.7*10^-3*x));
P=@(x)-(90071992547409920000*x-72826643121816530000)/(153122387330596864*x-167501279180178019);
c=0.85;
V=(5^2)*500*pi;
A=(1.4/2)^2*pi;
A_oil=(5/2)^2*pi;
H=5.8446;
Vmax=H*A_oil;
M=Vmax*density(0.5);
m0=density(100)*V;
delta_t=0.01;
p_wai=0.1;
theta=xlsread(‘附件1-凸轮边缘曲线’,’sheet1′,’A2:A629′);
r=xlsread(‘附件1-凸轮边缘曲线’,’sheet1′,’B2:B629′);
L0=0.7/tan(9/180*pi)*2.5/1.4;
A_zhen=(2.5/2)^2*pi;
A_kong=0.7^2*pi;
k=[];x=[];y=[];
al=xlsread(‘附件2-针阀运动曲线’,’sheet1′,’B2:B46′);
bl=2*ones(156,1);
cl=xlsread(‘附件2-针阀运动曲线’,’sheet1′,’E2:E46′);
h_zhen=[al;bl;cl;zeros(9755,1)]’;
yanchi=5000;
for i=10001:-1:yanchi+1
h_zhenl(i)=h_zhen(i-yanchi);
end
for i=1:1:yanchi
h_zhenl(i)=0;
end
for i=1:628
x(i)=sin(theta(i))*r(i);
y(i)=cos(theta(i))*r(i);
end
s=[];u=1;
for i=0:0.01:6.27
xk=[];
yk=[];
for j=1:628
xk(j)=x(j)*cos(i)-y(j)*sin(i);
yk(j)=x(j)*sin(i)-y(j)*cos(i);
end
s(u)=max(yk)-2.4130;
u=u+1;
end
s=interp1([0:0.01:6.27],s,[0:0.000001:6.27],’spline’);
w=0.0545;
j=1;
u=1;
p0=100;
rho0=0.85;
m_oil=M;
fy=1;
for i=0:delta_t:100000
fy=fy+int32(w*delta_t/0.000001);
while fy>6270001
fy=fy-6270001;
m_oil=M;
end
h_now=H-s(fy);
V_now=A_oil*h_now;
rho_now=m_oil/V_now;
P_now=P(rho_now);
vin=0;
if P_now>p0
vin=c*A*(2*(P_now-p0)/rho_now)^(.5);
end
m_in=rho_now*vin*delta_t;
m_oil=m_oil-m_in;
rl=(h_zhen(u)+L0)*tan(9/180*pi);
rl2=(h_zhenl(u)+L0)*tan(9/180*pi);
A_out=pi*rl^2-A_zhen;
if pi*rl^2-A_zhen-A_kong>0
A_out=A_kong;
end
if h_zhen(u)==0
A_out=0;
end
Al_out=pi*rl2^2-A_zhen;
if pi*rl2^2-A_zhen-A_kong>0
Al_out=A_kong;
end
if h_zhenl(u)==0
Al_out=0;
end
u=u+1;
if u>10001
u=u-10001;
end
vout=c*(A_out+Al_out)*(2*(p0-p_wai)/rho0)^(.5);
mout=vout*rho0*delta_t;
delta_m=m_in-mout;
rho0=(m0+delta_m)/V;
k(j)=p0;
p0=P(rho0);
j=j+1;
m0=m0+delta_m;
end
plot(k)
xlabel(‘时间(t/ms)’)
ylabel(‘压力(P/MPa)’)
官方培训：第七届“数维杯”数学建模夏令营倒计时11天
为了更好地服务参赛同学，提高参赛者的建模能力，暑期数维杯夏令营开班啦！数维杯数夏令营每年举办一届，已经连续举办七届，已成为高校数学建模培训示范性基地，属于系统性的数学建模学习，主要针对国赛美赛等其它区域赛，采用分班教学，知名专家教授授课，7天集训+3天模拟赛+赛后点评解析，并配有助教全程课后指导， 累计指导学生近1000余队（2800余人次），国赛荣获800余项省级以上数模奖项，参培学员获奖率超过80%，每期仅限招200人！
群内获取往年夏令营资料、授课内容、夏令营详情等最新资讯
好啦~数模国赛即将来临啦
数乐君继续给大家发福利啦
里面都是一些数学建模竞赛中
经常能用到的必备干货商品
可以根据自己的需求去购买哦
2019考研数学一(2019考研数学一参考及答案)