滑片泵進出口壓力對其內流特性影響的數(shù)值模擬

張群峰，閆盼盼，單建平，何萬發(fā)

(1. 北京交通大學 土木建筑工程學院，北京100044;2. 中航工業(yè)哈爾濱東安發(fā)動機集團有限公司，黑龍江 哈爾濱150066)

0 引言

空化［1］是潤滑油滑片泵工作過程中容易發(fā)生的一種物理現(xiàn)象?？栈l(fā)生時，滑片泵的工作性能受到很大的影響，往往導致滑片泵工作壽命縮短和效率下降，產生振動和噪聲?；脙鹊目栈饕獮槿~片表面空化和間隙空化［2］。葉片表面空化通常是由葉片背面的二次流引起的，間隙空化是潤滑油經(jīng)過葉片與定子壁面的徑向間隙以及葉片軸向端面間隙引起局部流速增加，壓力下降到潤滑油的當?shù)仫柡驼魵鈮憾纬傻摹Ｄ壳?，從公開發(fā)表的文獻來看，研究人員［3-10］通過理論分析和經(jīng)驗公式對滑片泵工作性能和空化現(xiàn)象做了一定的研究。由于滑片泵葉片的工作面積是變化的，僅靠理論分析和經(jīng)驗公式難以對滑片泵的空化特性作進一步的研究。隨著計算流體力學(CFD)的發(fā)展和計算機性能的提高，CFD 為進行滑片泵空化特性的研究提供了重要手段，已有研究人員對簡化滑片泵模型(簡化處理滑片頭部和定子之間的間隙形狀)進行了二維或三維的數(shù)值模擬［11-12］。

利用CFD 技術來模擬滑片泵的空化特性的關鍵問題是滑片泵在變容積過程中運動邊界的處理以及空化區(qū)氣液界面的處理。處理運動邊界的方法通常有基于貼體網(wǎng)格的任意拉格朗日歐拉(ALE)方法［13］、重疊網(wǎng)格法［14］和浸入邊界法［15］。由于浸入邊界法存在邊界上求解精度不高和存在體積誤差問題，目前沒有廣泛使用，而ALE 方法和重疊網(wǎng)格法則得到了較為廣泛的使用。空化區(qū)氣液界面的處理主要涉及空化模型和氣液界面的捕捉，目前在CFD軟件中廣泛采用的是流體體積(VOF)方法［16-18］和基 于Rayleigh-Plesset 方 程［1］的Rayleigh 空 化 模型［19-21］。本文通過自編程序與Star-cd 軟件相結合，采用ALE 方法準確模擬滑片運動過程中的相對位置，精確捕捉滑片頭部和定子之間的間隙形狀。并利用該功能，結合Menter SST k-ω 湍流模型［22-23］和Rayleigh 空化模型對某一航空發(fā)動機滑片泵的內部流場進行數(shù)值模擬，通過計算結果和實驗結果的比較確定合適的空化模型參數(shù)，分析滑片泵內的流動特性和空化產生、發(fā)展和潰滅的過程;研究滑片泵進、出口壓力對滑片泵工作性能和空化特性的影響。

1 控制方程與數(shù)值方法

1.1 控制方程和空化模型

潤滑油發(fā)生氣化時，氣液界面上的密度變化很大，而界面速度是連續(xù)的。為了保證數(shù)值計算的穩(wěn)定，通過求解三維N-S 方程并采用VOF 方法和Rayleigh 空化模型可以模擬空化的產生和潰滅等物理過程。

液態(tài)潤滑油和氣化潤滑油的連續(xù)方程分別為

式中:α 和β 分別為液態(tài)潤滑油和氣化潤滑油的體積分數(shù)，滿足α+β=1;Sβ為空化源項;ρl和ρv分別為液態(tài)潤滑油和氣化潤滑油密度;t 為時間;vgj為運動網(wǎng)格速度;uj為流動速度分量;xj為節(jié)點坐標值;j=1，2，3.

液態(tài)潤滑油和氣化潤滑油混合物的動量方程為

式中:ρ=ρlα+ρvβ;f 為液態(tài)潤滑油和氣化潤滑油界面動量輸運源項;為克羅內克算子，ui和uk為流動速度分量，xk為節(jié)點坐標值，i=1，2，3.

Sβ通過Rayleigh 空化模型來求解，該模型中，空化核的增長率利用Rayleigh-Plesset 方程來描述，該方程為

式中:R 為空化核半徑;pv為當?shù)仫柡驼魵鈮?σ 為表面張力系數(shù);p 為當?shù)仂o壓。令n 為單位液體體積內空化核數(shù)，對上述方程簡化后，可得

由以上各式可以看出，影響空化產生和空化區(qū)大小的因素主要是液體的飽和蒸氣壓、當?shù)仂o壓、空化核數(shù)和空化核的初始半徑。其中空化核數(shù)和氣泡初始半徑是需要人為給定的。研究表明，空化核的初始半徑對空化區(qū)域的大小影響較小，而空化核數(shù)目的大小對空化區(qū)域的大小影響較大，計算中需要通過實驗結果來確定。

模擬滑片泵的變容積工作過程時，網(wǎng)格發(fā)生變形或重構，為了避免控制方程在離散過程中產生人工質量源，應滿足空間守恒定律(SCL)［24］:

式中:nf為網(wǎng)格單元的表面數(shù)目;Sk為網(wǎng)格單元的第k 個面;Sj為表面Sk的面積矢分量。

1.2 湍流模型

式中:k 為湍流動能;μt為湍流粘性系數(shù)，采用Menter SST k-ω 湍流模型來計算湍流粘性系數(shù)。SST k-ω湍流模型在近壁處采用Wilcox k-ω 模型［25］，在邊界層邊緣和自由剪切層采用k-ε 模式，這兩個區(qū)域之間通過混合函數(shù)來過渡，因此SST k-ω 湍流模型兼有近壁面Wilcox k-ω 模型的穩(wěn)定性及邊界層外部k-ε模型高效性的優(yōu)點。具體形式為

1.3 計算方法

采用有限體積法在同位網(wǎng)格中對控制方程進行離散化，動量方程在空間上離散采用基于MUSCL格式［26-27］的二階精度差分格式。模擬空化和轉子轉動時，計算采用非穩(wěn)態(tài)計算，用PISO 方法［28］求解離散方程。引入空化模型后，計算容易發(fā)散，為了提高計算的收斂速度和穩(wěn)定性，計算初始時，不激活空化模型進行單相流計算，迭代到流場建立后，再激活空化模型進行計算。

2 網(wǎng)格生成與運動網(wǎng)格的實現(xiàn)

選用的滑片泵轉軸直徑為30.154 mm，葉片數(shù)為4 片，葉片轉動90°為一周期，因受安裝位置的限制，滑片泵的進口流道不是直接對著吸油腔，而是與吸油腔在轉動方向上近似呈90°夾角，圖1 為滑片泵的結構剖面圖。采用ICEM-CFD 軟件生成滑片泵進、出口流道、定子內流道的六面體計算網(wǎng)格，因滑片泵內部流道較為復雜，為了保證網(wǎng)格質量，分塊生成計算網(wǎng)格，各塊網(wǎng)格交接面非連續(xù)時，采用Couple 進行連接以保證各部分網(wǎng)格間的連通性，如圖2 所示。通過自編程序，根據(jù)定子輪廓線的曲線函數(shù)和滑片的運動規(guī)律，構建滑片泵轉動區(qū)域的內表面，采用代數(shù)網(wǎng)格生成法生成運動區(qū)域的網(wǎng)格并進行網(wǎng)格的光順處理，網(wǎng)格總數(shù)為57 萬左右。圖2 局部放大圖為端面間隙網(wǎng)格，這部分網(wǎng)格的分布隨著閥片的轉動而變化。圖3為轉動區(qū)域葉片附近的網(wǎng)格，自編網(wǎng)格生成程序可以實現(xiàn)葉片頭部與定子壁面間隙和端面間隙的網(wǎng)格生成，同時能保證間隙附近網(wǎng)格有較好的正交性。

圖1 結構剖面圖Fig.1 Cross-sectional view of vane pump

圖2 滑片泵計算網(wǎng)格Fig.2 Computational mesh of vane pump

計算時，根據(jù)物理時間確定滑片泵轉軸的位置，生成轉動區(qū)域的網(wǎng)格，并輸出網(wǎng)格節(jié)點號和坐標值，求解器利用用戶子程序功能，在進行當前時間步迭代之前讀入網(wǎng)格節(jié)點和坐標值，實現(xiàn)了當前時間步的網(wǎng)格重構，這種方法克服了傳統(tǒng)網(wǎng)格變形和重構時計算時間長和容易產生高畸變網(wǎng)格的缺點。

圖3 葉片附近的網(wǎng)格Fig.3 Mesh near vane

3 邊界條件和計算工況

實際運行中，滑片泵的進口和出口壓力是脈動的，加載時變的壓力值更為準確，但這些數(shù)據(jù)需要通過實驗測得，而且也難以將動態(tài)壓力與相應的閥片的轉動位置對應起來，因此計算時滑片泵的進口和出口設為定值壓力邊界條件?；玫乃斜诿娼o定無滑移、無滲透邊界條件，對于湍流問題，近壁區(qū)采用壁面函數(shù)。轉子和定子網(wǎng)格交接面為滑移面，設定Attached 邊界條件，保證轉子和定子網(wǎng)格的聯(lián)通?；玫念~定工況轉速為3 268 r/min，滑片泵的進口絕對壓力分別為47 kPa 、65 kPa、85 kPa 和101 kPa，滑片泵出口絕對壓力分別為301 kPa、501 kPa和1 001 kPa，溫度為140 ℃，動力粘度為0.002 36 Pa·s，密度為925.9 kg/m3. 考慮空化時，設油汽的密度為液態(tài)潤滑油的1/10 000、粘度為液態(tài)潤滑油的1/10，飽和蒸氣壓為100 Pa. 滑片泵葉片頭部和定子壁面間隙設為0.05 mm，軸向端面間隙設為0.02 mm.

4 空化模型參數(shù)的確定

為了確定空化模型中需要人為設定的單位液體體內空化核的數(shù)目，首先對某一在滑片泵專用實驗器上進行的實驗工況來進行模擬，該工況的溫度為140°，進口絕對壓力為47 kPa，出口絕對壓力為501 kPa，該工況的實驗測得的平均流量為38.1 L/min.因空化核的初始半徑對空化區(qū)域的大小影響較小，所以選取初始空化核半徑為10-6m，分別設單位液體體積的空化核個數(shù)為109、1010、1011和1012. 表1 分別為不同空化核數(shù)時，滑片泵進口和出口的總平均流量。

表1 空化核數(shù)對計算流量的影響Tab.1 Influence of nuclei number on calculated flow rate

根據(jù)比較，單位液體體積的空化核個數(shù)為1010時，流量的計算值和實驗值最接近，因此選擇該單位液體體積的空化核個數(shù)作為空化模型的設置參數(shù)。

5 計算結果分析

圖4 實驗工況下空化體積分數(shù)的分布Fig.4 Distribution of cavitation volume fraction on a section of vane pump under the experimental condition

5.1 滑片泵的流動特性

首先對上述實驗工況，且單位液體體積的空化核個數(shù)為1010的計算結果進行做進一步的分析。圖4 為不同轉角時，滑片泵某一截面上空化體積分數(shù)的分布，體積分數(shù)為1 表示空化后的潤滑油蒸氣，體積分數(shù)為0 表示液態(tài)潤滑油，0 與1 之間的值為液態(tài)潤滑和油潤滑油蒸氣的過渡區(qū)。圖5 為滑片泵進口和出口瞬時流量與轉子轉角的關系曲線，進口瞬時流量隨時間變化較小，而出口瞬時流量隨時間變化較大?？拷颓粋鹊倪M口流道內存在空化區(qū)，該區(qū)域的空化區(qū)域大小相對穩(wěn)定，使得進口流量隨轉角變化不明顯;在吸油腔內，葉片的背風側存在空化區(qū)，這些空化隨著葉片的轉動而不斷產生、發(fā)展和潰滅，并引起這些區(qū)域內質量的堆積(空化區(qū)減小)和釋放(空化區(qū)增大)，從而引起出口流量的脈動。當葉片轉動并打開封油區(qū)時，壓油腔內的壓力高于封油區(qū)內的壓力，潤滑油倒流回封油區(qū)內，此時，葉片背風區(qū)內的壓力升高，空化區(qū)減小，出口流量為負值，呈現(xiàn)吸油狀態(tài)。在周期時間內進口的平均流量為36.53 L/min，出口流量為37.66 L/min，二者的平均流量為37.10 L/min，進口和出口的平均流量基本一致，但由于引入空化模型后，數(shù)值上處理空化區(qū)域和液態(tài)潤滑油區(qū)域的界面時會存在一定的數(shù)值誤差，導致進、出口的平均流量有一定差別，但整體上計算是符合質量守恒定律的，計算值與實驗值非常接近。

圖5 滑油泵進、出口瞬時流量與轉子轉角的關系曲線Fig.5 Relation between inlet/outlet transient flow rates and rotor angle

圖6 轉動區(qū)域內空化體積分數(shù)的分布Fig.6 Distribution of cavitation volume fraction in rotating zone

圖6為某時刻的葉片頭部與定子壁面間隙和葉片軸向間隙的空化體積分數(shù)分布，與圖4 相比可知，葉片間隙引起的空化區(qū)域的大小和強度比葉片表面空化小。葉片端面間隙引起較大的泄露速度，使得對葉片端部背風側形成較大的空化區(qū)。

5.2 不同出口壓力對滑片泵流動特性的影響

為了研究滑片泵不同出口壓力對滑片泵流動特性的影響，在出口絕對壓力為501 kPa 工況的基礎上，增加了出口絕對壓力為301 kPa 和1 001 kPa 的工況。圖7 為不同出口壓力下，滑片泵出口瞬時流量與轉子轉角的關系曲線。由圖7 可知，出口壓力增加，封油區(qū)打開后，從壓油腔倒流回封油區(qū)的潤滑油增加，使得出油口的回流流量增加，但與回流區(qū)相對應的轉角范圍減小。與出口壓力301 kPa、501 kPa和1 001 kPa 相對應的平均流量分別為37.70 L/min、37.10 L/min 和36.55 L/min，可見出口壓力增加，滑片泵的平均流量減小，但減小幅度不大。圖8 為出口壓力1001 kPa 的工況，封油區(qū)打開前后的空化體積分數(shù)的分布，與圖4 相比，出口壓力提高，葉片背風區(qū)內空化潰滅加速，出口瞬時流量與轉子轉角的關系曲線變得更陡，出口瞬時流量波動增大。

圖7 出口壓力對滑片泵出口瞬時流量的影響曲線Fig.7 Curves of outlet pressures influencing the outlet transient flow rate

圖8 空化體積分數(shù)的分布Fig.8 Distribution of cavitation volume fraction

5.3 不同進口壓力對滑片泵流動特性的影響

為了研究滑片泵不同進口壓力對流動特性的影響，在進口絕對壓力為47 kPa 工況的基礎上，增加了進口絕對壓力為65 kPa、85 kPa、101 kPa 的工況。圖9 為不同進口壓力下，滑片泵出口瞬時流量與轉子轉角的關系曲線。由圖9 可知，進口壓力增加，封油區(qū)打開，從壓油腔倒流回封油區(qū)的潤滑油減少，出油口的回流減小，與回流相對應的轉角范圍減小，進口壓力101 kPa 時，出口回流基本消失。與進口壓力47 kPa、65 kPa、85 kPa 和101 kPa 的工況相對應的平均流量分別為37.10 min、42.33 L/min、48.07 L/min 和51.89 L/min，可見隨著進口壓力的提高，滑片泵平均流量增加，且增幅較大。圖10 為進口壓力101 kPa的工況，封油區(qū)打開前后的空化體積分數(shù)的分布，與圖4 相比，葉片背風區(qū)內空化區(qū)域減小，但由于出口壓力不變，因此葉片背風區(qū)內空化潰滅速率變化不大，使得不同進口壓力時，出口瞬時流量與轉子轉角的關系曲線的斜率基本一致，出口瞬時流量波動減小。

圖9 進口壓力對滑片泵出口瞬時流量的影響曲線Fig.9 Curves of inlet pressures influencing the outlet transient flow rate

圖10 空化體積分數(shù)的分布Fig.10 Distribution of cavitation volume fraction

6 結論

利用Star-cd 用戶子程序功能調用自編程序，實現(xiàn)滑片泵運動區(qū)域的網(wǎng)格重構和葉片與定子間隙的網(wǎng)格生成，以單位液體體積的空化核個數(shù)1010為空化模型的設置參數(shù)，模擬了滑片泵在不同進出口壓力狀態(tài)下的流動特性，得到如下結論:

1)在進、出口設定壓力邊界條件時，滑片泵工作過程中，由于進口流道內有一相對穩(wěn)定的空化區(qū)，瞬態(tài)進口流量波動較小，由于在吸油區(qū)和封油區(qū)內出現(xiàn)動態(tài)變化的空化區(qū)，而使得出口流量有較大的波動。

2)滑片泵進口壓力增加，葉片背風區(qū)內空化區(qū)域減小，但空化潰滅速率變化不大，使得出口瞬時流量波動減小;出口壓力增加，葉片背風區(qū)內空化區(qū)域變化不大，但空化潰滅加速，使得出口瞬時流量波動增大。

3)滑片泵進口壓力增加，出油口的回流減小，與回流相對應的轉角范圍減小，使得滑片泵的平均流量呈較大幅度的增加;滑油泵出口壓力增加，出油口的回流流量增加，與回流區(qū)相對應的轉角范圍減小，使得滑片泵的平均流量呈小幅減小。

References)

［1］Brennen C E. Cavitation and bubble dynamics［M］. New York:Cambridge University Press，2014.

［2］吳玉林. 流體機械及工程［M］. 北京:中國環(huán)境科學出版社，2005.WU Yu-lin. Fluid machinery engineering［M］. Beijing:China Environmental Science Press，2005. (in Chinese)

［3］侯訓波，孫玉清，宋昌平，等.葉片泵定子曲線及葉片對流量脈動的影響［J］. 液壓與氣動，2008，11(8):57 -59.HOU Xun-bo，SUN Yu-qing，SONG Chang-ping，et al. Analysis on cam ring and vane arousing flow fluctuation for hydraulic vane pump［J］. Chinese Hydraulics ＆ Pneumatics，2008，11(8):57 -59. (in Chinese)

［4］張海竹，盧勇，張薇，等.高壓葉片泵流體泄漏研究［J］.流體機械，2009，37(3):1 -5.ZHANG Hai-zhu，LU Yong，ZHANG Wei，et al. Analysis on the flow leakage in high pressure vane pump［J］. Fluid Machinery，2009，37(3):1 -5.(in Chinese)

［5］王力，權龍. 單作用葉片泵瞬時流量的分析與計算［J］. 振動、測試與診斷，2006，26(3):188 -191.WANG Li，QUAN Long. Transient flux analysis and calculation on single-acting vane pump［J］. Journal of Vibration，Measurement＆ Diagnosis，2006，26(3):188 -191. (in Chinese)

［6］李少年，魏列江，王崢嶸，等.油液中氣泡含量對高壓葉片泵工作腔油液壓力和葉片受力的影響［J］. 液壓與氣動，2010，13(9):91 -94.LI Shao-nian，WEI Lie-jiang，WANG Zheng-rong，et al. Influence of the bubble content in oil on working cavity pressure and force of vane in high pressure intra-vane type pump［J］，Chinese Hydraulics ＆ Pneumatics，2010，13(9):91 -94.(in Chinese)

［7］那焱青，王崢嶸，李少年，等.雙作用子母葉片泵瞬時流量的分析［J］.蘭州理工大學學報，2004，30(6):58 -60.NA Yan-qing，WANG Zheng-rong，LI Shao-nian，et al. Investigation of transient flow in double-acting pump with composite vanes［J］，Journal of Lanzhou University of Technology，2004，30(6):58 -60.(in Chinese)

［8］Osama A H. Theoretical modeling of sliding vane compressor with leakage［J］. International Journal of Refrigeration，2009，32(2):1555 -1562.

［9］Giuffrida A，Lanzafame R. Cam shape and theoretical flow rate in balanced vane pumps［J］. Mechanism and Machine Theory，2005，40(3):353 -369.

［10］Lu Y，Zhang W，Zhao Y Y，et al. Studies on several key problems of water hydraulic vane pump［J］. Industrial Lubrication and Tribology，2011，63(3):134 -141.

［11］Danardono D，Kim K S，Roziboyev E，et al. Design and optimization of an LPG roller vane pump for suppressing cavitation［J］.International Journal of Automobile Technology，2010，11(3):323 -330.

［12］Wang D M，Ding H，Jiang Y，et al. Numerical modeling of vane oil pump with variable displacement，2012-01-0637［R］. US:SAE，2012.

［13］Erzincanli B，Sahin M. An arbitary Lagrangian-Eulerian formation for solving moving boundary problems with large displacements and rotations［J］. Journal of Computational Physics，2013，255(15):660 -679.

［14］Liu J X，Akay H，Ecer A，et al. Flows around moving bodies using a dynamic unstructured overset-grid method［J］. International Journal of Computational Fluid Dynamics，2010，24(6):187-200.

［15］ildirim B，Lin S，Mathur S，et al. A parallel implementation of fluid-solid interaction solver using an immersed boundary method［J］. Computers ＆ Fluids，2013，86(5):251 -274.

［16］Hirt C W，Nichols B D. Volume of fluid(VOF)method for the dynamics of free boundary［J］. Journal of Computational Physics，1981，39:201 -225.

［17］Ubbink O，Issa R I. A method for capturing sharp fluid interfaces on arbitrary meshes［J］. Journal of Computational Physics，1999，153:26 -50.

［18］Jofre L，Lehmkuhl O，Castro J，et al. A 3-D volume-of-fluid advection method based on cell-vertex velocities for unstructured meshes［J］. Computers ＆ Fluids，2014，94(1):14 -29.

［19］Bakir F，Rey R，Gerber A G. Numerical and experimental investigations of the cavitating behavior of an inducer［J］. International Journal of Rotating Machinery，2004，10(1):15 -25.

［20］Ahuja V，Hosangadi A，Arunajatesan S. Simulation of cavitating flows using hybrid unstructured meshes［J］. Journal of Fluids Engineering，2001，123(2):331 -340.

［21］Fard M B，Nikseresht A H. Numerical simulation of unsteady 3D cavitating flows over axisymmetric cavitators［J］. Scientia Iranica B，2012，19(5):1258 -1264.

［22］Menter F R. Zonal two-equation k-w turbulence model for aerodynamic flows，1993-2906［R］. US:AIAA，1993.

［23］Langtry R B Menter F R. Transition modeling for general CFD application in aeronautics，2005-522［R］. US:AIAA，2005.

［24］Demirdzic I，Peric M. Space conservation law in finite volume calculations of fluid flow［J］. International Journal for Numerical Methods in Fluids，1988，8:1037 -1050.

［25］Wilcox D C. Turbulence modeling for CFD ［M］. California:DCW Industries Inc，2006.

［26］Van L B. Towards the ultimate conservative difference scheme V.a second order sequel to Godunov's method［J］. Journal of Computational Physics，1979，32:101 -136.

［27］Buffard T，Clain S. Monoslope and multislope MUSCL methods for unstructured meshes［J］. Journal of Computational Physics，2010，229(10):3745 -3776.

［28］Oliveira P J，Issa R I. An improved PISO algorithm for the computation of buoyancy-driven flows［J］. Numerical Heat Transfer，2001，40:473 -493.