葉頂間隙對葉輪內部流動與變形影響的流固耦合分析*

賀詩榕 董志強 劉 通

（太原科技大學機械工程學院 山西太原 030024）

離心空壓機是輸送和壓縮氣體的核心裝備， 近年來廣泛應用于航天航海以及汽車等領域［1］。 得益于電機技術、 密封技術與空氣軸承等關鍵技術的發(fā)展， 空壓機的工作壓比大幅度提高， 葉輪能承載的轉速達到了萬轉以上［2］。 半開式葉輪由于減少了輪罩的束縛，可以承載更高的許用圓周速度， 適合高轉速的工況，憑借這一優(yōu)勢廣泛應用于轉速高的工況環(huán)境［3］。 但沒有輪罩， 葉片與機殼間勢必存在一定的間隙， 工作時不可避免地出現氣體泄漏、 潛流損失與氣流分離等現象， 影響整個葉輪內部的氣體流動情況， 降低了壓縮機的效率。 另一方面， 空壓機在工作時處于高溫高壓高轉速的環(huán)境中， 會產生振動與結構熱變形， 若葉頂間隙過小， 很有可能使葉片與機殼發(fā)生碰撞， 影響空壓機的安全性。 因此， 在葉輪設計時確定合適的葉頂間隙尤為重要。

關于葉輪葉片數影響的研究方面， 靳亞峰等［4］用NUMECA 軟件分析了葉輪葉片數對離心壓縮機機組性能的影響， 發(fā)現葉片數過多會導致摩擦損失過大，機組效率下降， 葉片數過少， 葉片載荷分布增加， 二次流損失嚴重； 席光等人［5］從葉片厚度分布入手， 發(fā)現對等厚葉片沿展向和流向進行合理地削薄， 能顯著提升葉輪在大流量工況下的性能； 王藝達等［6］從葉輪葉片數、 進口傾角和子午流道型線控制點和葉片安裝角分布控制點等關鍵構型參數出發(fā)， 用數值計算方法對離心壓縮機的葉輪進行優(yōu)化。 關于葉頂間隙影響的研究方面， 姬田園等［7］對葉片前后緣的葉頂間隙進行尺寸波動干擾， 發(fā)現葉頂間隙主要對80%葉高以上部分的流場產生影響； 葉頂間隙存在偏差會產生泄漏流， 影響葉片的工作裕度； 劉玉文等［8］通過流場測量實驗和空化觀測實驗， 發(fā)現間隙流動中主要有2 種穩(wěn)定渦， 分別是泄漏渦和分離渦， 它們是引起葉頂間隙流動空化的主要因素。 針對葉輪機械在旋轉時發(fā)生的振動與變形問題， 方紹寒等［9］以離心風機葉輪為研究對象， 采用k-ε湍流模型與振動力學理論， 對葉輪結構及其流場進行分析， 得到葉輪的變形情況； 高亞飛［10］基于流固耦合理論， 運用數值模擬方法對離心泵葉輪進行仿真計算， 得到了應力應變云圖。

目前學者們針對壓氣機工作穩(wěn)定裕度范圍小、 效率低等問題已進行了深入研究， 但針對葉輪葉頂間隙大小對空壓機可靠性的影響研究較少。 在實際工程應用中， 需要保證葉輪運行穩(wěn)定， 避免產生安全隱患。所以葉輪的可靠性分析必不可少， 特別是在進行葉輪的設計及優(yōu)化時， 效率及穩(wěn)定性都要考慮在內。 董振等人［11］以軸向迷宮密封為例， 基于流固耦合理論，對迷宮密封的流場分布與結構進行數值分析與優(yōu)化。傳統的可靠性分析只考慮葉輪的離心載荷， 這與真實工況環(huán)境相差甚遠。 張衎和梁尚明［12］采用CFD 技術對葉輪內的流道進行數值模擬， 研究表明工作壓力越高， 流體載荷越大。 李彥啟等［13］以鑲嵌式機械密封為研究對象， 通過對結構受力分析和熱傳導方程， 建立熱力耦合仿真模型并進行求解， 發(fā)現熱應力對端面造成的變形量大于結構應力。

考慮到葉輪葉頂間隙的大小對空壓機可靠性的影響研究較少， 且在進行數值計算時容易忽略氣動載荷與熱載荷對葉輪的影響， 本文作者建立空壓機葉輪的三維模型， 在不同葉頂間隙下進行流場的數值模擬計算； 基于流固耦合理論建立熱流固模型， 進行葉輪的動力學分析； 綜合探究離心力、 氣動力和熱應力3 種載荷對葉片葉頂處應力和變形的影響， 為葉片的結構優(yōu)化提供理論參考。

1 計算方法

1.1 理論模型

1.1.1 流動方程

采用流體運動的控制方程， 即雷諾方程， 來描述流體的運動規(guī)律［14］。 質量方程、 動量方程、 能量方程分別為

式中：ρ為流體密度， kg／m3；t為時間， s；U為速度矢量；u、v、w分別為速度矢量在x、y、z方向的分量；μ為動力黏度， Pa·s；Su、Sv、Sw為動量守恒方程的廣義源項；cp為流體比定壓熱容， J／（kg·K）；T為流體的溫度， K；k1為流體的傳熱系數，W／（m·K）；ST為黏性耗散項。

1.1.2 湍流模型

由于CFX 具有豐富的物理模型， 文中采用CFX進行葉輪的流場模擬實驗， 用葉輪流體域的雷諾數來確定其流動情況。 雷諾數定義為

式中：ρ為流體密度， kg／m3；μ為流體動力黏度， MPa·s；v為速度， m／s；L為流體域的長度， m。

將葉輪物理模型幾何參數代入式（6） 計算可得，Re＝1.293×30×0.024 4／（18.16×10-6）＝51 189?2 300， 由此可確定流體的運動狀態(tài)為湍流。

由于k-ε湍流模型具有很好的預測能力， 又易于收斂， 適合絕大多數的工程湍流模型［15］， 所以文中的流動模型選用CFX 默認的k-ε模型。 其中，k為湍動能， m2／s2；ε為湍動能耗散， m2／s3。

1.1.3 結構力學方程

由經典力學可知， 物體動力學通用方程為

式中：M是質量矩陣；C是阻尼矩陣；K是剛度矩陣；是加速度矢量；是速度矢量；x是位移矢量；F（t）是力矢量。

1.2 物理模型

利用Cfturbo 軟件對離心空壓機葉輪進行三維建模， 葉片幾何參數如表1 所示， 葉輪實體模型如圖1所示。

圖1 葉輪實體模型Fig.1 Solid model of impeller

表1 葉輪幾何參數Table 1 Geometrical parameters of impeller

葉頂間隙指葉片的頂部與機殼之間的距離， 靠近葉片前緣的間隙稱為徑向葉頂間隙， 靠近后緣的間隙稱為軸向葉頂間隙。 文中設徑向與軸向的間隙值相同， 即S徑向＝S軸向＝S， 離心空壓機葉輪橫截面如圖2所示。

圖2 葉輪橫截面Fig.2 Impeller cross section

1.3 流固耦合

采用分離求解（不同的求解器計算各自的物理變量， 共同變量采用異步傳遞的方式進行更新） 的計算方式， 以及單向的數據傳遞方式， 將流體求解器（CFX） 計算出的壓力、 溫度數據作為載荷傳遞到固體上， 由固體求解器（Mechanical） 計算出位移應力等信息。 對離心式壓氣機葉輪進行單向流固耦合求解分析如圖3 所示。

圖3 單向流固耦合流程Fig.3 Flow of one way fluid-structure coupling

1.3.1 固體控制方程

由牛頓第二定律得：

溫差引起的熱變形為

式中：ρs為固體密度， kg／m3；as為固體域當地加速度矢量；σs為柯西應力張量；fs為體積力矢量；αT為熱膨脹系數，℃-1。

1.3.2 流固耦合方程

由于流固耦合數值計算時流體模型和結構模型是分開建立的， 所以耦合面上須滿足以下條件［16］：

式中：τ為應力， Pa；d為位移， m；Φ為熱流量， W； 下標f 表示流體， s 表示固體。

1.4 網格無關驗證及邊界設置

1.4.1 網格劃分與無關性驗證

由于葉輪幾何模型復雜， 曲率變化劇烈， 為了獲得較高的網格質量， 文中采用ICEM 對葉輪的流體區(qū)域進行非結構化網格劃分。 選用四面體劃分網格， 整體網格單元尺寸因子設為3， 葉片、 輪緣以及輪轂處的網格大小設為1 mm。 由于文中主要研究葉頂處的流動情況， 故將葉片頂部的網格進行加密處理， 設置單元大小為0.1 mm。 考慮到葉頂間隙較小， 葉輪出口處的網格也進行加密， 設置單元大小為0.1 mm。圖4 所示為葉輪流體域網格分布。

圖4 離心葉輪流體域網格分布Fig.4 Grid distribution of fluid area of centrifugal impeller

為保證結果的準確性， 以S＝0.5 mm 的葉輪為例， 進行網格無關性分析。 如圖5 所示， 當網格數量達到100 萬后， 網格數對葉輪的等熵效率和多變效率影響較小。 在保證計算精度的基礎上， 考慮到節(jié)約計算資源， 最終選取網格單元數為966 759， 總節(jié)點數為17 338。

圖5 網格無關性驗證Fig.5 Grid independence verification

1.4.2 邊界條件

工質選用可壓縮性理想空氣， 湍流模型為k-ε模型， 參考壓力為0.1 MPa， 壁面與葉輪表面采用絕熱無滑移條件， 內部傳熱設置為Total Energy。 進口總壓為一個標準大氣壓 （0.1 MPa）， 進口溫度為288.15 K， 出口壓力為0.3 MPa。 求解設置采用物理時間步長， 迭代步數為200。 葉片選用高強度的鈦合金材料， 密度為4 620 kg／m3， 屈服應力為930 MPa。

2 氣動力學分析

2.1 氣動性能分析

計算中保持葉輪的基本形狀不變， 僅改變葉頂間隙的大小。 設定葉輪的旋轉速度為120 000 r／min， 分別選取S＝0.1、 0.3、 0.5、 0.7、 0.9、 1.1 mm 進行計算， 得到不同葉頂間隙下的葉輪壓比與效率， 如圖6 所示。

圖6 葉輪氣動性能Fig.6 Aerodynamic performance of impeller

從圖6 中可以看出， 葉頂間隙在由0.1 mm 增大到1.1 mm 的過程中， 葉輪的壓比呈現上升態(tài)勢， 由3.006 增長到3.27， 1.1 mm 時的壓比值比設計壓比增加了9%， 并且還有向上增長的趨勢； 葉輪的等熵效率呈逐步下降的趨勢， 由92.42%降低到81.25%。間隙從0.1 mm 增加到0.7 mm 的過程中， 效率下降幅度比較明顯， 而從0.7 mm 增大到1.1 mm 時， 葉輪的效率下降曲線較之前相對平緩， 說明效率下降程度與葉頂間隙大小不是呈線性的變化關系。 上述結果表明葉頂間隙的大小對空壓機效率的影響非常大， 隨著間隙值的增大， 氣動性能對其不再敏感， 即S存在最優(yōu)解。

2.2 流動特性分析

2.2.1 流線分析

設計工況下離心葉輪通道內流線如圖7 所示。 可見， 紊流首先出現在葉輪的進口處， 隨著間隙值的增大， 葉輪入口處的氣流流動越發(fā)不均勻， 致使入口處流場紊亂， 這部分氣流進入主流道內， 造成主流流場失穩(wěn)。 這對壓氣機的效率以及運行時的可靠性都產生了一定的威脅。S＝0.1 mm 時的流動最為平穩(wěn)， 速度分布均勻。 與之相比，S＝0.3 mm 的流場圖入口處的流速相對均勻， 但是原本靠近葉片前緣的高速氣流變?yōu)槲蓙y的低速氣流， 這不利于整機的流動。 間隙越大， 入口紊流的區(qū)域也隨之變大。S＝1.1 mm 時， 入口葉片前緣處出現了較多的渦流， 渦流區(qū)占據了葉輪入口流體區(qū)域的50%。

圖7 葉輪整體流線Fig.7 Integral streamline of impeller： （a） S＝0.1 mm； （b） S＝0.3 mm； （c） S＝0.5 mm；（d） S＝0.7 mm； （e） S＝0.9 mm； （f） S＝1.1 mm

由于間隙的存在， 隨著葉片的旋轉方向， 在葉頂處產生泄漏流， 該氣流分為兩部分。 一部分低速氣流從葉片前緣至中部的壓力面透過葉頂進入葉片吸力面的前半部分， 形成渦流， 在吸力面一側形成很明顯的低速區(qū)。 隨著間隙值增大， 葉片頂端附近的泄漏流動影響范圍擴大， 即低速區(qū)越大。 這是因為間隙越大，泄漏流的流量越大造成的。 另一部分從葉片中后緣透過葉頂進入相鄰流道內， 與其主流線發(fā)生摻混。 隨著間隙增大， 泄漏流與主流的摻混更為劇烈， 范圍也越大。 這兩部分的泄漏流在葉輪主流道內造成了氣流堵塞， 產生了較大的流動損失， 但在一定程度上提高了主流區(qū)氣流的壓力， 提高了壓比。 當S＝0.1 mm 時，葉片頂部幾乎無泄漏流。 葉輪內流線平穩(wěn)， 且都為高速氣流， 流動狀態(tài)良好， 這也是該葉輪效率最高的原因。 間隙值增大到1.1 mm 時， 泄漏量最大， 這兩部分的泄漏流都通過葉頂間隙流入相鄰流道， 與主流摻混， 混流區(qū)域占據總流道的1／2 左右。

2.2.2 葉片載荷分析

葉片吸力面與壓力面的壓差為葉片的載荷， 它決定了葉輪的做功能力。 圖8 所示為不同葉高的葉片靜壓分布。

圖8 不同葉高處葉片靜壓分布Fig.8 Static pressure distribution of blade at different leaf heights： （a） 20% leaf height； （b） 50% leaf height； （c） 80% leaf height

由圖8 （a） 可知， 不同葉頂間隙時的葉片載荷，在20%葉高處幾乎相同， 說明葉頂間隙對葉高20%以下的氣流影響不大。 由圖8 （b） 可知， 在50%葉高處，S＝0.3 mm 時， 葉片前半部分（0.0 ～0.25 弦長） 的葉片載荷比S＝0.1 mm 時的葉片載荷小很多，這是因為間隙0.3 mm， 葉輪入口處發(fā)生紊流， 而這部分氣流隨葉片進入主流道， 擾亂了葉片壓力面上的流動， 使得載荷降低。 在0.5 弦長之后， 葉片的載荷緩慢上升， 葉輪整體的做功能力略微減弱。 由圖8 （c）可知， 80%葉高處的載荷情況與前者相比能更加直觀準確地反映出葉頂間隙的大小對葉片做功能力的影響， 隨著葉頂間隙的增大， 葉片壓力面和吸力面的壓差整體都在減小， 葉輪的做功能力明顯下降。 特別是S＝1.1 mm 時， 葉片上出現負載荷， 不利于葉輪的做功。

2.2.3 速度分析

80%葉高處不同間隙的葉片附近絕對速度的分布云圖如圖9 所示。 當間隙為0.1 mm 時， 在葉片吸力面后緣發(fā)現小范圍低速區(qū)， 這個現象一部分是由于附面層分離所導致的， 另一部分是間隙泄漏流與主流區(qū)摻混在流道前中部形成的低速渦流。 隨著間隙的增大， 低速區(qū)面積范圍擴大， 通道內的有效通流面積逐漸減小。 以上分析表明， 間隙越大， 泄漏量越多， 徑向間隙過大造成的低速氣流積聚在吸力面， 增大了能量損失， 降低了葉輪效率， 所以減小徑向葉頂間隙可以有效減少流道內渦流， 降低流動損失。

圖9 80%葉高處葉片絕對速度分布Fig.9 Blade absolute velocity distribution on 80% blade height

2.2.4 馬赫數分析

圖10 所示為葉輪子午面馬赫數的分布情況。 可知， 所有葉輪的馬赫數都不超過1， 所以在葉輪中不會產生激波， 使得邊界層分離而影響空壓機的可靠性。 當S＝0.1 mm 時， 流體運動相對平緩， 馬赫數分布均勻。 葉片前緣處馬赫數較大， 是因為氣流在靠近葉片前緣時出現膨脹加速。 在輪緣中部出現由于急轉彎， 流動擴壓度突然增大導致的氣流分離從而產生的低速區(qū)， 而在輪轂處則出現由于附面層分離產生的低速角區(qū)。 隨著間隙值的增加， 角區(qū)逐漸消失。 說明泄漏流在一定程度上可以有效帶動輪轂處的低速氣流運行， 緩解此處的流動情況。S＝0.3 mm 時， 由于泄漏流的影響， 輪緣側流動分離區(qū)域顯著擴大， 并且隨著間隙的增大， 分離區(qū)也隨之增大，S＝1.1 mm 時的分離區(qū)約占整個子午面流道的35%。 這表明葉頂處的泄漏流由葉片壓力面通過葉頂間隙向吸力面擴散， 且間隙越大， 擴散程度越大。

圖10 子午面馬赫數分布Fig.10 Mach number distribution on meridianal plane

3 結構力學分析

3.1 載荷對葉片形變以及應力大小的影響

為研究離心力、 氣動力以及熱載荷對葉片形變的影響程度， 取S＝0.5 mm 時的葉片進行仿真， 得到如表2 所示的數據。

表2 葉片的形變與應力Table 2 Deformation and stress of blade

由表2 可知， 氣動力與熱載荷造成的形變分別約占葉片總形變量的5%和10%， 所以在數值計算中，求解結構的變形以及應力時這二者的大小不可忽略，也體現了文中計算葉片變形時利用流固耦合方法的必要性。 對比第1 組和第4 組數據， 發(fā)現離心力和氣動力的載荷導致的葉片形變量還沒有單個離心力造成的形變大， 是因為氣流運動方向與葉輪旋轉方向相反，氣動力導致的形變方向與離心力相反。 由于葉片的最大變形量為0.187 2 mm， 遠小于葉片的幾何尺寸，并且氣動力和熱載荷產生的等效應力很小， 分別為23.6 和30.7 MPa， 說明雙向耦合結果與單向耦合結果相差不大， 這也是文中選用單向流固耦合的原因。葉片的形變所引起的流場變化對模擬計算的結果影響較小， 可以忽略。 考慮到單向耦合問題不涉及流體域的變化， 所以文中沒有采用動網格。

3.2 間隙對葉片形變以及應力大小的影響

圖11、 12 所示為葉輪葉片在3 種載荷共同作用下的總變形以及等效應力云圖。 由圖11 可知， 不同葉頂間隙的葉片變形趨勢一致。 較大變形主要集中在葉片前緣葉頂處， 說明葉片的高度對葉片的形變位置沒有影響。 由于輪盤的固定支撐作用， 變形程度沿著周向位置向葉根部分逐漸降低。S＝0.1 mm 時的最大變形量達到了0.195 68 mm， 占前緣總葉高的1.25%， 超過了設定的葉頂間隙值， 葉輪運行時會與機殼發(fā)生碰撞， 所以葉輪在設計時要充分考慮葉頂間隙的大小。

圖11 葉片的位移分布Fig.11 The deformation distribution of blade： （a） S＝0.1 mm； （b） S＝0.3 mm； （c） S＝0.5 mm；（d） S＝0.7 mm； （e） S＝0.9 mm； （f） S＝1.1 mm； （g） new blade

從圖12 可看到， 葉片前緣葉頂處的應力最小，葉高30%～40%處的應力最大， 出現明顯的應力集中部位。 鈦合金的屈服極限強度為930 MPa， 根據強度屈服理論：

圖12 葉片的等效應力分布Fig.12 Equivalent stress distribution of blade： （a） S＝0.1 mm； （b） S＝0.3 mm； （c） S＝0.5 mm；（d） S＝0.7 mm； （e） S＝0.9 mm； （f） S＝1.1 mm； （g） new blade

式中：σmax為最大應力， MPa； ［σs］ 為屈服極限，n為安全系數， 一般取2～5。

所以葉片上的最大應力小于等于465 MPa 較為安全可靠。

圖12 中S＝0.1 mm 與S＝0.3 mm 時， 葉片的最大應力分別為481.73 和492.7 MPa， 都大于465 MPa的安全值， 說明該類葉輪設計不合理。 由此可以在應力集中處對葉片表面進行強化處理或適當地增加葉片的厚度。 為了不增加葉片制造成本， 文中將在應力集中處適當地增大和減小葉片的厚度， 以達到對葉輪優(yōu)化設計的目的。

為了分析流固耦合作用下葉片頂部的變形特征，提取不同葉高葉片頂部的變形與應力信息繪制成曲線圖， 如圖13、 14 所示。 由圖13 可以看出， 葉片頂部變形最大的區(qū)域在弦長0.35 左右， 葉片吸力面與壓力面的變形趨勢相同， 發(fā)生較大變形的位置相同， 但最大變形量不同， 吸力面比壓力面變形大3.7%左右； 隨著間隙值的增大， 頂部的變形隨之減少。 由圖14 可以看出， 在葉片頂部， 前緣的應力最小， 其壓力面與吸力面的等效應力分布存在差異， 吸力面應力最大在弦長0.75 處， 壓力面應力最大值在弦長0.55處； 隨著間隙值的增大， 頂部的應力值整體逐漸降低。 以上分析表明， 在設計葉片時應適當降低葉片0.3～0.4 弦長處的葉高， 削薄葉片前緣的厚度， 增加0.55～0.75 弦長處的葉高， 增加30%～40%葉高處的葉片厚度。

圖13 葉片頂部變形分布曲線Fig.13 The deformation distribution curves on blade tip

圖14 葉片頂部等效應力分布曲線Fig.14 Equivalent stress distribution curves on blade tip

4 線型優(yōu)化

在保證葉輪氣動性能的基礎上， 從葉輪的安全性出發(fā)， 結合上文結構力學分析得到的結論， 對原始葉片進行模型的改進： （1） 將葉片徑向間隙值設為0.3 mm， 軸向間隙值設為0.2 mm； （2） 將葉片輪轂側前后緣的厚度設為0.8 mm， 輪緣側前后緣厚度分別設為0.5、 0.7 mm； （3） 具體厚度分布情況由具體點用貝塞爾曲線連接而成， 如圖15 所示。 將改進后的葉片重新建模， 再進行仿真分析。

圖15 葉片厚度分布曲線Fig.15 Blade thickness distribution curves

數值計算得到新葉片的等熵效率為90.416%， 壓比為3.08， 達到了工程要求。 由圖11—14 可以看出，新葉片上的最大變形為0.189 16 mm， 且主要分布在前緣， 在給定的間隙值0.3 mm 之內， 運行時不會與機殼發(fā)生碰撞， 保證了機器的安全性； 葉片上的最大應力為382.25 MPa， 與S＝0.3 mm 葉片上的應力相比， 降低了22.42%， 大大提高了機器的可靠性。

5 結論

（1） 葉頂間隙的大小對空壓機效率影響非常大，徑向間隙會造成葉輪進口處的氣流紊亂， 也會產生葉片前緣處的泄漏流， 二者共同作用使主流流場失穩(wěn)；軸向間隙產生葉片后緣處的泄漏流， 與相鄰流道的主流線發(fā)生摻混， 進一步增加了流動損失。

（2） 氣動力與熱載荷造成的形變分別約占葉片總形變量的5%和10%， 所以求解結構變形以及應力時二者不可忽略； 氣動力造成的形變方向與離心力相反。

（3） 不同葉高的葉片變形趨勢基本一致， 較大變形區(qū)主要集中在葉片前緣葉頂弦長0.35 處， 其變形量的大小隨葉高減小而減小， 變形程度沿著周向位置向葉根逐漸降低； 葉片葉高30%～40%處的應力最大， 出現應力集中。 葉頂前緣處的應力最小， 吸力面與壓力面的應力最大值分別在弦長0.75 和0.55 處，且頂部的應力值隨間隙的增大逐漸降低。

（4） 通過對葉片頂部的流場和結構力學分析，確定了最佳間隙值， 給出了葉片的線型優(yōu)化方案。 優(yōu)化模型的結果顯示葉輪氣動性能達到了工程要求， 葉片變形量在允許范圍內， 整體應力值得到顯著降低，為空壓機葉片的工業(yè)設計提供了參考。