高壓油管壓力控制問題的微分模型與求解

徐興波

（淮陰工學院 數理學院，江蘇 淮安）

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

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

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

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

1 函數擬合

1.1 燃油密度

數值發現，即便用關于P 的三次多項式或分段函數來擬合彈性模量E(P)，所得結果也不十分理想。根據經驗公式，燃油的壓力變化量與密度變化量滿足，得到關于E(P)的定積分，這啟發用多項式來擬合E-1(P)，并且可以得到實用公式：

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

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

1.2 凸輪運動

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

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

圖3 凸輪曲線和擬合曲線

凸輪上邊緣頂著柱塞。以凸輪中心為坐標原點，豎直向上方向為正方向，建立一維坐標系。令θ (t) =ωt且 θ (0) =0，凸輪的最長極徑在初始時刻豎直向上。凸輪的轉動使得豎直方向的極徑周期地變化，進而影響著柱塞的周期運動。柱塞的運動幅度即上下止點的距離為4.826 mm,上止點的速度為，這也導致柱塞腔體積Vh有所變化。關于Vh的計算，將在1.4 節給出。

1.3 針閥運動

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

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

在噴油時間2.45 ms 范圍內的針閥運動數據繪圖及其擬合曲線在圖4 中可見。針閥運動函數Z(t,t0)可以按五種情形來分段定義。(1) 在開始噴油之前，有t＜t0，得Z(t,t0) ≡ 0。在噴油時間段內，有t≥t0，令= mod(t-t0,Tz)和args=(c1,c2,c3)，其中Tz=100 ms 為針閥運動周期;（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 針閥運動數據繪圖及其擬合曲線

1.4 燃油流量

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

當單向閥關閉時，Qin≡0。

2 微分模型

2.1 子問題一

質量等于密度乘以體積，那么單位時間內高壓油管內燃油質量的變化等于進入和噴出燃油的流量差，又由于高壓油管體積恒定，便可以得到密度對應的微分方程，再根據壓力與密度之間的微分方程關系，最后得到關于高壓油管內燃油密度和壓力的一階常微分方程組：

初值為 ρ0= 0.850，P0=100。第一個問題是如何選取未知參數Ta和t0使得壓力盡可能穩定在100 MPa。為此，需要對方程組(5)進行數值積分，得到節點數為N 的解向量，目標函數為

第二個問題是，如何控制Ta和t0，使得經過T 毫秒時間，壓力盡可能穩定在150 MPa?？梢允褂脙呻A段的調整策略。首先，調整參數，對方程組(5)數值積分T 毫秒，使得

其次，使用針對第一問的方法，在目標函數(6)中令P0=150，采用數值優化方法求出參數，使得壓力盡量維持在150 MPa。

2.2 子問題二

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

2.3 子問題三

3 最優解

3.1 子問題一的解

起初，本文嘗試在固定t0的前提下用二分法來求Ta的最優值，接著固定Ta，再利用二分法來求t0的最優值，可求出近似值，代碼執行效率雖不高，但能給出最優值的估計。后來利用邊值問題打靶法的思路，在參考文獻[12]中尋找優化代碼，經過驗證，方法奏效。關鍵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 為目標函數，通過對ODE 初值問題(5)做數值積分，并計算目標函數值；接著給出兩個參數的初值；所采用的優化方法是對非線性多元函數求局部極小值有效的Nelder-Mead 方法；options 中的內容很關鍵，它決定了程序何時停止，如果缺失options，程序運行可能不收斂，maxiter 乘以參數的個數給出了最大迭代次數，xtol 和ftol 控制精度，bounds 給出了解的搜索范圍，比如T? [0.25,0.3]，t0?[20,60]。數值積分采用了scipy 庫函數，變步長四階算法dopri5，選擇步長為0.01，這保證了數值積分的精度。利用積分時長2000 ms，迭代運算53 次，第一小問的最優解為T= 0.28799396ms，t= 42.5590721ms，目標函數(6)對應的最優值約2.6312 MPa，圖5 給出了對計算結果可靠的驗證。

圖5 初始壓力為100 MPa 時選擇Ta≈0.288 和t0≈42.56 使得高壓油管內的壓力盡量穩定在100 MPa

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

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

圖6 初始壓力為150 MPa 時選擇Ta≈0.754 和t0≈45.73 使得高壓油管內的壓力盡量穩定在150 MPa

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

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

3.2 子問題二的解

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

圖8 通過凸輪轉速和噴油時間延遲控制高壓油管壓力在100 MPa 附近效果

3.3 子問題三的解

假設甲乙兩個噴油器噴油時間間隔t1，待優化的參數有三個，即凸輪轉ω速、甲的噴射時間t0、乙的相對甲延遲噴射的時間t1。問題的難度相對子問題二增加了優化計算量。本文給出較優的解，即 ω=0.054 0722，t0=74.1114，t1=44.3691，目標函數值 (6) 的值約6.3647。高壓油管內的壓力波動見圖9。

圖9 雙噴油嘴情況下進行參數控制使得高壓油管內的壓力在100 MPa 附近波動

考慮再增加一個減壓單向閥。如果假設減壓閥每工作t3時間便關閉時間t4，這樣并不能求出最優解。

而假設減壓閥在壓力大于102 時打開，否則關閉，那么仍需優化3 個參數。經仿真，得到最優解5.67 ×10 ,75.1114, 45.25，目標函數值 (6) 的值約4.56，從圖10 可以看出優化效果明顯。

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

4 結論

本文以2019 年大學生數學建模題為問題來源，完整地給出了高壓油管壓力控制問題的微分方程模型，基于python 優化命令做數值計算，得到了關于燃油密度、針閥運動、凸輪曲線的精確擬合函數，最終給出了可靠的最優參數。由于噴油嘴外界的壓力未知，本文采用標準大氣壓。所用方法和所得結論為燃油噴射系統精準控制提供有意義的參考。