液體火箭發動機推進劑泵誘導輪與離心輪的匹配

楊寶鋒，李斌，陳暉，劉占一

1. 西安航天動力研究所 液體火箭發動機技術重點實驗室，西安 710100 2. 航天推進技術研究院，西安 710100

作為液體火箭發動機核心部件之一，渦輪泵主要用于推進劑的輸送及增壓。隨著中國大推力補燃循環火箭發動機的研制以及發動機性能提升的迫切需求，渦輪泵的性能以及運行穩定性受到越來越多的重視。

與普通民用離心泵相比，火箭發動機推進劑泵轉速較高、結構更為復雜，在葉輪入口處容易發生汽蝕，通常采用加裝前置誘導輪來改善泵組的汽蝕性能，然而誘導輪與離心輪匹配不佳將會引起泵性能惡化以及流動不穩定現象的發生，對發動機工作的可靠性產生威脅[1]。近年來，針對離心泵性能及流動不穩定問題，國內外學者已做了大量的研究工作，主要集中在不同工況及結構參數對離心泵性能及壓力脈動的影響，獲得了豐富的研究成果。郭曉梅等[2]對有無誘導輪以及誘導輪結構變化對離心泵汽蝕性能的影響進行了研究，得到了誘導輪汽蝕嚴重性與離心葉輪汽蝕嚴重性并非成正比的結論。王洪杰等[3]對渦輪泵0.75倍額定流量工況下的壓力脈動特性進行了研究，指出誘導輪與離心輪間隙為5 mm左右時能改善該工況下的異常振動。Stel等[4]采用數值方法研究了不同流量以及轉速對兩級離心泵性能的影響，并給出了揚程與轉速、流量之間的關系式。Al-Qutub等[5]試驗研究了葉輪葉片V型出口對壓力脈動的影響，指出采用V型出口葉片后壓力脈動降低30%以上，而揚程降低5%；Barrio等[6]利用數值方法研究了葉輪蝸殼間隙大小對離心泵壓力脈動以及徑向力的影響，表明間隙減小使得壓力脈動及徑向力顯著增加。Zhang等[7-9]采用數值及試驗方法對某型低比轉速離心泵內不穩定流動與壓力脈動關系進行了研究，并對不同葉片尾緣形狀對壓力脈動影響進行了分析，指出壓力脈動幅值與渦量分布密切相關。Long等分別采用數值[10]和試驗[11]方法對非均勻入流下核反應堆冷卻泵非定常特性進行了研究，表明入口不均勻流動對泵揚程以及壓力脈動具有顯著影響，在離心泵的設計中應當給予考慮。然而上述研究主要集中在低轉速普通離心泵，對高速離心泵尤其是時序效應對其性能及穩定性方面的研究還比較匱乏。

時序效應的研究始于渦輪及壓縮機領域，主要研究轉子-轉子或靜子-靜子之間周向相對位置變化對渦輪以及壓縮機氣動性能的影響，獲得了較多的研究成果[12-14]。然而時序效應在水力機械領域的研究起步較晚，主要集中在導葉/隔舌時序效應[15-19]以及多級泵級間葉輪時序效應這兩方面的研究[20-22]，對離心泵誘導輪與離心輪周向匹配產生的時序效應研究很少見到。徐成波[23]和潘中永[24]等對高速離心泵誘導輪離心輪匹配關系進行了研究，但其針對的是能量匹配以及葉片安裝角的汽蝕性能的影響，并未對兩者周向相對位置變化引起的時序效應進行研究。盧金玲等[25]對誘導輪離心輪時序效應進行了初步研究，但其研究模型較為簡單，且轉速較低，對于復雜高速離心泵(如火箭發動機渦輪泵)，其參考價值不高。

本文以中國某型液體火箭發動機渦輪氧泵為研究對象，基于計算流體力學(CFD)全流道數值仿真結果，對誘導輪離心輪匹配的時序效應對泵外特性及壓力脈動的影響進行了研究，研究結果可為高速離心泵減振以及性能提升提供指導。

1 計算模型和數值方法

1.1 物理模型

本文研究模型為全尺寸渦輪氧泵，其幾何參數如表1所示，轉速為18 000 r/min，流體介質為低溫液氧，溫度為90 K，密度為1 086.9 kg/m3，黏性系數為1.5×10-4Pa·s。為保證仿真結果準確可靠，考慮前后泄漏流域，建立離心泵全流場仿真模型如圖1所示，共包含入口域、誘導輪域、葉輪域、擴壓器域、蝸殼域和前后泄漏域以及出口管道8個流域。此外，為消除進出口邊界擾動的影響，將泵入口及出口管道沿直線延長一段距離。

表1 泵幾何參數Table 1 Geometrical parameters of pump

圖1 計算域Fig.1 Computational domain

為研究誘導輪離心輪時序效應的影響，定義誘導輪葉片尾緣與離心輪主葉片前緣夾角沿軸向投影為匹配角度θ。考慮到周向匹配的循環對稱性，在60°對稱周期內，平均選取8個周向位置(θ=0°, 7.5°, 15°, 22.5°, 30°, 37.5°, 45°, 52.5°)研究時序效應對泵外特性的影響，其中0°、 15°、 30°和45°這4個角度用于研究時序效應對非定常壓力脈動的影響。

1.2 網格生成與計算設置

利用ICEM軟件對各流域進行六面體結構化網格劃分，以提高計算精度及收斂性，最終獲得的離心泵全流域如圖2所示。對各壁面區域進行加密，使得離心輪葉片以及擴壓器葉片等關鍵壁面平均y+<10，其余壁面平均y+<300，以滿足計算要求。采用4套網格方案進行網格無關性驗證，各流域網格數見表2。額定工況下，各方案下泵效率計算結果如圖3所示，可以看出，當網格數超過3 000萬時，計算效率基本一致，其中方案2與方案4計算結果誤差僅0.25%。為能夠更為準確地捕捉流場細節，最終選取方案3網格(網格數為6 655.65萬)進行非定常仿真計算。

圖2 全流域網格軸向截面圖Fig.2 Grid of whole domain in axial cross-section

表2 各流域網格數Table 2 Grid number of each part

流域網格數/106方案1方案2方案3方案4入口0.62810.93131.37822.2114誘導輪2.76146.798011.096116.7832離心輪7.13939.913322.326432.1822擴壓器2.10854.79709.533914.2709蝸殼2.81663.426411.128118.3977前泄漏1.72352.25414.22208.9803后泄漏2.44253.16126.571213.3173出口管0.30060.30060.30060.3006總計19.920531.581966.5565106.4436

圖3 網格無關性驗證Fig.3 Grid independence validation

采用商業軟件ANSYS CFX 17.2對三維全流道進行數值計算。對于定常仿真，采用雷諾平均Navier-Stokes(RANS)方法進行求解，湍流模型選取SST (Shear Stress Transport)k-ω模型，壁面處采取Automatic Wall Function進行處理，動靜耦合交界面采用Frozen Rotor模型進行數據傳遞，收斂精度設定為1×10-5。針對復雜模型難收斂現象，首先以一階格式進行計算，控制物理時間尺度使計算結果收斂；以收斂結果為初值，選取高精度格式繼續計算，直至結果收斂。對于非定常仿真，以定常收斂結果作為初始邊界，采用DES (Detached Eddy Simulation)方法進行求解，壁面處采用RANS方法求解，以避免實際工程問題中LES (Large Eddy Simulation)方法在壁面處網格量要求過大的限制，主流區域采用LES方法進行求解，以更好地模擬復雜流動，捕捉流場細節，動靜耦合面采用Transient Rotor Stator模型進行模擬。為保證非定常仿真結果的可靠性，時間步設置為Δφ= 1°，即每個旋轉周期對應360個時間步，計算進行20圈以獲得可靠的收斂結果，取最后5圈結果用于非定常結果分析。邊界條件根據渦輪泵真實工作狀態測量值分別定義為總壓入口和質量流量出口，各壁面給定無滑移邊界條件。

仿真計算在航天推進技術研究院高性能仿真平臺進行，每個算例使用5個節點共80個 CPU核數，非定常仿真每個算例耗時約1 000 h。

1.3 仿真結果驗證

仿真結果通過兩部分試驗進行驗證。其中泵外特性仿真結果通過渦輪泵水力試驗結果進行對比驗證；非定常壓力脈動仿真結果與發動機熱試車相應測點(圖1測點)測量結果(采樣頻率為25 600 Hz)進行對比驗證，其中誘導輪離心輪安裝角度接近30°。

渦輪泵水力試驗在西安航天動力研究所水力試驗中心進行，試驗對象為全尺寸渦輪泵，試驗介質為常溫水，水試轉速為9 000 r/min，測量7種不同流量下泵的揚程及效率，并通過相似準則轉換到額定轉速下與仿真結果進行對比。

仿真與試驗獲得的氧泵外特性曲線如圖4所示(圖中Q/Qd為泵內流量與額定工況下的流量之比)。可以看出，整個流量范圍內，仿真結果與試驗結果吻合較好，總體上仿真結果高于試驗結果，在低工況時誤差較大，在額定工況附近誤差較小，其中額定點揚程、效率誤差分別為1.52%和2.73%，表明額定工況點泵外特性仿真結果的可靠性。

圖5給出了仿真以及試驗中氧泵出口管道測點(OD1)壓力脈動時域結果。其中仿真結果為5個周期數據，試驗結果為50個周期數據，可以看出兩者時域結果吻合較好，壓力系數峰-峰值誤差<5%。對壓力信號進行快速傅里葉(FFT)變換，其頻域結果如圖6所示(圖中f為頻率,fr為轉子轉頻)。由圖6(a)可以看出，仿真結果中壓力脈動由6倍頻和12倍頻主導，這兩個頻率由離心輪擴壓器之間動靜干涉效應引起，其中6倍頻為離心輪主葉片的通過頻率(fMBPF)，12倍頻為離心輪總葉片的通過頻率(fBPF)。而熱試車由于環境復雜，其所得壓力頻譜組成也較為復雜，除了動靜干涉的主導頻率外，還出現了1倍頻(fr)、3倍頻(3fr)等其他幅值相對較高的頻率。其中1倍頻是由于真實產品裝配誤差導致轉子偏移軸線從而破壞轉子軸對稱性導致；3倍頻的出現可由本文后續分析解釋，是由于誘導輪離心輪匹配引起。總體來說仿真結果能夠準確地捕捉動靜干涉主導頻率及幅值，其中6倍頻幅值誤差為60.4%，12倍頻幅值誤差為30.9%，考慮到試車測量環境復雜性，該誤差可以接受。證明本文后續研究匹配效應對動靜干涉壓力脈動的影響具有可靠性。

圖4 數值與試驗結果對比Fig.4 Comparison between numerical and experimental results

圖5 出口管道測點壓力脈動時域結果Fig.5 Time-domain pressure pulsation of monitor point at outlet duct

圖6 出口管道測點壓力頻譜Fig.6 Pressure spectrum of monitor point at outlet duct

2 外特性分析

2.1 匹配角度對泵性能的影響

圖7給出了額定工況下誘導輪離心輪不同匹配角度下泵揚程系數以及效率的變化曲線。可以看出，周向匹配產生的時序效應對泵揚程及效率具有一定的影響，隨著匹配角度的增加，揚程、效率均呈現先降低后緩慢增加的趨勢，兩者變化幅度分別為0.8%、1.2%。其中0°時具有最高的揚程及效率，30°時達到最低值，這表明當誘導輪葉片尾緣與離心輪葉片前緣相對時，可獲得最高的揚程及效率。

圖7 不同匹配角度下泵性能變化情況Fig.7 Pump performance variation at different matching angles

2.2 熵產分析

為闡釋匹配效應對泵性能的影響，引入熵產理論對泵內部能量損失進行分析。流場熵產分析方法由Kock和Herwig[26]提出，近兩年來才逐漸應用到水力機械領域中，取得較好的效果[27-28]。其將湍流流場損失分為湍流平均運動引起的損失以及脈動運動引起的損失兩大部分(湍流耗散損失)，兩者的計算表達式為

(1)

(2)

泵內各部件能量損失可通過對式(1)和式(2)進行體積分來獲得，即

(3)

式中：V為流體體積。

由熵產理論獲得各部件損失隨匹配角度變化曲線如圖8所示，可以看出，葉輪及擴壓器流域損失遠大于其他流域損失，隨著匹配角度增加，呈現出先增加后減小的趨勢，這與外特性曲線變化趨勢相對應。其他流域隨匹配角度變化損失變化不大，由此可以得出，誘導輪離心輪匹配的時序效應對外特性的影響主要來自離心輪及擴壓器流域流動狀態的變化。

圖9給出了誘導輪、離心輪Blade to Blade (B2B)中截面的局部熵產率(Local Entropy Production Rate, LEPR)以及流線分布圖。由圖9(a) 可以看出，離心輪內損失遠大于誘導輪內損失，當匹配角度θ為0°和15°時，高熵產區域主要分布在靠近誘導輪葉片的3個流道內(圖中橢圓標注)，而在匹配角度為30°和45°時，其余3個流道也出現較高熵產分布，因此也導致更高的損失的發生，這與圖8的損失變化趨勢相符。由圖9(b) 流線分布可知，較高的熵產分布區域對應較強的流動分離以及由此產生的分離渦，在匹配角度為0°和15°時，靠近誘導輪的3個葉片通道出現明顯的流動分離，其余3個通道流動狀態相對較好，而在匹配角度為30°和45°時，6個葉片通道均出現不同程度的流動分離以及相應的分離渦，幾乎堵塞了整個流道，最終導致揚程效率的下降。此外，由離心輪內流動狀態可知該離心輪具有較大的優化空間。

圖8 不同匹配角度下各流域的能量損失Fig.8 Energy loss of each domain at different matching angles

圖9 誘導輪、離心輪B2B中截面LEPR及流線分布Fig.9 Distribution of LEPR and streamlines in B2B cross-section of inducer and impeller

圖10給出擴壓器內熵產率及流線分布圖，由圖10(a)可知，不同匹配角度下熵產分布模式相似，其中高熵產區域主要分布在擴壓器入口處葉輪尾跡區域以及靠近隔舌處的葉片通道內(通道A，虛線橢圓標注)。隨著角度變化，高熵產區域也在變化，θ=30°時擴壓器入口處以及葉片通道A中高熵產區域顯著增大，對應著較大的損失產生。由圖10(b)可以看出，在通道A中出現明顯的回流現象，匹配角度為0°時，回流現象不明顯；在θ=30°時，回流現象顯著增強，并發展成較強的回流渦，堵塞了整個通道A，這也是高熵產區域增大，損失增加的一個重要原因；此外，在其他通道也出現流動分離等不穩定現象(圖10(b)中實線橢圓標注)，然而此區域能量損失很小，這表明傳統的利用流線分析確定流場損失的方法存在一定的缺陷。

通過熵產分析可知，誘導輪離心輪匹配對外特性的影響主要由離心輪及擴壓器內流動狀態決定，其形成機制為：不同匹配角度下，葉輪通道內分離渦、葉輪尾跡效應以及靠近隔舌的擴壓器葉片通道內回流渦的變化共同作用導致。

圖10 擴壓器B2B中截面LEPR及流線分布Fig.10 Distribution of LEPR and streamlines in B2B cross-section of diffuser

3 壓力脈動分析

3.1 壓力脈動強度分布

為了對離心泵整個旋轉周期內流場的壓力脈動強度進行評估，通過引入壓力脈動標準差與葉輪出口處動壓進行無量綱化，定義相應的壓力脈動強度系數Cpsd如式(4)所示，其優勢在于能夠獲得流場內壓力脈動強度分布情況，準確定位流場內高壓力脈動發生的位置。

Cpsd=

(4)

式中：N為一個計算周期內時間步數，即壓力采樣數；p(x,y,z,ti)為節點(x,y,z)在第i個時間步的靜壓大小；U2為葉輪出口圓周速度。

圖11為誘導輪流域B2B平面內壓力脈動強度分布，可以看出誘導輪內壓力脈動水平較低，較高的壓力脈動主要分布在葉片出口壓力面處；此外，不同匹配角度下壓力脈動強度分布相同，表明時序效應對誘導輪內壓力脈動影響較小。

圖12給出了離心輪、擴壓器以及蝸殼中截面壓力脈動強度分布情況。與圖11相比，這3個流域的壓力脈動強度水平遠高于誘導輪流域。此外，高壓力脈動區域主要集中在動靜干涉區域以及擴壓器導葉入口吸力面附近。由圖12可知，不同匹配角度對壓力脈動強度分布影響顯著，當匹配角度為0°時，壓力脈動水平最高，各導葉入口處高壓力脈動區域大小相當；隨著匹配角度變化，壓力脈動水平開始下降，尤其表現在擴壓器導葉入口處。當匹配角度為30°時，即誘導輪葉片尾緣位于離心輪相鄰主葉片中間位置時，壓力脈動水平達到最低狀態，主要表現為擴壓器導葉入口處高壓力脈動區域以及蝸殼區域壓力脈動的顯著減小，并且越遠離隔舌的導葉入口處高壓力脈動區域減小幅度越大。

圖11 誘導輪B2B中截面壓力脈動強度分布Fig.11 Distribution of pressure pulsation intensity in B2B surface of inducer

圖12 離心輪、擴壓器、蝸殼中截面壓力脈動強度分布Fig.12 Distribution of pressure pulsation intensity in impeller, diffuser and volute domain at mid-span section

3.2 壓力脈動頻譜分析

根據壓力脈動強度分析結果，對壓力脈動水平較高的動靜干涉區域以及擴壓器導葉附近壓力脈動頻譜進行分析。圖13給出了相應的壓力脈動測點分布，在動靜干涉區域沿周向平均布置10個 測點，各測點位于導葉葉片前緣附近；在靠近隔舌處的導葉布置4個測點，分別位于導葉前緣吸力面處(DF1)、導葉吸力面(DF2)與壓力面(DF4)中部以及導葉尾緣吸力面(DF3)處。此外，在出口管道處設置測點如圖1(b)所示，該測點與熱試車測點保持一致，用于數值結果與試驗結果的對比。

對各測點壓力信號進行FFT變換，為了評估壓力脈動能量在特定頻域變化趨勢，定義不同頻率離散幅值的RMS值[8]為

(5)

式中：Ai為不同頻率下壓力脈動幅值。

表3給出了匹配角度為0°和30°時動靜干涉區域各測點在0～20 760 Hz頻域內壓力脈動RMS值(無量綱結果)。可以看出，靠近隔舌的測點RS1壓力脈動水平最高，遠離隔舌的測點壓力脈動水平顯著下降，這與3.1節壓力脈動強度分析結果相符。當匹配角度為30°時，動靜干涉區域壓力脈動水平顯著降低，各測點RMS幅值均下降10%以上，平均下降14.50%。其中測點RS10降幅最大，達到18.87%。

表4給出了擴壓器表面測點在0～20 760 Hz下的壓力脈動RMS值(無量綱結果)。可以看出壓力脈動最大水平出現在DF1測點，這與3.1節壓力脈動強度分析結果高脈動區域分布在導葉入口吸力面附近相符。當匹配角度為30°時，壓力脈動水平顯著降低，RMS幅值平均下降16.7%。其中測點DF2降幅最大，達到34.76%。

圖13 壓力脈動測點設置Fig.13 Arrangement of pressure pulsation monitor points

表3 動靜干涉區域壓力脈動RMS值Table 3 RMS values of pressure pulsation in rotor-stator interaction region

表4 擴壓器葉片表面壓力脈動RMS值

Table 4 RMS values of pressure pulsation on diffuser blade surface

匹配角度/(°)壓力脈動RMS值DF1DF2DF3DF405.28581.27873.30062.4459304.86610.83422.86022.1784Difference/%7.9434.7613.3410.94

根據RMS值分析結果，選取脈動水平最大以及降幅最大的測點(動靜干涉區域RS1、RS10；擴壓器表面DF1、DF2)進行頻譜分析。

圖14給出了4種不同匹配角度下測點RS1壓力脈動頻譜。可以看出，由于動靜干涉效應，4種 匹配角度下離心輪葉片通過頻率(fMBPF、fBPF)及其倍頻起主導作用。此外，由于誘導輪3個 葉片的影響，在0°時，3倍轉頻(3fr)非常突出，其幅值與主葉片通過頻率幅值相當；隨著匹配角度變化，3倍頻逐漸較小，在30°時，3倍頻基本消失，而其他頻率幅值基本保持不變。

圖15為RS10測點在不同匹配角度下的壓力脈動頻譜。可以看出，各匹配角度下，葉片通過頻率及其倍頻起主導作用，當匹配角度為0°時，3倍頻較大，其幅值已超過葉片通過頻率幅值，隨著匹配角度變化，3倍頻逐漸減小，并在30°時消失，這與RS1測點變化規律一致。

圖14 不同匹配角度下RS1測點壓力脈動頻譜Fig.14 Pressure pulsation spectra of RS1 at different matching angles

圖15 不同匹配角度下RS10測點壓力脈動頻譜Fig.15 Pressure pulsation spectra of RS10 at different matching angles

圖16給出擴壓器葉片前緣吸力面附近DF1測點壓力脈動頻譜。該測點靠近動靜干涉區域，因此其頻譜呈現出與前述RS1、RS10相似的頻譜特性。但由于其位于較高壓力脈動區域(如圖12所示)，因此各主導頻率幅值顯著高于RS1、RS10測點。

圖17給出了擴壓器吸力面中心DF2測點壓力脈動頻譜。該測點在不同角度下主導頻率為主葉片通過頻率(fMBPF)，其諧頻成分逐漸消失。同樣當匹配角度為0°時，3倍頻成分也起到主導作用，30°時，3倍頻成分消失。由此可以得出結論，誘導輪離心輪周向匹配參數對離心泵動靜干涉效應影響顯著，當誘導輪出口葉片位于離心輪相鄰主葉片中間位置時，能夠顯著降低壓力脈動水平，其本質為降低了壓力脈動的3倍頻成分。

為對誘導輪與離心輪周向匹配引起的壓力脈動3倍轉頻成分的出現與抑制現象進行解釋，圖18(a)和圖18(b)給出了周向匹配角度為0°以及30°時離心輪內部流線分布情況。由圖18(a)可以看出，當匹配角度為0°時，靠近誘導輪葉片出口吸力面的3個流道出現了嚴重的流動分離(圖中橢圓標注部分)，幾乎堵塞了整個葉片通道，而其余3個葉片通道流動狀態相對較好，整個離心輪流域被平均分成3股循環對稱流動。因此在離心輪出口與擴壓器入口區域的動靜干涉過程中，這3股對稱流動引起了較高的3倍轉頻成分。而當誘導輪葉片尾緣位于離心輪相鄰主葉片中間位置，即匹配角度為30°時，這種對稱效應被破壞，離心輪6個主葉片通道出現了不同程度的分離現象，因此在隨后的動靜干涉過程中，3倍頻消失，離心輪葉片通過頻率(6倍頻)及其倍頻起到主導作用。

圖16 不同匹配角度下DF1測點壓力脈動頻譜Fig.16 Pressure pulsation spectra of DF1 at different matching angles

圖17 不同匹配角度下DF2測點壓力脈動頻譜Fig.17 Pressure pulsation spectra of DF2 at different matching angles

圖18 葉輪流域瞬時流線分布Fig.18 Distribution of instantaneous streamlines in impeller

4 結 論

本文采用基于DES的離心泵三維CFD全流道數值仿真方法對中國某型液體火箭發動機渦輪泵誘導輪離心輪周向匹配產生的時序效應進行了研究，得到以下結論：

1) 誘導輪離心輪周向匹配的時序效應對泵外特性有一定的影響。隨著匹配角度的增加，泵揚程及效率均呈現先減小后緩慢增大的趨勢，兩者變化分別達到0.8%和1.2%。當誘導輪葉片后緣與離心輪葉片前緣正對時，可獲得最大的揚程及效率。

2) 通過熵產分析可知，誘導輪離心輪時序效應對能量損失的影響主要集中在葉輪以及擴壓器流域。時序效應對外特性的影響機制由葉輪通道分離渦、葉輪葉片尾跡以及靠近隔舌處擴壓器葉片通道回流渦的變化所決定。

3) 誘導輪離心輪時序效應對葉輪擴壓器動靜干涉效應影響顯著。當誘導輪葉片后緣位于離心輪相鄰主葉片中間位置時，可有效消除3倍頻成分，顯著降低泵內壓力脈動水平。其中動靜干涉區域以及隔舌處擴壓器葉片表面壓力脈動幅值平均下降14.5%和16.7%。