陳益宏， 徐剛強(qiáng)， 李孝祿， 陳 源， 李運(yùn)堂

(1.中國(guó)計(jì)量大學(xué)機(jī)電工程學(xué)院， 浙江杭州 310018;2.浙江省智能制造質(zhì)量大數(shù)據(jù)溯源與應(yīng)用重點(diǎn)實(shí)驗(yàn)室， 浙江杭州 310018;3.浙江博星工貿(mào)有限公司技術(shù)中心， 浙江金華 321016)

引言

液體傳動(dòng)系統(tǒng)廣泛應(yīng)用于機(jī)械裝備中[1-2]。由于傳動(dòng)液混入空氣、溫度過高等因素，傳動(dòng)系統(tǒng)內(nèi)經(jīng)常產(chǎn)生復(fù)雜的氣液兩相流流動(dòng)[3-4]。當(dāng)傳動(dòng)壓力引起傳動(dòng)介質(zhì)壓力低于空氣分離壓時(shí)，溶解在傳動(dòng)液中的空氣析出，產(chǎn)生氣泡。同時(shí)，傳動(dòng)管內(nèi)產(chǎn)生氣泡也與系統(tǒng)溫度的變化有關(guān)。氣泡的存在影響傳動(dòng)系統(tǒng)傳遞動(dòng)力。另外，氣泡的破裂除了對(duì)部件造成損害外，還會(huì)產(chǎn)生噪聲。為此，建立傳動(dòng)管內(nèi)氣液兩相流含氣率模型，模擬空化過程，可以預(yù)測(cè)傳動(dòng)管內(nèi)的壓力和含氣率的變化情況，為精確控制液體傳動(dòng)系統(tǒng)的性能提供理論指導(dǎo)，并有助于防止沖蝕損傷的發(fā)生[5-6]。

國(guó)內(nèi)外學(xué)者建立了不同的管道流動(dòng)模型：CATANIA等[7-8]建立了正壓空化模型，采用二階精確數(shù)值格式進(jìn)行求解，在摩擦產(chǎn)生的熱量完全通過管壁傳遞的假設(shè)下，對(duì)壓力管內(nèi)聲空化瞬態(tài)過程進(jìn)行數(shù)值預(yù)測(cè)，結(jié)果表明，預(yù)測(cè)的壓力時(shí)程與實(shí)驗(yàn)的壓力時(shí)程幾乎一致。李汝寧等[9]基于氣泡模型建立了氣液兩相條件下的柴油機(jī)噴油系統(tǒng)模型，進(jìn)行數(shù)值仿真，對(duì)比了仿真結(jié)果和實(shí)驗(yàn)數(shù)據(jù)，結(jié)果表明：壓力大于45 MPa后，氣泡體積減小速度增大，造成壓力振蕩，仿真結(jié)果和實(shí)驗(yàn)數(shù)據(jù)吻合。唐東林等[10]依據(jù)含氣油液中氣相成分隨壓力變化過程，預(yù)測(cè)了有效體積彈性模量的變化，結(jié)果表明：壓力越小，有效體積彈性模量受初始含氣量影響越明顯。

另外，求解方法是提升數(shù)值解穩(wěn)定性和準(zhǔn)確性的關(guān)鍵。為此，水鴻壽[11]介紹了一維流體力學(xué)的差分方法，其精確度和穩(wěn)定性也各不相同。李松晶等[12-13]針對(duì)傳動(dòng)管中的數(shù)學(xué)模型，采用MATLAB Simulink中SELECTOR模塊對(duì)偏微分方程進(jìn)行空間上的積分，并與特征線法的仿真結(jié)果進(jìn)行對(duì)比，結(jié)果表明：傳動(dòng)閥突然關(guān)閉后管道中出現(xiàn)壓力脈動(dòng)變化與特征線法的仿真結(jié)果幾乎一致。方桂花[14]等通過一般迎風(fēng)差分方法對(duì)氣動(dòng)管路進(jìn)行分析后，再利用緊致差值方法對(duì)管路內(nèi)氣柱的狀態(tài)進(jìn)行參數(shù)分析，結(jié)果表明：利用緊致插值方法所分析的狀態(tài)值與差分格式理論值基本吻合，證明了此方法能應(yīng)用于氣體流體領(lǐng)域。

從上述研究可以看出，大多研究集中在傳動(dòng)介質(zhì)的壓力變化，較少研究含氣率變化帶來的影響。本研究分析了聲空化條件下傳動(dòng)液中空氣析出與溶解過程，并應(yīng)用于發(fā)動(dòng)機(jī)燃油噴射系統(tǒng)模型中。首先，在傳動(dòng)管中建立傳動(dòng)液質(zhì)量守恒方程和動(dòng)量守恒方程，得到空氣溶解和析出的氣泡模型，并考慮傳動(dòng)液壓縮帶來的溫度變化，編寫了MATLAB計(jì)算程序進(jìn)行求解，得到氣液兩相流主要參數(shù)的變化，并分析其對(duì)含氣率的影響。

1 數(shù)學(xué)模型

1.1 傳動(dòng)管道內(nèi)一維斜壓流流動(dòng)模型

連接柱塞式噴油泵和噴油器之間的傳動(dòng)管為一段長(zhǎng)度為L(zhǎng)=0.6 m的等截面圓管，如圖1所示。其中d為傳動(dòng)管直徑，取為2 mm，x為沿管道的軸向變量，dx為管內(nèi)取得的一小段微元。

圖1 傳動(dòng)管結(jié)構(gòu)圖Fig.1 Transmission tube structure diagram

在傳動(dòng)管中，可壓縮的純液體或氣液均勻混合物的一維斜壓流流動(dòng)通常采用質(zhì)量和動(dòng)量守恒方程進(jìn)行建模，利用熱力學(xué)能的變化建立能量守恒方程。模型建立在流體靜壓傳動(dòng)的基礎(chǔ)上，根據(jù)傳動(dòng)管內(nèi)壓縮前后傳動(dòng)液的質(zhì)量守恒以及流動(dòng)對(duì)管壁剪切應(yīng)力的考慮，得到了一維流體力學(xué)方程組：

(1)

式中，ρ—— 混合流體密度，kg/m3

u—— 傳動(dòng)液流速，m/s

p—— 傳動(dòng)管內(nèi)壓力，Pa

τw—— 管壁剪切應(yīng)力，Pa

t—— 時(shí)間變量，s

(2)

其中，

(3)

(4)

式中，vs—— 聲速，m/s

f—— 傳動(dòng)液流動(dòng)阻力系數(shù)

當(dāng)雷諾數(shù)Re≤2000時(shí)，流動(dòng)為層流，f=64/Re；當(dāng)雷諾數(shù)Re>2000時(shí)，流動(dòng)為紊流，f=0.316/Re0.25。

對(duì)于斜壓流模型，必須增加流體的狀態(tài)和轉(zhuǎn)換方程。本研究中，將流體看作是純液體和一定量的空氣的均質(zhì)混合物。液態(tài)和氣態(tài)成分的混合物在宏觀上被認(rèn)為是均質(zhì)和同性的系統(tǒng)，因此將空氣質(zhì)量分?jǐn)?shù)和體積分?jǐn)?shù)作為變量引入：

(5)

式中，μ—— 質(zhì)量分?jǐn)?shù)，%

g —— 代表空氣

m—— 混合物質(zhì)量，kg

V—— 混合物體積，m3

空氣密度ρg由理想氣體的狀態(tài)方程來定義：

(6)

式中，Rg—— 空氣氣體常數(shù)，J/(kg·K)

T—— 溫度，K

將比體積的定義用于混合物各組分的表達(dá)，并引入式(5)的質(zhì)量分?jǐn)?shù)，均勻混合物的密度可表示為壓力、溫度和未溶解空氣質(zhì)量分?jǐn)?shù)的函數(shù)：

(7)

式中， 下標(biāo)l —— 傳動(dòng)液

對(duì)于斜壓流流動(dòng)，通常在空化區(qū)，溫度T是壓力p的函數(shù)，T=T(p)，在方程式(7)中傳動(dòng)液密度是壓力p和溫度T的函數(shù)，ρl=ρl(T,p)。

對(duì)于可壓縮流體的一維流動(dòng)，能量守恒方程可以表示為：

(8)

式中，h—— 系統(tǒng)單位質(zhì)量的焓，J/kg

v—— 混合流體的比體積，m3/kg

Q—— 從邊界傳遞到系統(tǒng)的單位質(zhì)量的熱功率，J/(kg·s)

q—— 系統(tǒng)內(nèi)單位質(zhì)量的黏性功率耗散，

J/(kg·s)

將焓的變化用溫度的變化直觀表示，引入熵s以及熱力學(xué)第一定律：

Tds=dh-vdp

(9)

考慮純液體流動(dòng)時(shí)傳動(dòng)管內(nèi)的壓縮帶來溫度的變化，可以推導(dǎo)出以下方程：

(10)

式中,Cp—— 定壓比熱容，J/(kg·K)

vsT—— 等溫聲速，m/s

α—— 體積膨脹系數(shù)，1/K

ET—— 等溫體積彈性模量，Pa

在等熵條件下，式(10)將任何過程中液體的壓力和密度變化與溫度變化聯(lián)系起來，得到以下關(guān)系式：

(11)

1.2 氣泡模型

傳動(dòng)管內(nèi)的傳動(dòng)液中本身溶解一部分空氣，且以氣態(tài)存在的空氣通過吸附溶解，直到達(dá)到某種飽和平衡狀態(tài)。在一定的壓力和溫度下，傳動(dòng)液中可以溶解的空氣的溶解度為常數(shù)，且空氣在各種液體中的溶解度的經(jīng)驗(yàn)關(guān)系是已知的，而在高壓條件下溶解度隨壓力變化不大。由于本研究的選擇壓力范圍較大，通過對(duì)傳動(dòng)液不同溫度下的溶解度擬合溶解度系數(shù)隨溫度變化的函數(shù)：

S=4e-5T2-2.743e-2T+4.838

(12)

由亨利定律可知，在等溫條件下溶解在傳動(dòng)液中的初始?xì)馀蒹w積與初始?jí)毫Τ烧汝P(guān)系：

(13)

式中，Vc—— 初始溶解于傳動(dòng)液的氣泡體積，m3

pc—— 初始?jí)毫Γ琍a

p1—— 大氣壓，Pa

S—— 溶解度系數(shù)

V1—— 壓力為p1時(shí)溶解于傳動(dòng)液的氣泡體積，m3

當(dāng)傳動(dòng)管內(nèi)的壓力降低到空氣分離壓以下時(shí)，空氣從傳動(dòng)液中釋放出來，形成氣泡。當(dāng)壓力達(dá)到空氣分離壓以上，空氣重新溶解，氣體體積減小。若初始?xì)馀葸_(dá)到溶解與析出平衡狀態(tài)，則t時(shí)刻的傳動(dòng)液中的氣泡體積為：

Vt=Vc+V0-Vtr

(14)

式中，Vt——t時(shí)刻氣泡體積，m3

V0—— 初始?xì)馀蒹w積，m3

Vtr——t時(shí)刻溶解于傳動(dòng)液氣泡體積，m3

當(dāng)傳動(dòng)管內(nèi)處于壓力p時(shí)，空氣析出與溶解處于平衡狀態(tài)下，溶解的空氣體積為：

(15)

t時(shí)刻溶解于傳動(dòng)液氣泡體積為：

(16)

式中，XC—— 空氣析出時(shí)間常數(shù)，s

RJ—— 空氣溶解時(shí)間常數(shù)，s

pj—— 空氣分離壓，Pa

混合物的總密度ρ、體積彈性模量E和聲速vs的公式為：

(17)

式中，φ—— 傳動(dòng)液含氣率，%

El—— 傳動(dòng)液體積彈性模量，Pa

Eg—— 空氣體積彈性模量，Pa

1.3 求解方法

對(duì)式(2)采用特征線法：

(18)

(19)

圖2 x-t平面上的特征線Fig.2 Characteristic lines in the x-t plane

這樣就把瞬態(tài)偏微分式(2)轉(zhuǎn)化成2個(gè)獨(dú)立的常微分式(19)，簡(jiǎn)化了管道內(nèi)部流動(dòng)的求解過程。

采用數(shù)值求解的方法，將管道沿長(zhǎng)度方向離散成若干管段，將瞬變過程離散成若干個(gè)時(shí)間步長(zhǎng)，劃分情況如圖3所示。將管道分成i段，每段為dx，時(shí)間分為n段，每段為dt，以管長(zhǎng)為橫坐標(biāo)，以時(shí)長(zhǎng)為縱坐標(biāo)，組成x-t平面網(wǎng)格節(jié)點(diǎn)分布圖。

圖3 x-t平面網(wǎng)格節(jié)點(diǎn)分布Fig.3 Distribution of grid nodes in the x-t plane

將兩組特征方程進(jìn)行離散化，用三點(diǎn)差分法得到一階差分方程：

[p(i-1,n-1)+p(i+1,n-1)]}

(20)

(21)

結(jié)合式(8)～式(11)以及傳動(dòng)液的物性參數(shù)，對(duì)溫度進(jìn)行差分得到一階差分方程：

<1|p(i,n)-p(i,n-1)>1|

(22)

2 求解參數(shù)

采用ISO 4113試驗(yàn)油模擬傳動(dòng)液的體積彈性模量、密度和體積膨脹系數(shù) (用來替代柴油物理特性[7])。通過對(duì)該數(shù)據(jù)的比較分析，擬合出溫度T和壓力p相關(guān)的函數(shù)，如式(23)～式(25)所示，并將該函數(shù)引用到程序中使用：

El=(-4.712e-4T2+0.2935T-35.35)p-

9.177e6T+4.226e9

(23)

ρl=(1.737e-11T2-9.149e-9T+1.65e-6)p-

0.6675T+1020

(24)

α=(6e-23T+2.36e-20)p2-p·(8.967e-15T+

5.411e-12)+6.85e-7T+8.098e-4

(25)

偏微分方程組式(2)的初值條件為：

(26)

壓力邊界條件考慮到模擬傳動(dòng)液壓縮和膨脹的過程，即從進(jìn)口通入的壓力波和出口處的壓力波有一定的壓差。采用文獻(xiàn)[15]中的壓力邊界條件，如圖4所示。

圖4 壓力邊界條件[15]Fig.4 Pressure boundary condition[15]

相應(yīng)地，模擬所用其他參數(shù)如表1所示。

表1 主要參數(shù)Tab.1 Main parameters

3 數(shù)值結(jié)果與分析

3.1 傳動(dòng)管內(nèi)壓力變化及對(duì)含氣率的影響

圖5為傳動(dòng)管內(nèi)的壓力變化。從圖中可知，傳動(dòng)介質(zhì)初始?jí)毫?.09 MPa，最高達(dá)到50 MPa左右。在管道內(nèi)壓力是波動(dòng)的，仿真結(jié)果很好地模擬了壓力波動(dòng)。由于傳動(dòng)管路中常摻入空氣，使壓力波傳播速度減慢，導(dǎo)致系統(tǒng)壓力不穩(wěn)定。這與文獻(xiàn)[16]中的壓力波動(dòng)情況類似。

圖5 傳動(dòng)管內(nèi)壓力變化Fig.5 Pressure variation in transmission tube

圖6為傳動(dòng)管內(nèi)不同位置的含氣率、壓力變化。結(jié)果表明，傳動(dòng)管進(jìn)口附近的含氣率要高于出口附近，且壓力峰值出現(xiàn)延時(shí)現(xiàn)象。這是由于壓力波在可壓縮性驅(qū)動(dòng)的聲速下向進(jìn)口傳播，傳動(dòng)管進(jìn)口附近的壓力上升較慢，達(dá)到峰值較慢，但壓力下降較快，迅速下降至空氣分離壓以下，而傳動(dòng)管出口附近的壓力甚至未降至空氣分離壓。這與文獻(xiàn)[7,12]的延遲現(xiàn)象相一致。因此，在圖7中可以看到含氣率也相應(yīng)出現(xiàn)了延遲情況，而進(jìn)口附近的空化現(xiàn)象相比出口附近更加劇烈。

圖6 傳動(dòng)管內(nèi)不同位置的含氣率、壓力變化Fig.6 Variation of air volume fraction and pressure at different positions in transmission tube

圖7 傳動(dòng)管內(nèi)不同位置的含氣率變化Fig.7 Variation of air volume fraction at different positions in transmission tube

3.2 傳動(dòng)管內(nèi)流速變化

圖8為傳動(dòng)管內(nèi)流速變化。從圖中可以看出傳動(dòng)介質(zhì)的正向流速最高在50 m/s左右，逆向流速在40 m/s左右。流速的起伏是由于在變化過程中，壓差的不斷改變，影響傳動(dòng)液流速。圖9為不同壓差下流速隨時(shí)間的變化。結(jié)果表明，初始?jí)毫β孕∮诔隹趬毫?，開始階段傳動(dòng)液尚未流動(dòng)，隨著壓差的變化，產(chǎn)生了正向或逆向的流動(dòng)。這也證實(shí)了傳動(dòng)系統(tǒng)中傳動(dòng)液的壓縮和膨脹過程。

3.3 傳動(dòng)管內(nèi)含氣率變化

圖10為傳動(dòng)管內(nèi)含氣率變化。給定管內(nèi)初始含氣率為1%。在空化發(fā)生前，含氣率隨壓力升高而降低；空化發(fā)生后的含氣率高于1%。這是因?yàn)槿芙庥谝后w中的空氣多于液體中析出的空氣；空化后快速產(chǎn)生氣泡，使含氣率高于初始值。圖11所示為不同析出和溶解空氣時(shí)間常數(shù)下含氣率的變化。結(jié)果表明，含氣率隨空氣溶解時(shí)間常數(shù)增大而減小，但在空化區(qū)域，含氣率未發(fā)生改變；空氣析出時(shí)間常數(shù)小于1 s在空化區(qū)域的含氣率變化幅度要高于空氣析出時(shí)間常數(shù)大于1 s時(shí)的含氣率變化幅度。這是由于氣泡模型將空化區(qū)和非空化區(qū)分開，空化區(qū)的流體特性受到氣液相變化的影響，壓力波的傳播速度會(huì)明顯降低；含氣率主要受氣體析出時(shí)間常數(shù)影響，小于1 s的氣體析出時(shí)間常數(shù)時(shí)，模型對(duì)參數(shù)的變化敏感。這與文獻(xiàn)[17]中對(duì)氣體析出時(shí)間常數(shù)的研究結(jié)果相符。氣泡生成會(huì)抑制流體傳動(dòng)，為使傳動(dòng)液流動(dòng)時(shí)含氣率控制在較小范圍內(nèi)，降低空氣析出速率能夠有效減少氣泡。

圖10 傳動(dòng)管內(nèi)含氣率變化Fig.10 Variation of air volume fraction in transmission tube

圖11 析出、溶解時(shí)間常數(shù)對(duì)含氣率的影響Fig.11 Effect of precipitation and dissolution time constants on air volume fraction

3.4 傳動(dòng)管內(nèi)溫度變化及對(duì)含氣率的影響

圖12為傳動(dòng)管內(nèi)溫度變化。從圖中可以看出，溫度隨時(shí)間變化幅值在3 ℃以內(nèi)，正向流動(dòng)時(shí)溫度增速比逆向流動(dòng)時(shí)要快。這是由于壓力下降時(shí)間較短，對(duì)于管壁的摩擦較為劇烈，在此過程中管道內(nèi)的溫度迅速增高。文獻(xiàn)[8]中的壓力在100 MPa左右，溫度變化在7 ℃以內(nèi)，因此，3 ℃的溫度變化是合理的。圖13為溫度對(duì)含氣率的影響。結(jié)果表明，溫度在35 ℃以下時(shí)含氣率最小，60～80 ℃時(shí)的含氣率最大，80～90 ℃時(shí)比35～60 ℃時(shí)的含氣率要大。由于傳動(dòng)系統(tǒng)的工作溫度要高于室溫，但溫度過高會(huì)使部件性能變差?？紤]含氣率與運(yùn)行溫度的因素，宜將傳動(dòng)系統(tǒng)承受的溫度控制在一定范圍以內(nèi)。

圖12 傳動(dòng)管內(nèi)溫度變化Fig.12 Variation of temperature in transmission tube

圖13 溫度對(duì)含氣率的影響Fig.13 Effect of temperature on air volume fraction

3.5 含氣率對(duì)體積彈性模量的影響

圖14和圖15為傳動(dòng)管內(nèi)聲速和體積彈性模量的變化。從圖中可以看出，初始聲速在930 m/s左右，隨時(shí)間而增大，在0.12 s左右達(dá)到峰值，隨后減??；在管內(nèi)出口附近的聲速略高于進(jìn)口附近；體積彈性模量與聲速的變化趨勢(shì)相近。這是因?yàn)椴ㄊ菙_動(dòng)區(qū)和未擾動(dòng)區(qū)的分界面，穿過此界面，擾動(dòng)使傳動(dòng)介質(zhì)的壓強(qiáng)升高，后面的波的傳播速度比前面的波快，最終疊在一起形成一道突躍的壓縮波，使聲速增大；進(jìn)口附近的含氣量要高于出口附近，氣液相態(tài)變化影響聲波傳遞；初始含氣量使得初始聲速處于較低值，聲速峰值在1200 m/s左右，這與文獻(xiàn)[7]中的聲速相近。圖17為x=0.2 m處的不同含氣率下體積彈性模量和聲速隨時(shí)間的變化。結(jié)果表明，傳動(dòng)液流速始終低于聲速(即處于亞聲速的狀態(tài))；在圖16中，隨著傳動(dòng)管內(nèi)含氣量的增加，體積彈性模量明顯減小，而空化區(qū)的含氣率增幅不大，對(duì)體積彈性模量影響較小。這是由于稀疏波對(duì)管內(nèi)壓力傳遞產(chǎn)生了影響，降低了傳動(dòng)液的抗壓縮能力，導(dǎo)致傳動(dòng)性能變差，這與文獻(xiàn)[10]中的體積彈性模量變化趨勢(shì)一致。

圖14 傳動(dòng)管內(nèi)聲速變化Fig.14 Variation of sound velocity in transmission tube

圖15 傳動(dòng)管內(nèi)體積彈性模量變化Fig.15 Variation of bulk elastic modulus in transmission tube

圖16 含氣率對(duì)體積彈性模量的影響Fig.16 Effect of air volume fraction on bulk elastic modulus

圖17 x=0.2 m處的不同含氣率下體積彈性模量和聲速隨時(shí)間的變化Fig.17 Variation of bulk elastic modulus and sound velocity with time at x=0.2 m with different air volume fraction

4 結(jié)論

(1) 在本研究給定的條件下，傳動(dòng)管道內(nèi)的傳動(dòng)液壓縮和膨脹過程帶來的溫度變化在3 ℃以內(nèi)，且空化區(qū)域的溫度變化較?。?/p>

(2) 在非空化區(qū)，壓力越大，含氣率越?。辉诳栈瘏^(qū)，含氣率短時(shí)間內(nèi)會(huì)迅速增加，隨后緩慢增加；傳動(dòng)管內(nèi)壓力越大，空化區(qū)域越少，含氣率較低；

(3) 體積彈性模量隨含氣率的增大而減少；當(dāng)發(fā)生空化現(xiàn)象時(shí)，含氣率變化不大，對(duì)體積彈性模量的影響較小。