黃偉峰 王偉達 劉 瑩 劉向鋒 何 強

(清華大學摩擦學國家重點實驗室 北京 100084)

氦氣透平在高溫氣冷堆核電站[1]、航天等領域[2-3]具有重要應用。對于氦氣透平，軸封技術是其關鍵技術和難點之一，而干氣密封由于其低泄漏、高可靠性特性，是一種非常理想的氦氣透平軸密封型式。干氣密封是一種非接觸式氣體潤滑密封[4]，其工作性能與基于流場、傳熱、變形的綜合仿真分析設計能力密切相關。

自干氣密封于20世紀60年代出現(xiàn)以來，大批專家學者在干氣密封的理論建模及仿真計算方面開展了研究。相關研究大致可分為兩類。其中之一是基于等溫假設的純流場分析計算[5-8]，由于流場分析是干氣密封性能分析的基礎，因此相關研究較多。采用的方法包括用于螺旋槽密封的解析法[5]，基于有限元法[6]、有限差分法[7]等的全數值法，以及結合了解析法和離散方法特征的半解析法[8]等。另一類是考慮了傳熱過程、固體變形兩者之一或者全部的多場耦合分析，這類方法考慮的因素完整，對干氣密封性能預測更加準確，是干氣密封分析方法的發(fā)展方向。LEBECK[9]在其著作中從求解機械密封環(huán)端面位移的圓環(huán)理論出發(fā)，構建了基于線性位移分布的多場耦合模型。THOMAS等[10]通過數值求解流場、傳熱、變形場，建立了高壓干氣密封的多場耦合模型。BAI等[11]對高壓及高轉速工況下的干氣密封進行了耦合分析。FAIRUZ等[12]面向工程中的密封結構，在較為準確的幾何模型和邊界條件下進行了單向熱-流固分析。約翰克蘭等代表性國際密封企業(yè)，開發(fā)了自己的機械密封多場耦合分析軟件，但這些軟件僅在企業(yè)內部使用，軟件中采用的算法和模型的具體內容未發(fā)表在公開文獻中。總體而言，可以看到目前已有一些關于干氣密封的多物理場模型的研究，為干氣密封耦合分析方法的發(fā)展奠定了很好的基礎。但將這些模型用于干氣密封設計尚存在一些問題：大部分耦合計算的研究采用簡化邊界條件，如補償環(huán)完全固定；部分研究未納入干氣密封軸向開啟力和閉合力的平衡，在給定膜厚條件下進行計算；部分研究的能量方程求解尚存在一些不足之處等。相關簡化處理方法往往可以大大降低耦合分析的收斂難度，但與干氣密封實際受力狀態(tài)并不一致，使這些模型在工程設計中的應用受到限制。

本文作者以高溫氣冷堆用氦氣干氣密封為研究對象，在氦氣介質條件下，建立流場、傳熱、變形及耦合過程模型，以及更接近干氣密封實際力熱狀態(tài)的邊界條件，發(fā)展干氣密封的多場耦合分析方法，并對氦氣干氣密封在不同工況參數和槽型參數下的性能變化規(guī)律進行分析，為氦氣干氣密封的機制研究和結構設計提供一定的依據和幫助。

1 數學模型和數值方法

1.1 幾何模型

文中以螺旋槽干氣密封為研究對象，圖1(a)為螺旋槽干氣密封的端面結構示意圖，各幾何參數分別為：密封環(huán)帶內徑ri=96.5 mm，外徑ro=120.3 mm，槽根半徑rg=108.4 mm，平衡半徑rb=100.5 mm，螺旋角α=15°，槽深hg=10 μm，槽數Ng=30，堰槽寬度比γ=1。密封介質為氦氣。從干氣密封結構特征出發(fā)，取旋轉環(huán)與靜止環(huán)密封環(huán)進行軸對稱幾何建模，如圖1(b)所示，其他結構件的作用以邊界條件體現(xiàn)。

圖1 螺旋槽干氣密封的端面結構示意圖及固體域分析模型

1.2 流體域方程

流場分析采用螺旋槽機械密封的半解析式方法[8]。由MUIJDERMAN解析模型[13]，螺旋槽密封的壓力分布如下：

(1)

式中：p(r)為氣膜壓力在密封環(huán)帶上的徑向分布；ωa為角速度；η為流體的動力黏度；ρ為流體密度；st為質量泄漏率；h1為槽區(qū)膜厚；h2為非槽區(qū)膜厚；C為積分常數。

干氣密封運行時，密封環(huán)間隙中氣體的密度、黏度及氣膜厚度均會沿徑向發(fā)生變化，MUIJDERMAN模型無法直接使用。為此，文中采用了基于該模型發(fā)展出的半解析式模型[8]。將密封端面流場在徑向進行離散化處理，可用于介質黏度、介質密度以及膜厚沿徑向位置發(fā)生變化等條件下的機械密封流場分析?？紤]氦氣介質的可壓縮性，引入真實氣體狀態(tài)方程[14]：

(2)

式中：ρ為氦氣的密度；p為氣體的絕對壓力；T為氣體的絕對溫度。

氦氣的動力黏度可由黏溫公式[14]確定：

η=3.674×10-7T0.7

(3)

在忽略氣膜內沿徑向和周向的熱傳導、忽略氣體體積力條件下，流體域能量方程[15]如下：

(4)

式中：e為比內能；v為速度；p為壓力；τ為黏性應力張量；k為熱傳導系數；S為體積源項。

當密封運行達到穩(wěn)定狀態(tài)時，補償環(huán)的軸向工作位置可由力平衡方程確定：

Fo=Fc

(5)

式中：Fo為開啟力；Fc為閉合力。

開啟力可通過將氣膜壓力在密封環(huán)帶上積分進行求解，其表達式為

(6)

閉合力由密封入口壓力、出口壓力與彈簧力共同組成，其表達式為

(7)

式中：po為密封入口壓力；pi為密封出口壓力；psp為彈簧比壓。

1.3 固體域方程

旋轉環(huán)和靜止環(huán)中，均采用軸對稱穩(wěn)態(tài)熱傳導方程[16]：

(8)

用于求解熱彈性固體域變形問題的控制方程[17]如下：

平衡方程為

(9)

式中：σrr、σθθ與σzz分別為徑向、周向與軸向正應力；τrz為圓柱面上的切應力；br與bz分別為徑向與軸向的單位體積力。

幾何方程為

(10)

式中：εrr、εθθ與εzz分別為徑向、周向與軸向正應變；γrz為徑向與軸向間的切應變；ur與w分別為徑向與軸向的位移分量。

物理方程為

(11)

式中：E為材料的彈性模量；ν為材料的泊松比；G=E/[2(1+ν)]為材料的切變模量；αT為材料的熱膨脹系數；Tref為熱應變的參考溫度,文中取Tref=293.15 K。

1.4 邊界條件

對于流體域，壓力邊界條件為

p(r=ri)=pi;p(r=ro)=po

(12)

溫度邊界條件為

T(r=ro)=To

(13)

式中：To為密封入口溫度。

在流體域與固體域界面處，滿足熱流密度和溫度連續(xù)條件：

Tg=Ts;qg=qs

(14)

式中：Tg和Ts分別為界面處氣體與固體的溫度；qg和qs分別為界面處氣體與固體的法向熱流密度。

對于固體域，圖2(a)所示為固體域的受力與約束邊界條件。其中密封環(huán)帶部分施加由流場分析求解得出的密封氣體壓力分布，其余部分施加入口壓力或出口壓力，動環(huán)與靜環(huán)背部靠近內側處施加軸向位移約束。動環(huán)以一定的速度隨軸旋轉，因此對其施加離心力載荷。

圖2(b)所示為固體域熱邊界條件。其中密封環(huán)帶部分施加由流體域求解得出的熱流密度分布；兩密封環(huán)背部和底部與周圍部件形成的間隙較小，氣體流動速度慢，導熱能力差，因此采用絕熱邊界條件；動環(huán)和靜環(huán)外徑處部分取對流換熱邊界條件(由于動環(huán)旋轉引發(fā)的外徑處氣體流動會作用于靜環(huán)，因此圖2(b)中靜環(huán)外徑處與動環(huán)外徑處所取對流換熱系數相同)，即

(15)

圖2 密封環(huán)邊界條件示意

式中：h*為對流換熱系數;ks為固體的熱傳導系數。

目前常用于機械密封對流換熱系數求解的經驗公式有BECKER[18]公式、DOANE[19]公式以及TACHIBANA[20]公式等。文中計算采用經典的BECKER公式，其表達式為

(16)

式中：kf為流體的熱傳導系數；D為旋轉環(huán)相應邊界處的直徑；Re=ρuD/η為以D為特征長度的雷諾數，其中ρ為流體密度，η為流體的動力黏度，u為旋轉環(huán)相應邊界處的線速度；Pr=cpη/kf為普朗特數，其中cp為流體的定壓比熱容。

1.5 耦合過程和數值方法

干氣密封運行中，流體場的壓力分布與黏性剪切生熱會導致密封環(huán)產生變形，進而改變氣膜厚度分布；而氣膜厚度分布的改變又會對流場與溫度場產生影響。為準確反映干氣密封的這種耦合效應，文中利用上述理論對干氣密封的流場、傳熱、變形進行計算，并通過界面間的物理量傳遞對三者之間進行反復迭代。其中主要包含3個迭代過程：開啟力閉合力平衡迭代，流體域固體域傳熱迭代，流體壓力與密封環(huán)變形迭代。整體計算流程和3個迭代的嵌套關系如圖3所示。

圖3 耦合計算流程

在文中研究中，流體域方程采用半解析法求解；固體域傳熱方程和變形方程均采用有限元法進行數值離散和求解。

2 結果與討論

基于文中建立的熱-流固耦合模型，首先對氦氣干氣密封在典型工況下的基本工作特性進行分析；進而分別以轉速和槽深作為工況參數和幾何參數的代表，研究熱-流固耦合條件下工況參數和幾何參數變化對密封產生的影響；對所建立的模型的能力進行檢驗，并通過分析仿真結果深化對密封工作機制和特性的認識。

2.1 基本工況分析

文中計算設定的基本工況為：入口壓力po為7 MPa；出口壓力pi為0.2 MPa；彈簧比壓psp為0.01 MPa；入口溫度To為320.15 K，轉速ω為4 000 r/min。旋轉環(huán)采用碳化硅材料，靜環(huán)采用石墨材料，其性能參數如表1所示。

表1 密封環(huán)材料性能參數

為了方便地評估密封環(huán)變形，以密封環(huán)帶內徑為參考點，定義密封環(huán)帶上任意一點的軸向相對變形量為

Δz=z-zri

(17)

式中：z為密封環(huán)帶上任意一點處的軸向位移分量；zri為密封環(huán)帶內徑處的軸向位移分量。

圖4所示為密封環(huán)端面的軸向相對變形量分布情況?？梢钥闯觯瑒迎h(huán)端面最大軸向變形量為0.222 μm，其中熱致變形沿徑向分布使密封間隙趨于收斂形(這里指由外徑向內密封間隙由大變小)，最大變形量為0.217 μm；力致變形沿徑向呈曲線分布，與熱致變形量綜合作用后的總變形沿徑向呈近似線性收斂分布。靜環(huán)端面總變形總體上使密封間隙趨于收斂形，最大軸向變形量為1.607 μm，其中熱致變形沿徑向呈近似線性分布，最大變形量為0.748 μm；力致變形為曲線分布，導致總變形在密封環(huán)帶內徑附近有輕微的彎曲，其余部分為近似線性分布。由此可見，動環(huán)的力致變形量小于靜環(huán)，這是因為碳化硅材料的彈性模量較大，因此動環(huán)抵抗壓力變形的能力更強。熱致變形會形成有利于密封穩(wěn)定性的收斂錐度，并且在密封環(huán)端面的軸向變形中占據重要地位。

圖4 密封環(huán)端面力致軸向相對變形量、熱致軸向

圖5所示為熱-流固耦合模型與純流場模型計算得到的氣膜厚度沿徑向分布情況。熱-流固耦合模型中考慮了密封環(huán)變形對膜厚的影響，最小膜厚出現(xiàn)在內徑處，大小為3.238 μm，最大膜厚出現(xiàn)在外徑處，大小為5.066 μm；而純流場模型由于無法考慮密封環(huán)變形的影響，因此氣膜厚度沿徑向保持不變，大小為2.960 μm。通過熱-流固耦合模型計算得到的密封泄漏率為5.9 m3/h(標況下體積泄漏率，下同)；通過純流場模型計算得到的密封泄漏率為4.1 m3/h。由2種模型的計算結果對比可知，密封環(huán)產生的變形會使得平均膜厚增加，從而增大密封氣體的泄漏。

圖5 2種模型的氣膜厚度分布對比

圖6所示為熱-流固耦合模型計算得到的溫度分布?？梢钥闯觯芊猸h(huán)端面溫度自氣體入口至出口先升高后降低，于內側達到極大值。這是因為在氣體強制對流換熱的作用下，密封環(huán)外徑處具有較好的散熱性能，端面產生的熱量可以及時傳出，因此溫度較低；距離入口越遠的位置散熱性能越差，因此溫度不斷升高。氣體出口附近壓力梯度較大，壓力驟降使得氣體產生明顯的膨脹吸熱現(xiàn)象，導致該區(qū)域溫度下降。密封環(huán)端面遠離入口的部分平均溫度相對較高，熱膨脹效應更加明顯，因此密封環(huán)的熱致變形呈收斂分布。碳化硅材料與石墨材料相比熱傳導系數較高，因此密封運轉中產生的剪切熱量主要由動環(huán)導出。

圖6 密封環(huán)溫度場分布

2.2 轉速對密封性能的影響

以轉速為代表量，研究了工況參數變化時密封性能和其他參數變化的情況。采用文中模型分別對密封轉速為3 000、4 000、5 000、6 000、7 000、8 000 r/min時的工況進行仿真分析，計算過程中采用的其他參數同基本工況。

圖7所示為不同轉速下的氣膜沿徑向的溫度分布?？梢钥闯?，轉速為3 000 r/min時氣膜的最大溫升為2.68 K，而轉速為8 000 r/min時最大溫升達到了10.27 K。雖然轉速升高后密封環(huán)外側的對流換熱條件會得到一定改善，但是氣體黏性剪切效應更加明顯，生熱量急劇增大，導致氣膜整體溫度顯著升高，可見黏性剪切生熱是影響干氣密封溫升的主要因素。

圖7 不同轉速下氣膜溫度的徑向分布

圖8所示為不同轉速下的動環(huán)與靜環(huán)變形分布。

圖8 不同轉速下密封環(huán)端面軸向相對變形量

從圖8可以看出，轉速為3 000 r/min時動環(huán)的最大變形量為0.146 μm，靜環(huán)的最大變形量為1.480 μm；而轉速為8 000 r/min時動環(huán)的最大變形量為0.606 μm，靜環(huán)的最大變形量為2.159 μm。這是因為在轉速較高時密封溫升較大，會造成更大的熱致變形，從而導致動環(huán)與靜環(huán)形成較大的收斂變形。

圖9所示為不同轉速下的氣膜厚度分布。轉速升高主要可能帶來幾個方面的效應：(1)螺旋槽產生的流體動壓效應變大，導致氣膜厚度增加；(2)螺旋槽對氣體的泵送能力增強，也會使泄漏增大；(3)氣膜剪切生熱增大(不考慮膜厚變化)，使得溫度和熱變形發(fā)生變化；(4)壓力分布變化，力致變形隨之發(fā)生變化。從圖9可以看出，氣膜厚度整體隨轉速升高而增大，變形形成的錐度也同步增大，這是以上4種效應綜合產生的結果。隨之帶來的密封泄漏和剛度的變化如圖10所示，可見，泄漏率隨轉速升高顯著增加，但氣膜剛度有所下降，后者主要是因為膜厚增大導致的結果。

圖9 不同轉速下的氣膜厚度分布

圖10 不同轉速下的泄漏率與氣膜剛度

2.3 槽深對密封性能的影響

以槽深為代表量，研究了槽型參數變化時密封性能和其他參數變化的情況。采用文中模型在槽深為7～12 μm時，如2.1節(jié)所述的基本工況條件下進行了仿真分析。

圖11所示為不同槽深下的氣膜厚度分布?？梢姎饽ず穸入S槽深單調增加，這是因為槽深的增加會使流體動壓效應、泵送能力、平均膜厚均增大(在文中選取參數范圍內)。泄漏率對槽深變化十分敏感，隨著槽深增加，密封泄漏率顯著增大，如圖12所示，但氣膜剛度隨轉速增高而變小，其原因也是氣膜厚度增加。

圖11 不同槽深下的氣膜厚度分布

圖12 不同槽深下的泄漏率與氣膜剛度

圖13所示為不同槽深下的氣膜沿徑向的溫度分布?？梢钥闯?，槽深為7 μm時氣膜的最大溫升為5.66 K，而槽深為12 μm時最大溫升僅2.91 K。這是因為增加槽深會導致氣膜整體厚度的增大，氣體黏性剪切生熱量隨之降低，使得氣膜溫度下降；同時泄漏量的增加也會帶走更多熱量。

圖13 不同槽深下氣膜溫度的徑向分布

圖14所示為不同槽深下的動環(huán)與靜環(huán)變形分布?？梢钥闯觯瑒迎h(huán)變形隨槽深增大逐漸減小，靜環(huán)變形隨槽深增大不斷增大，但動、靜環(huán)的整體變形量總體變化較小，對膜厚形狀影響不大。

圖14 不同槽深下密封環(huán)端面軸向相對變形量

3 結論

(1)考慮溫度分布、密封環(huán)變形與氣體性質、膜厚分布等因素之間的相互作用關系，建立了螺旋槽干氣密封熱-流固耦合模型，可用于泄漏率、氣膜剛度、膜厚分布以及溫度分布等密封性能參數的預測，可為干氣密封結構設計提供有力的技術支撐。

(2)利用熱-流固耦合模型與純流場模型分別對典型工況下的密封性能進行了計算分析，結果表明，在文中所給參數范圍內，密封環(huán)的熱、力綜合變形會明顯改變間隙形狀，從而對泄漏率等密封性能造成顯著影響；而純流場模型無法考慮密封環(huán)變形及溫度梯度的影響，會導致計算結果產生較大誤差。

(3)對不同轉速及槽深下的密封性能進行了計算分析，結果表明，在文中所給參數范圍內，轉速及槽深的變化會直接對流場動壓及傳熱性能產生影響，并通過熱、流、固三場之間的相互耦合作用影響密封性能。