2019年高教社杯全国大学生数学建模竞赛A题代码和思路

基于数值模拟的高压油管压力控制模型

摘要

高压油管是燃油发动机供油系统中连接喷油泵与喷油器的重要零部件，其内部压力的变化会影响发动机的工作效率。对于给定参数规格的高压油管，本文通过建立数学模型，研究了不同进油、喷油方式下油管内压力变化情况。

对于问题一，首先根据题目给定的压力变化量与密度变化量的对应关系，通过数值积分的方式，得到了一定精度下，不同燃油密度对应的燃油压力。建立进出油过程中油管内部压力变化的微分方程，即高压油管内部压力变化的微分模型，在此基础上研究了油管系统随时间的波动情况和稳态性质。将单向阀开启时间作为决策变量，以稳态阶段油管燃油压力与期望压力偏差的平方和为目标函数，建立基于油管压力变化微分模型的优化模型，并将微分模型改写成差分形式，通过编程遍历求解得到最优单向阀开启时长为0.288ms。通过对模型稳态的研究，将期望压力设定为150MPa，结合上述微分模型，通过数值模拟求解得出，调整时间为 2s 的策略为 0 到 2s 内单向阀每次开启时长为 0.913ms，随后开启时长变为0.752ms; 调整时间为5s的策略为0到5s内单向阀每次开启时长为0.755ms,随后开启时长变为0.752ms； 调整时间为10s的策略为0到4.1s内单向阀每次开启时长为0.288ms, 随后开启时长变为0.752ms。

对于问题二，针对重新定义的高压油泵供油方式和喷油嘴的出油方式，首先研究柱塞腔内压力变化与凸轮转动的关系，建立柱塞腔内压力变化的微分模型。接着分析针阀运动，确定喷油嘴有效面积与针阀运动的关系。最后，建立油管系统压力变化的微分模型，并将其转化为差分形式，以凸轮转动角速度作为决策变量建立优化模型。通过编程遍历求解得出凸轮转动角速度的最优值为27.5rad/s。

对于问题三，将两个喷油嘴工作时间的间隔和凸轮转动角速度共同作为决策变量，重新定义出油量，建立了两个喷油嘴的高压油管系统二元最优化模型。通过数值模拟的方法得到凸轮角速度最优值为 55rad/s，喷油嘴工作时间间隔为50ms。进一步考虑安装减压阀后的情况，分别研究泄压阀式和电磁阀式的减压阀。对于含泄压阀的油管系统，将泄压阈值作为新的决策变量，改进上一小问的模型，通过数值模拟求解得到一个最优解为压力阈值取 100MPa， 凸轮角速度为440 rad/s。对于含电磁阀的油管系统，将电磁阀开启时长和关闭时长作为决策变量建立多元优化模型，通过数值模拟求解得到一个最优解为每 50ms 泄压阀工作2ms，且距离上一次喷油嘴出油间隔为25.225ms，此时凸轮角速度为110 rad/s。随后，对两种减压阀的减压效果做出评价。

最后，给出了模型的优缺点和改进方法。

关键词：常微分方程 差分法 数值模拟 数值积分 最优化模型

一、 问题重述

1  问题背景

高压油管是柴油发动机供油系统中连接喷油泵与喷油器的重要零部件，它作为发动机工作时的"大动脉"，承担着给发动机输送高压燃油的任务。燃油进入和喷出高压油管是许多燃油发动机工作的基础，其间歇性工作过程会导致高压油管内压力的变化，使得所喷出的燃油量出现偏差，从而影响发动机的工作效率。

2  需要解决的问题

某型号高压油管的内腔长度为500mm， 内直径为10mm，供油入口A处小孔的直径为1.4mm。

问题1：高压油泵在入口A处提供的压力恒为160 MPa，高压油管内的初始压力为100 MPa。喷油器每秒工作10次，每次喷油时间为2.4ms，已知其每一时刻的喷油速率，要求通过单向阀开关控制供油时间的长短，且单向阀每打开一次后就要关闭10ms，设置单向阀每次开启的时长使高压油管内的压力尽可能稳定在100 MPa。若要将高压油管内的压力从100 MPa增加到150 MPa, 且分别经过约2 s、5 s和10 s的调整过程后稳定在150 MPa，确定调整单向阀开启时长的策略。

问题2：在实际工作过程中，高压油管A处的燃油来自高压油泵的柱塞腔出口，喷油由喷油嘴的针阀控制。当凸轮运动使柱塞腔内的压力大于高压油管内的压力时，柱塞腔与高压油管连接的单向阀开启，燃油进入高压油管内。柱塞腔内直径为5mm，柱塞运动到上止点位置时，柱塞腔残余容积为 20mm³。柱塞运动到下止点时，低压燃油会充满柱塞腔(包括残余容积)，低压燃油的压力为0.5 MPa。喷油器针阀直径为2.5mm、密封座是半角为9°的圆锥，最下端喷孔的直径为1.4mm。针阀升程为0时，针阀关闭； 针阀升程大于0时，针阀开启，燃油向喷孔流动，通过喷孔喷出。在问题1中给出的喷油器工作次数、高压油管尺寸和初始压力下，确定凸轮的角速度，使得高压油管内的压力尽量稳定在100MPa左右。

问题3：在问题2的基础上，再增加一个喷油嘴，每个喷嘴喷油规律相同，喷油和供油策略应如何调整? 为了更有效地控制高压油管的压力，现计划在D处安装一个单向减压阀。单向减压阀出口为直径为1.4mm的圆，打开后高压油管内的燃油可以在压力下回流到外部低压油路中，从而使得高压油管内燃油的压力减小。请给出高压油泵和减压阀的控制方案。

二、问题分析

问题的研究对象是高压油管，研究内容为其中燃油的压力变化情况。该问题描述了该系统在不同的进油和出油方式下的压力变化特点，并在不同的进油和出油方式下提出了不同的要求。

1  问题1

第一小问：问题定义了恒定的出油方式，同时其定义的进油方式只和系统压力与单向阀每次开启时长相关。其次提出了“将高压油管内的压力尽可能稳定在100MPa左右”的要求。这首先要求我们对“尽可能稳定在100MPa附近”做出描述，其次由问题描述易知，系统唯一的可改变变量为单向阀开启时长， 这就要求我们将其作为决策变量，研究上述描述的最优化问题。

第二小问：问题要求我们在不同的时间内将系统压力从100MPa调整到150MPa并使其稳定。这首先要求我们先将前一小问中期望的压力设定为150MPa，研究如何设定单向阀开启时间使其稳定，然后另外确定一个单向阀开启时间使得在规定时间内， 系统压力能够从100MPa上升到150MPa。

2  问题2

该问题重新定义了进油出油方式，此时进油方式将与凸轮的角速度有关，这要求我们在问题一的基础上重新描述进出油方式，并将决策变量设定为凸轮的角速度。同时由于引进了柱塞腔与凸轮的概念，这要求我们将其作为一个与高压油管相关联的子系统，并对其做出描述。

3  问题3

第一小问：该问题在问题二的基础上，给系统新加了一个喷油嘴，由此出油方式将会发生改变，由于新加的喷油嘴的喷油时间是否与原来的时间重合未得到描述，所以这需要我们将这一问题作为一个决策变量。

第二小问：在前一小问的基础上，又定义了一个泄压阀。 由于问题未对其工作方式进行准确的描述，所以需要我们对泄压阀工作方式进行定义，并且将其作为一个决策因素。

三、模型假设与约定

1  高压油管是刚体，不考虑其形变；

2  高压油管的体积是理想圆柱体体积；

3  不考虑高压油管内壁对油的粘滞力；

4  不考虑液体压力波动和压力传播的时间，腔内液压处处相等；

5  整个系统温度恒定；

6  不考虑供油入口和喷油口的阀门开关的延迟时间。

四、符号说明及名词定义

P(t)	t时刻高压油管的压力
m(t)	t时刻高压油管储油总质量
p(t)	t时刻高压油管内燃油密度
I(t)	t时刻的进油速率
E(t)	t时刻的出油速率
Q₁	进油时的速率
Q₂	出油时的速率
T	单向阀的开启时间
T₁00	稳态为100MPa时， 单向阀开启时间
T₁50	稳态为150MPa时， 单向阀开启时间
te	期望的时间
Pe	期望的稳态压力
Plow	喷油嘴及减压阀外界的压力
tn	系统达到稳态所需的时间
[tn₁,tn₂]	系统稳态后的某段周期整数倍的时间
Tlcm	喷油周期与供油周期的最小公倍数的某一倍数
Pc, Pc, mc, Vc	柱塞腔内燃油密度， 压力， 质量， 体积
ω	凸轮转动的角速度
Tc	凸轮转动的周期
f₁(t), f₂(t)	进出油的状态变量
PR	泄压阈值压力

五、模型准备与性质分析

1  基本模型：高压油管压力变化的微分模型

将高压油管内的压力作为研究对象，从其变化的角度研究其进油与出油情况。

假设每次单向阀开启的时间为T； 记初始时刻为0，且此时单向阀开启； t₀时刻喷油嘴开始喷油， t₀∈097.6;

另设t时刻高压油管的压力为 P(t)、储油总质量为 m(t)、燃油密度为 ρ(t)；设高压油管的容积恒为 V=12500πmm3,ρI为 160MPa 下燃油的密度，其值恒定。

描述高压油泵及喷油嘴周期性所产生的工作进油时间段与出油时间段，记：进油时间段

 Ω1,1,Ω1,2,,Ω1,n1

 Q1,i=0+i-1⋅T+10T+i-1⋅T+10,i=1,2,,n

出油时间段

 Ω2,1,Ω2,2,,Ω2,n1

 22,i=t0+100⋅i-12.4+t0+100⋅i-1;i=1,2,,n2于是可以得到t时刻的进油速率I(t)和出油速率E(t)

                 It={Q1t,t∈Ω1,i},i=1,2,,n1          (1)

                         (2)

其中

                      Q1t=CA2APtρI                   (3)

                Q2t=100t',t'∈00.2240-100t',t'∈2.22            (4)

                     t'=t-t₀+100⋅i-1                (5)

考虑 t 时刻到t+ dt时刻高压油管中燃油的质量变化量，应有 mt+dt-mt=ρI⋅It-ρt⋅Etd于是

                      dmdt=ρI⋅I-ρt⋅E                  (6)

另考虑系统的初始时刻状态

                         m(0)=ρ(0)·V                     (7)

其中ρ(0)已知，其值为( 0.850mg/mm³;

综上所述，高压油管进出油模型为以下初值问题

                                      (8)

其中

                                Q₁(t)=C 22₂P₄₀                    (9)

Q₂(t)=(2020100220)24. Y-t- to+100.(-1)

0. x-10.0(-1)(-1-1));i=12.1=10+100(-1));i=12.1=10+100(-1)(1+100.

系统的波动和稳态

设系统中单位时间内液体压力的均值随时间变化情况的描述量为 hₚt。

                              hpt=t-T/2t+T/2 PtdtT                     (10)

总体上，考虑从系统初始状态到无限远的时间内的变化，由式(6) 可知，当供油端压力恒定时，若一段时间内高压油管的燃油整体压力升高，则相应的I将会减小，而整体出油速率不发生改变，所以 dmdt逐渐减小，即高压油管燃油质量增速放缓甚至呈负增长。而若高压油管的燃油整体压力下降，则相应的I将会增大，于是 dmdt逐渐增大，即高压油管燃油质量减速放缓甚至呈正增长。

所以无论系统初始状态如何，无论开阀时长多少，系统中单位时间内压力的均值随时间的增长将趋近于一个常数，即

 hp(t)→ 常数

当 hp(t)等于一个常数时，系统达到稳态。如图1所示。

此时系统应满足以下方程组：

                  dmdt=ρI⋅I-ρt⋅EmnTlcm=mn-1Tlcm,n>n0               (11)

其中 Tlem为喷油周期与供油周期的最小公倍数的某一倍数。

另考虑单个喷油周期，因为若要使高压油管的压力稳定在期望压力附近，其压力应当随时间在期望压力附近波动。而综合比较单向阀开启时间和喷油嘴开启时间，单次喷油量应远大于单次供油量。于是理想情况下单个喷油周期应有如图2的特征，当系统达到稳态时，单个喷油周期内，期望压力以下曲线面积与以上面积相等。

3  燃油压力P和燃油密度ρ关系的确定

进油与出油决定了高压油管内燃油的总质量变化，于是由密度计算公式可以确定燃油密度的变化情况

 ρt=mtV

根据题目注解1的内容，可以确定压力P(t)与密度ρ(t)之间的关系，记为

P = g(ρ)

设燃油的压力变化量为 dP，密度变化量为 dρ； 由于燃油的压力变化量与密度变化量成正比，比例系数为 Eρ,有

                          dP=Eρ⋅dρ                      (12)

其中E=E(P);

于是

                        ∫1EdP=∫1ρdρ+c                  (13)

求出上述不定积分，且利用条件： ρ=0.85mg/mm³时, P =100MPa,即可得到ρ与P的关系表达式。

但由于E=E(P)没有具体表达式，且该积分难以计算。根据积分的定义，我们对该公式进行了简单的数值积分，得到不同步长下ρ与P的值，如图1所示。具体代码与结果见文件“密度压力关系. xlsx”.

六、模型的建立与求解

1  问题1：简化结构下基于基本模型的规划模型

1.1 第一小问模型的建立与求解

1.1.1  建立规划模型

       版权归全国大学生数学建模竞赛所有，如有侵权请联系删除                                                 转载自中国大学生在线

题目要求设定高压油管中燃油的期望压力为 Pₑ=100MPa,设计单向阀开启时长T，使得一定时间里高压油管燃油实际压力的偏差达到最小； 假设研究的时间段为 t₁t₂;

于是设决策变量为

T

目标函数为

 minZ=t1t2 Pt-Pe2dt

约束条件有：

(1)  模型满足基本模型的方程组(8)(9)

(2)  模型满足燃油密度和燃油压力的关系P=g(ρ)

即

                        dmdt=ρI⋅I-ρt⋅Em0=ρ0.V                (14)

其中

 ρt=mtV

 Q1t=CA2ΔPtρt

 Q2t=100t',t'∈00.220,t'∈0.22240-100t',t'∈2.22.4t'=t-t0+100⋅i-1

 

1.1.2 求解规划模型

由于该模型难以从得到解析解的角度上进行求解，即方程组(14)难以写成显性的解析解。这里给出方程组(14) 的差分形式，并对其进行数值计算得到一列不同T取值下对应 Z 的数值，然后对[T，Z]进行多项式的拟合，即用多项式的形式对其实际关系 Z=ht进行描述，于是模型转化为无约束条件的优化问题，对其进行求解即可得到问题的近似解。

首先，设时间步长为 dt，以 dt为间隔将时间分为

 0=t0<t1<t2<⋯<tN

假设 tn时刻，高压油管的燃油压力为 Pₙ、燃油密度为( ρₙ、、燃油总质量为 mₙ、

进油速率为 In、出油速率为1

Eₙ。

 

表示系统达到稳态后的某段周期整数倍的时间。

这里具体数值计算由 Matlab实现，具体代码见附录2，程序流程如下：

1.1.3  模型的结果与检验

取期望压力 Pₑ=100MPa,通过遍历T得出以下公式所给出的偏差平方和 Z

                        Z=∑i=tn1tn2Pn-Pe2⋅dt                 (18)

程序模拟足够长的时间 50S，取最后一秒即 10 个喷油周期的以上偏差平方和Z，如图所示。

图中可以看出，当T =0.288ms时，高压油管内的压力与100 MPa偏差最小，最为稳定。

对于该结果的所取时间段是否达到稳态，用该段时间压力的偏差和 Z。进行检验：

                         Ze=∑i=tn1tn2Pn-Pe⋅dt                   (19)

当所取时间段，高压油管内压力达到稳定状态时，Z。的值应趋于0，检验结果如下：

30000   20000     y = -2E+06x² + 1E+06x- 249497 R² = 0.9998 10000 0         0.125              0.3               0.35 -10000 -20000

6000 y = -2E+06x² + 1E+06x- 263735   4000                  R² = 1 2000 0       0. 28    0.285      0.29    0.295      0.3 -2000 -4000

   图7 T=0.25:0.35时, 最后一秒内压力的偏差和     图8 T=0.28:0.3时, 最后一秒内压力的偏差和

可以看出当T =0.288ms，50S后压力偏差和近似为0，即此时达到稳态。

综上可以得出结论，当单向阀每次开启的时长T=0.288ms时，高压油管内的压力尽稳定在100 MPa左右，且偏差最小。

 

   版权归全国大学生数学建模竞赛所有，如有侵权请联系删除                                             转载自中国大学生在线

1.2 第二小问模型的建立与求解

1.2.1  稳态性质的进一步研究

在五.2中，我们论述了不同初始状态、不同开阀时长 T 的系统，最终都要到达稳态，压力值在一个常数附近波动。

进一步分析，猜想不同开阀时长对应不同的稳态压力值，即当T确定时，从不同的初始状态出发都会到达同一个稳态，不过到达稳态的时间不同。下面我们利用第一小问的模型, 绘制了T=0.288ms, 初始压力分别为100MPa 和130MPa的压力变化图。见图7、8。

由图示以及计算稳态时压力的均值，发现两个过程最后都在100MPa附近波动，证明猜想正确。

另外，易知到达稳态所用的时间与T有关，T越大，到达稳态越迅速。

这表明了不同T的取值将对应唯一的系统稳态的压力均值 hₚt;喷油开始时刻t₀和初始压力的取值，对最终系统状态无影响。

1.2.2 模型的分析

由模型稳态的描述，稳定时150MPa的期望压力应当对应唯一的单向阀开启时间。所以，可以先将模型期望压力设定为 150MPa，得到对应单向阀开启时长 T₁₅₀。然后研究该T₁₅₀时长下系统到达 150MPa 的稳定压力需要的时间。最后对2S、5S、10S的要求分别设计不同的调整方案。

1.2.3 模型的建立

在前一小问的基础上，将期望压力设定为 Pₑ=150MPa,初始状态不改变；期望时间为 te，喷油嘴开启时刻取 t₀=0;利用前一小问的模型，得到对应' T₁₅₀。利用一定时间内压力的均值 hₚt(如式(10))来描述整体压力，找到系统到达150MPa的稳定压力所需要的时间 tn。分以下情况设计调整方案。

3.1.2 模型的求解

模型的差分模型与程序与前一问基本一致，不再赘述，具体代码见附录。

由于该模型最终转化为一个无约束条件的多元规划问题，结合问题一所描述问题的性质与稳态与决策变量的唯一对应关系，可知该问题有唯一的全局最优解，且目标函数是下凸的，如下图所示，则可以通过以下方式找到最优解。

设目标函数为

 minZ=fx1x2xn

固定x₂，…，xn为任意值，于是问题转化为无约束条件下单变量的规划问题

 minZ=f₁₁x₁

找到该问题的最优解，例如  x₁=x₁₀; 固定 x₁=x₁₀, 然后固定x₃,…, xn, 求解

 minZ=f₁₂x₂

…………

固定 x1=x10,x2=x20,⋯,xn-1=xn-1,求解

 minZ=f₁ₙx₁

固定 x2=x20,⋯,xn=xn0,求解

 minZ=f₂₁x₁

如此往复，当重复次数足够多时，即可找到问题的近似最优解。

3.1.3 模型的结果

通过数值模拟的方法，固定 t₀'=50ms,取凸轮角速度ω为54.5到55.5rad/s,得到压力的偏差平方和以及偏差和随ω的变化情况如图，得到ω的最优取值为55rad/s

 

接着固定ω=55rad/s, 取凸轮角速度t0为0到95ms, 如图所示, t0的最优值为45ms。

最后固定 t₀'=45ms,取不同ω的值，重复以上过程，得到问题的近似最优解为

 t₀'=50ms

ω= 55rad/s

代码

clear;

2. tic;             %程序计时

3. %数据载入

4. load ' Result';     %密度压强关系

5. load ' rad';       %凸轮极径

6. load 'S'        %升程与小孔面积关系

7. %以上S 的来源

8. % load '针阀运动. mat';

9. % 50= pi*0.7^2;

10. %%生成针阀运动对应的出油口面积

11. % for i=1:1:246 %0-2.45ms

12. %   S(i)= pi*f(2,i)* sind(9)*(2.5+f(2,i)* sind(9)* cosd(9));

13. %     if   S(i)>S0

14. %       S(i)=S0;

  15. %    end

16. % end

17. %可调节参数

18. w=27.5;         %转速，弧度制

19. dt=0.01;        %精度，时间步长

20. Talls=50;       %总运行时间，秒

21. Tall= Talls*1000;%总运行时间, 毫秒

22. rou=0.85;       %初始密度(程序对该数据精度敏感)

23. t0=0;           %喷油开始时间

24. pc=0.5;         %初始油泵压强

25. %恒常参数及衍生参数

26. v=12500* pi;    %油管体积

    27.m=rou*√;        %油管油量

28. v0=20;          %油泵残余油量

 29.vc=vΘ+pi*2.5*2.5*7.239-2.413;%油泵最大油量

   30. rc=2.413;       %升程

31. rouc=0.80255;       %油泵密度

32. mc= rouc* vc;       %油泵质量

33. p(1)= Result(1, floor(10^5*( rou)-80229));%油管初始压强

34. %程序用参数

35. pt=p;

36. for i=1:1:( Tall/ dt)

37.     p(i)= pt;

38.     t=i* dt;           %当前时间, ms

39.     %循环累积量初始化 (程序用变量)

40.       dQ=0;

41.       dm=0;

42.       dQc=0;

43.      dvc=0;

44.       dmc=0;

45.     %%变化量累积

46.     %进出油

47.     %当前时间处于出油时间, 即t∈[t 0+100·(i-

1),2.4+t 0+100·(i-1)]

48.       if t0/100- floor(t0/100) <= t/100-

 floor(t/100) && t/100- floor(t/100) < (t0+2.45)/100-

 floor((t0+2.45)/100)

49.          t1=100*((t-t0)/100- floor((t-t0)/100));%于

t 0+100·(i-1)时刻

50.          dQ=-0.85*S( round(100*t1)+1)* sqrt(2* pt/ rou)* dt;

51.            dm= dm+dQ* rou;

52.      end

53.     %进油

54.      if pc> pt

55.          dQ=0.85*0.7*0.7* pi* sqrt(2*( pc- pt)/ rouc)* dt;

56.             dm= dm+dQ* rouc;

57.          dQc=-0.85*0.7*0.7* pi* sqrt(2*( pc- pt)/ rouc)* dt;

58.            dmc= dmc+dQc* rouc;

59.      end

60.     %凸轮转动

61.      seta=w*t/1000;

62.      while seta>=2* pi          %限制在0到2* pi

63.           seta= seta-2* pi;

64.      end

65.      if set a <= pi              %升程小孔面积由分段拟合得到

66.           dr=fun1( seta)- rc;

67.           rc= funl( seta);

68.      else

69.           dr=fun2( seta)- rc;

 

70.           rc=fun2( seta);

71.       end

72.      %%核算高压油管状态

73.      m=m+ dm;

74.       rou=m/v;

75.       if rou<=0.8025

76.           pt=0.016182;

77.       else

78.           pt= Result(1, floor(10^5* rou-80229));

79.       end

80.     %%核算栓塞腔状态

81.      dv= pi*2.5*2.5* dr;

82.      vc= vc- dv;

83.      mc= mc+ dmc;

84.       rouc= mc/ vc;

85.     %判断油泵是否回油

86.      if rouc <= 0.80255               %(即对应压强小于0.5Mpa)

87.           pc=0.495479034926858;

88.          rouc=0.80255;

89.           mc= rouc* vc;

90.       else                          %核算压强

91.           pc= Result(1, floor(10^5* rouc-80229));

92.       end

93. end

94. clear dm dmc dQ dQc dr dt dv dvc i m mc pc pt rc rad Res ult rou rouc S set a t t0 t1 Talls v v0 vc w;

95. %% 结果作图

96. A=0:0.01*0.001: Tall/1000-0.01*0.001;

97. plot(A,p);

98. toc

99. clear A Tall;

100.     function [r]= funl( setat)

101.     r=0.00003* setat^6 + 0.018* setat^5 - 0.1429* setat^4 + 0.0504* setat^3 + 1.178* setat^2 + 0.0028* setat + 2.4126;

102.     end

103.     function [r]=fun2( setat)

104.       r=-1E-

05* setat^6 - 0.018* setat^5 + 0.4272* setat^4 - 3.6351* seta t^3 + 13.38* setat^2 - 20.644* setat + 16.659;

105.     end