高壓油管壓力控制問(wèn)題的微分模型與求解

徐興波

（淮陰工學(xué)院 數(shù)理學(xué)院，江蘇 淮安）

噴射系統(tǒng)是燃油發(fā)動(dòng)機(jī)的重要部件，如何保持高壓油管內(nèi)的壓力穩(wěn)定是提高燃油效率的關(guān)鍵。在噴射系統(tǒng)中，霧化低壓燃油經(jīng)由凸輪控制的高壓油泵壓縮為高壓燃油，再在壓力差的作用下通過(guò)單向閥涌入高壓油管，而高壓油管有若干個(gè)由針閥控制的錐型噴油器來(lái)進(jìn)行噴油，這一動(dòng)態(tài)過(guò)程被描述為高壓油管壓力控制問(wèn)題，即2019 年高教社杯全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽A 題，見(jiàn)圖1。建模競(jìng)賽題有三個(gè)子問(wèn)題，已有不少文章就該賽題進(jìn)行了研究[1-2]，但該問(wèn)題仍然值得進(jìn)一步細(xì)化討論。

圖1 具有減壓閥和兩個(gè)噴油嘴時(shí)高壓油管示意（源自賽題）

需要厘清三個(gè)方面的函數(shù)關(guān)系：(1) 彈性模量E(P)和壓力P 之間的關(guān)系，進(jìn)而精確給出燃油密度ρ(P)和壓力P 的關(guān)系；(2) 針閥運(yùn)動(dòng)hz(t)與時(shí)間t 的關(guān)系，進(jìn)而給出噴油截面積Sp(t)與時(shí)間t 的關(guān)系；(3) 凸輪極徑r(θ ) 與極角 θ之間的關(guān)系，進(jìn)而描述柱塞上下運(yùn)動(dòng)速度與時(shí)間的關(guān)系，也能導(dǎo)出柱塞腔的體積Vh及其變化率。對(duì)E(P)的擬合，陳興志等(2019)[3]用指數(shù)函數(shù)，何雅寧等(2019)[4]、宋常修等(2021)[5]、王一等(2021)[6]、趙澤宇等(2020)[7]均用多項(xiàng)式函數(shù)。楊曉琪等(2021)[8]、毛睿昕等(2022)[9]用數(shù)值定積分計(jì)算燃油密度，侯超鈞等(2021)[10]用了線性插值的方法計(jì)算E(P)。本文用二次多項(xiàng)式擬合E-1(P)，擬合精度好于之前結(jié)果，而且可以給出燃油密度的計(jì)算公式。何雅寧等(2019)對(duì)針閥升降程運(yùn)動(dòng)曲線采用高斯函數(shù)擬合，本文稍有改進(jìn)。蔡志杰(2020)[11]給出了賽題評(píng)價(jià)，但連續(xù)模型忽略了補(bǔ)充的低壓燃油。本文完善了連續(xù)模型并基于python 編程給出了合理的最優(yōu)解。

在模型建立前，做如下假設(shè)：(1) 燃油在高壓油管中近似為一維無(wú)粘性可壓縮理想流體，對(duì)高壓油管的形狀不產(chǎn)生影響；(2) 不考慮能量損失，燃油質(zhì)量守恒；(3) 噴油器將燃油噴向氣缸，氣缸壓力不低于標(biāo)準(zhǔn)大氣壓。

1 函數(shù)擬合

1.1 燃油密度

數(shù)值發(fā)現(xiàn)，即便用關(guān)于P 的三次多項(xiàng)式或分段函數(shù)來(lái)擬合彈性模量E(P)，所得結(jié)果也不十分理想。根據(jù)經(jīng)驗(yàn)公式，燃油的壓力變化量與密度變化量滿足，得到關(guān)于E(P)的定積分，這啟發(fā)用多項(xiàng)式來(lái)擬合E-1(P)，并且可以得到實(shí)用公式：

其中P=100 MPa， ρ= ρ (P) =0.850mg/mm3，且。對(duì)E(P)的二次多項(xiàng)式擬合可以采用基于最小二乘法的思想，在基于python 的具體實(shí)現(xiàn)上，可以scipy 庫(kù)中的optimize.leastsq 命令，擬合參數(shù)a≈ 6.492165×10-4，a≈-2.005 134×10-6，a≈1.180684×10-9，擬合值與原數(shù)據(jù)誤差均值為10-20量級(jí)、方差為10-14量級(jí)，擬合效果可以從圖2 中看到。另外，在幾種壓力值下的燃油密度是需要求出的，計(jì)算結(jié)果為 ρ (P= 0.5)≈ 0.804536，ρ(P= 150)≈0.867920，ρ (P= 160)≈0.871110。

圖2 彈性模量E(P)與壓力P 的函數(shù)關(guān)系擬合效果

1.2 凸輪運(yùn)動(dòng)

可以利用附件數(shù)據(jù)繪出凸輪極徑r 和極角θ 的函數(shù)圖像，另外根據(jù)x=rcosθ 和y=rsinθ可以繪出凸輪的形狀，可以看出凸輪的輪廓是光滑且軸對(duì)稱的，同時(shí)可以計(jì)算出極徑最大值rmax和最小值rmin之差為4.826 mm。利用python 命令optimize.leastsq 來(lái)對(duì)r(θ )進(jìn)行6 次多項(xiàng)式擬合，擬合誤差的均值為10-11量級(jí)而方差為10-6量級(jí)，盡管這樣的擬合精度已經(jīng)不錯(cuò)，但仍不夠理想，用更高階多項(xiàng)式又顯得復(fù)雜。凸輪曲線類似三角函數(shù)余弦曲線，利用函數(shù)

擬合可得到更好擬合精度，擬合誤差均值為10-16量級(jí)且方差為10-10量級(jí)。公式(3)既精確又簡(jiǎn)潔。凸輪曲線和擬合曲線可見(jiàn)圖3。

圖3 凸輪曲線和擬合曲線

凸輪上邊緣頂著柱塞。以凸輪中心為坐標(biāo)原點(diǎn)，豎直向上方向?yàn)檎较?，建立一維坐標(biāo)系。令θ (t) =ωt且 θ (0) =0，凸輪的最長(zhǎng)極徑在初始時(shí)刻豎直向上。凸輪的轉(zhuǎn)動(dòng)使得豎直方向的極徑周期地變化，進(jìn)而影響著柱塞的周期運(yùn)動(dòng)。柱塞的運(yùn)動(dòng)幅度即上下止點(diǎn)的距離為4.826 mm,上止點(diǎn)的速度為，這也導(dǎo)致柱塞腔體積Vh有所變化。關(guān)于Vh的計(jì)算，將在1.4 節(jié)給出。

1.3 針閥運(yùn)動(dòng)

針閥運(yùn)動(dòng)函數(shù)Z(t,t0)可以寫(xiě)成分段函數(shù)，升程和降程函數(shù)雖然可以利用關(guān)于時(shí)間的多項(xiàng)式來(lái)擬合，但利用高斯形式函數(shù)來(lái)擬合更有效。高斯函數(shù)在概率論中是用來(lái)描述正態(tài)分布密度函數(shù)的，本文采用如下形式的函數(shù)來(lái)擬合針閥升降程運(yùn)動(dòng)：

其中t0≥0表示針閥開(kāi)始升程運(yùn)動(dòng)起始時(shí)間，也是開(kāi)始噴油的時(shí)刻。賽題所提供的數(shù)據(jù)對(duì)應(yīng)于t0=0，利用python 命令optimize.leastsq 擬合升程階段函數(shù)z1，得args 為(2.02894068,0.45660788,0.02822825,-0.0067 0365)；擬合降程階段函數(shù)z2，得atgs 為(2.02932375,1.99283185,0.02823941,-0.00670779)。這兩個(gè)擬合函數(shù)的誤差均值、方差分別為10-10和10-5量級(jí)。當(dāng)然，令c2分別為0.45 和2.0, 然后再做最小二乘擬合也是可以接受的。

在噴油時(shí)間2.45 ms 范圍內(nèi)的針閥運(yùn)動(dòng)數(shù)據(jù)繪圖及其擬合曲線在圖4 中可見(jiàn)。針閥運(yùn)動(dòng)函數(shù)Z(t,t0)可以按五種情形來(lái)分段定義。(1) 在開(kāi)始噴油之前，有t＜t0，得Z(t,t0) ≡ 0。在噴油時(shí)間段內(nèi)，有t≥t0，令= mod(t-t0,Tz)和args=(c1,c2,c3)，其中Tz=100 ms 為針閥運(yùn)動(dòng)周期;（2） 如果? [0,0.45)，Z(t,t0)=z1(,args)；（3）如果 ～t?[0.45,2.0)，Z(t,t0) ≡2；（4） 如果?[2.0,2.45)，Z(t,t0) =z2(,args)；（5）如果?[2.45,100.0)，則Z(t,t0) ≡0。

圖4 針閥運(yùn)動(dòng)數(shù)據(jù)繪圖及其擬合曲線

1.4 燃油流量

燃油注入和噴出的流量實(shí)際是單位時(shí)間(ms)內(nèi)進(jìn)出高壓油管的燃油體積，賽題提供了經(jīng)驗(yàn)公式，也可參考蔡志杰(2020)。從單向閥注入高壓油管的高壓燃油壓力為Ph、密度為 ρh，而高壓油管內(nèi)的燃油壓力為Pg、密度為 ρg。記單向閥入口處的截面積為Sa= 0.49π。那么，當(dāng)單向閥開(kāi)啟時(shí)注入高壓油管的燃油流量為

當(dāng)單向閥關(guān)閉時(shí)，Qin≡0。

2 微分模型

2.1 子問(wèn)題一

質(zhì)量等于密度乘以體積，那么單位時(shí)間內(nèi)高壓油管內(nèi)燃油質(zhì)量的變化等于進(jìn)入和噴出燃油的流量差，又由于高壓油管體積恒定，便可以得到密度對(duì)應(yīng)的微分方程，再根據(jù)壓力與密度之間的微分方程關(guān)系，最后得到關(guān)于高壓油管內(nèi)燃油密度和壓力的一階常微分方程組：

初值為 ρ0= 0.850，P0=100。第一個(gè)問(wèn)題是如何選取未知參數(shù)Ta和t0使得壓力盡可能穩(wěn)定在100 MPa。為此，需要對(duì)方程組(5)進(jìn)行數(shù)值積分，得到節(jié)點(diǎn)數(shù)為N 的解向量，目標(biāo)函數(shù)為

第二個(gè)問(wèn)題是，如何控制Ta和t0，使得經(jīng)過(guò)T 毫秒時(shí)間，壓力盡可能穩(wěn)定在150 MPa?？梢允褂脙呻A段的調(diào)整策略。首先，調(diào)整參數(shù)，對(duì)方程組(5)數(shù)值積分T 毫秒，使得

其次，使用針對(duì)第一問(wèn)的方法，在目標(biāo)函數(shù)(6)中令P0=150，采用數(shù)值優(yōu)化方法求出參數(shù)，使得壓力盡量維持在150 MPa。

2.2 子問(wèn)題二

初值滿足 ρ= ρ= 0.850， = =100。該方程組含有 ω和t0兩個(gè)未知參數(shù)，問(wèn)題是確定這兩個(gè)參數(shù)使得高壓油管內(nèi)的壓力穩(wěn)定在100 MPa 左右。對(duì)方程組進(jìn)行數(shù)值積分，時(shí)間為T，節(jié)點(diǎn)數(shù)為N，在選取目標(biāo)函數(shù)時(shí)，用Pg[i] 替換式(6)中的Pi。

2.3 子問(wèn)題三

3 最優(yōu)解

3.1 子問(wèn)題一的解

起初，本文嘗試在固定t0的前提下用二分法來(lái)求Ta的最優(yōu)值，接著固定Ta，再利用二分法來(lái)求t0的最優(yōu)值，可求出近似值，代碼執(zhí)行效率雖不高，但能給出最優(yōu)值的估計(jì)。后來(lái)利用邊值問(wèn)題打靶法的思路，在參考文獻(xiàn)[12]中尋找優(yōu)化代碼，經(jīng)過(guò)驗(yàn)證，方法奏效。關(guān)鍵python 代碼如下：

x_cons_opt=optimize.minimize(f,[0.287,40.28],met hod='nelder-mead',

options={'maxiter':30,'xtol':0.01,'ftol':0.01,'disp': True},bounds=[bnd_ta,bnd_tj]).x

其中f 為目標(biāo)函數(shù)，通過(guò)對(duì)ODE 初值問(wèn)題(5)做數(shù)值積分，并計(jì)算目標(biāo)函數(shù)值；接著給出兩個(gè)參數(shù)的初值；所采用的優(yōu)化方法是對(duì)非線性多元函數(shù)求局部極小值有效的Nelder-Mead 方法；options 中的內(nèi)容很關(guān)鍵，它決定了程序何時(shí)停止，如果缺失options，程序運(yùn)行可能不收斂，maxiter 乘以參數(shù)的個(gè)數(shù)給出了最大迭代次數(shù)，xtol 和ftol 控制精度，bounds 給出了解的搜索范圍，比如T? [0.25,0.3]，t0?[20,60]。數(shù)值積分采用了scipy 庫(kù)函數(shù)，變步長(zhǎng)四階算法dopri5，選擇步長(zhǎng)為0.01，這保證了數(shù)值積分的精度。利用積分時(shí)長(zhǎng)2000 ms，迭代運(yùn)算53 次，第一小問(wèn)的最優(yōu)解為T= 0.28799396ms，t= 42.5590721ms，目標(biāo)函數(shù)(6)對(duì)應(yīng)的最優(yōu)值約2.6312 MPa，圖5 給出了對(duì)計(jì)算結(jié)果可靠的驗(yàn)證。

圖5 初始?jí)毫?00 MPa 時(shí)選擇Ta≈0.288 和t0≈42.56 使得高壓油管內(nèi)的壓力盡量穩(wěn)定在100 MPa

在計(jì)算初值的時(shí)候，也可以分析估計(jì)。在一個(gè)噴油周期內(nèi)對(duì)噴油流量進(jìn)行積分得44 毫升每百毫秒，進(jìn)油流量15.35 毫升每毫秒，假設(shè)在100 毫秒內(nèi)供噴平衡，得44 ≈ 15.35Ta×1 00/(10 +Ta)，而噴油次數(shù)為整數(shù)，可以估計(jì)出Ta?(0.2866,0.2955)，在不考慮噴油時(shí)間延遲的情況下Ta取0.2951 較為合理，但在短期來(lái)說(shuō)高壓油管壓力波動(dòng)太大并不合乎題意。

在上述基礎(chǔ)上，考慮壓力穩(wěn)定在150 MPa 的情形，式(5)的初值 (0.867920,150)。在分析方面，仍考慮100 毫秒內(nèi)的供噴平衡方程44 ≈ 6.27T×1 00/(10 +T)，可以估計(jì)出在初始?jí)毫?50 MPa、噴油不延時(shí)、長(zhǎng)期工作條件情況下，單向閥最佳開(kāi)啟時(shí)長(zhǎng)約0.7547 ms。在數(shù)值方面，先固定t0在40 附近，利用兩步二分法，也可以得到較好的結(jié)果，但所得結(jié)果T≈0.7537 、t0≈ 39.673，對(duì)比使用優(yōu)化算法所得結(jié)果Ta=0.753 84499、t0=45.7308504 還是稍遜一籌，后者為最優(yōu)解。在前者作為初值、積分時(shí)長(zhǎng)2000 ms、搜索范圍為[0.7,0.8]× [30,50]時(shí)，迭代次數(shù)23 次，最優(yōu)值約為3.1801，圖6 可作為對(duì)最優(yōu)解的驗(yàn)證。

圖6 初始?jí)毫?50 MPa 時(shí)選擇Ta≈0.754 和t0≈45.73 使得高壓油管內(nèi)的壓力盡量穩(wěn)定在150 MPa

考慮2 秒后實(shí)現(xiàn)將壓力P 從100 增加到150，需要調(diào)整參數(shù)Ta和t0，使得式(7)的值最接近0，利用二分法可以給出較好的估計(jì)值T≈ 0.73125、t≈1004.88，目標(biāo)函數(shù)值約7.2×10-3。在該參數(shù)值附近進(jìn)行優(yōu)化，更優(yōu)的目標(biāo)函數(shù)值約為1.32×10-4，最優(yōu)解為Ta=0.730357993、t0=1003.04631。在2 秒后采用壓力穩(wěn)定在150 MPa 的參數(shù)，數(shù)值圖像見(jiàn)圖7。同理，考慮5 秒后P=150 情形，得最優(yōu)解Ta=0.70251709，t0≈3.0×10-5≈0，目標(biāo)函數(shù)值式(7)約5.03×10-4；考慮10 秒后P=150 情形，得最優(yōu)解T≈ 0.700625，t0≈0，目標(biāo)函數(shù)值式(7)約0.150。

圖7 兩階段控制策略使得高壓油管內(nèi)的壓力在2 秒后盡量穩(wěn)定在150 MPa

3.2 子問(wèn)題二的解

問(wèn)題二是個(gè)強(qiáng)非線性問(wèn)題?？梢苑e分求出噴油器在2.45 毫秒也就是一個(gè)噴油周期100 毫秒時(shí)間內(nèi)的噴油量約38.7086，而在凸輪轉(zhuǎn)動(dòng)一周的時(shí)間內(nèi)供油量約為67.706，供噴平衡方程近似為67.706×100ω /(2π )≈38.7086，得ω≈0.0359rad/ms，可以估計(jì)出凸輪角速度的數(shù)量級(jí)。數(shù)值計(jì)算得最優(yōu)解是ω=0.0271478506、t0=109.194263，圖8 給出了高壓油管壓力控制在(97.475,102.55)的效果。

圖8 通過(guò)凸輪轉(zhuǎn)速和噴油時(shí)間延遲控制高壓油管壓力在100 MPa 附近效果

3.3 子問(wèn)題三的解

假設(shè)甲乙兩個(gè)噴油器噴油時(shí)間間隔t1，待優(yōu)化的參數(shù)有三個(gè)，即凸輪轉(zhuǎn)ω速、甲的噴射時(shí)間t0、乙的相對(duì)甲延遲噴射的時(shí)間t1。問(wèn)題的難度相對(duì)子問(wèn)題二增加了優(yōu)化計(jì)算量。本文給出較優(yōu)的解，即 ω=0.054 0722，t0=74.1114，t1=44.3691，目標(biāo)函數(shù)值 (6) 的值約6.3647。高壓油管內(nèi)的壓力波動(dòng)見(jiàn)圖9。

圖9 雙噴油嘴情況下進(jìn)行參數(shù)控制使得高壓油管內(nèi)的壓力在100 MPa 附近波動(dòng)

考慮再增加一個(gè)減壓?jiǎn)蜗蜷y。如果假設(shè)減壓閥每工作t3時(shí)間便關(guān)閉時(shí)間t4，這樣并不能求出最優(yōu)解。

而假設(shè)減壓閥在壓力大于102 時(shí)打開(kāi)，否則關(guān)閉，那么仍需優(yōu)化3 個(gè)參數(shù)。經(jīng)仿真，得到最優(yōu)解5.67 ×10 ,75.1114, 45.25，目標(biāo)函數(shù)值 (6) 的值約4.56，從圖10 可以看出優(yōu)化效果明顯。

圖10 有減壓閥的雙噴油嘴情況下高壓油管壓力控制效果

4 結(jié)論

本文以2019 年大學(xué)生數(shù)學(xué)建模題為問(wèn)題來(lái)源，完整地給出了高壓油管壓力控制問(wèn)題的微分方程模型，基于python 優(yōu)化命令做數(shù)值計(jì)算，得到了關(guān)于燃油密度、針閥運(yùn)動(dòng)、凸輪曲線的精確擬合函數(shù)，最終給出了可靠的最優(yōu)參數(shù)。由于噴油嘴外界的壓力未知，本文采用標(biāo)準(zhǔn)大氣壓。所用方法和所得結(jié)論為燃油噴射系統(tǒng)精準(zhǔn)控制提供有意義的參考。