非線性動態特性系數對汽輪機轉子-軸承-密封系統運動特性的影響

曹麗華, 薛 川, 司和勇, 高路路, 李 想, 郝德成

(東北電力大學 能源與動力工程學院, 吉林 吉林 132012)

密封可以大幅度提高超超臨界汽輪機的經濟性,但密封腔內也會產生非線性特征明顯的汽流激振力影響轉子-軸承-密封系統的運動特性與穩定性[1-3]。為更真實的反應轉子的運動特性,系統內的非線性動態特性以及轉子工作過程中熱、動載荷引起的運動特性變化不容忽視。

Alford[4-5]首次發現汽流激振現象,并提出Alford力計算公式,早期的相關研究主要是圍繞著Alford力公式分析。為更準確地得到計算模型,研究發現建立密封控制體模型可以提高計算的準確性[6-7]。Muszynska模型雖融入了流體慣性但忽略了其軸向速度,且非線性不強難以計算偏心渦動較大的情況[8-9]。國內學者先后提出攝動算法與振蕩流體力學法對轉子動力特性的研究提供較大幫助[10-11]。隨著數值模擬的發展,模擬精度大幅度提高可以更真實地反映密封腔內流場運動狀態。丁學俊等[12]利用靜偏心模型考察了轉子偏心距、轉速對激振力的影響。陳堯興等[13]利用動網格模型分析了汽輪機非均勻進氣對轉子動力特性系數的影響。通過理論求解與數值模擬的研究,汽流激振的機理已十分清晰,研究重點已轉移到轉子動力特性的分析上。曹麗華等[14-15]結合試驗法對汽輪機末級葉片模型進行模態分析,得到葉片有無裂紋對轉子固有頻率以及振型的影響,以及高壓缸的運動特性。Kirk等[16-17]利用CFX計算了密封-轉子系統泄漏量、壓力分布以及密封激振力,并給出了密封齒的結構參數對系統動力特性的影響。司和勇等[18]分析了汽輪機工作狀態中高溫密封齒膨脹變形和直接剛度與阻尼的變化對轉子運動特性的影響。Cangioli等[19]建立了一種新的迷宮式密封的熱彈性體流模型,并分析了高溫下邊界壓力、轉速、進口預旋比等對轉子運力特性和穩定性的影響。在轉子-軸承系統的研究中為讓模擬更接近實際情況,部分學者提出單純線性理論的考察有一定局限性和考慮非線性因素是十分必要的并推導了包含非線性因素的動力學方程[20-22]。Elliott等[23]分析了非線性阻尼系統討論了計算方法以及實際例子。Hou等[24]利用多尺度方法歸納飛機旋翼系統非線性特征得到并求解運動方程分析了其運動特性。Iskakov等[25]利用線性阻尼與非線性阻尼綜合表達系統中的總阻尼,分析了其對系統穩定性的影響。Sayed等[26]嘗試利用二階微擾方確定軸承非線性剛度、阻尼系數并確定其有效范圍。Civera等[27]利用具有非線性剛度和阻尼項的解析SDoF模型,計算了高度靈活機翼的大振蕩,并討論了其在頻域中的穩態響應。Amer等[28]研究了具有三自由度的運動系統中非線性阻尼對其穩定性的影響。目前針對汽輪機轉子-軸承-密封系統的運動特性的研究基本上都基于線性理論分析,在轉子渦動過程中僅用線性剛度、阻尼來表達轉子的彈性恢復力,為更真實的反應轉子的運動狀態,轉子的非線性動態特性不能忽略。

為準確地考察超超臨界汽輪機轉子-軸承-密封系統的運動特性,本文充分考慮了系統中非線性動態特性以及熱、動載荷對其的影響,推導了不同條件下的轉子運動微分方程并通過試驗驗證了其準確性。在此基礎上建立轉子-軸承-密封系統Jeffcott模型,將高精度數值模擬求得的汽流激振力,建立密封激振力的擬合公式,將其耦合到轉子運動方程中。利用龍格-庫塔法求解對應的運動方程,分析不同條件下轉子的運動特性與穩定性。為抑制汽流激振、提高汽輪機的運行穩定性提供可靠理論依據。

1 數學模型

1.1 物理模型

根據某超超臨界汽輪機高壓缸第二壓力級設計參數,建立1.5級轉子-軸承-密封全周物理模型,如圖1所示,考察靜葉、葉頂、動葉運行過程中產生的激振力。模型具體結構參數如表1所示。工作條件如表2所示。采用曹麗華等的動網格計算方法建立渦動模型,將動靜間設為變形面,使密封與葉片流域相通,可準確反應出密封內部流動,有效的解決數值仿真中的網格變形問題。根據不同負荷下的蒸汽參數來設置計算時對應的邊界條件,可獲得8組不同偏心量,每組12個不同頻率的汽流激振力,將其按直角坐標系分解,其中x方向的激振力如圖2所示。

圖1 1.5級物理模型

表1 結構參數

表2 計算模型的工作條件

1.2 油膜力模型

利用Capone短軸承油膜力模型[29],獲得x、y方向的非線性油膜力。短軸承是指長徑比較小的軸承,采用此模型是由于其具有良好的計算精度與收斂性,通過改進長徑比可以獲得更加適用于超臨界汽輪機的軸承特性。

圖2 x方向不同偏心量下的激振力

(1)

其中：

1.3 激振力模型與運動微分方程

通過Fluent和動網格技術所獲得的不同偏心量激振力具有十分明顯的非線性特征難以用現有公式描述,將其進行多頻渦動計算,利用MATLAB冪次擬合多項式的方法表達激振力函數,如式(2),可準確的描述任何頻率與偏心量變化的非線性汽流激振力,其具體計算公式如式(3)。圖3為x方向激振力MATLAB擬合函數圖。

Fax=f(Ω,Cr)

Fay=f(Ω,Cr)

(2)

Fax=3 212+2 846Cr-1 384f-1 718Cr2-511Crf-
721.8f2+68.94Cr3+199.8Cr2f-934.1Crf2+
1 818f3+775.4Cr4+139.6Cr3f+52.13Cr2f2+
1 014Crf3+321.3f4-134.2Cr5-72.18Cr4f+
60Cr3f2-45.59Cr2f3+385.2Crf4-210.4f5

Fay=3 381+1 895Cr+85.31f-3 170Cr2+1 471Crf-
1 634f2+2 011Cr3-214.6Cr2f-1 667Crf2+
845.9f3+2 366Cr4-1 199Cr3f+368.3Cr2f2+
514.7Crf3+720.9f4-1 170Cr5+546.5Cr4f+
35.13Cr3f2-187.5Cr2f3+848.4Crf4-62.01f5

(3)

式中：Fax、Fay分別為x,y方向的非線性汽流激振力;Cr為偏心量;f為渦動頻率。

圖3 x方向不同偏心量下激振力擬合函數圖

圖4 Jeffcott模型示意圖

具有一般線性單自由度機械系統的運動微分方程為

(4)

表3 系統模型的主要參數

圖5 轉子轉動示意圖

式中：m為系統集中質量;c、k分別為系統線性阻尼、剛度系數。

將其擴展到轉子-軸承-密封系統可得

(5)

為使方程能更準確地反映轉子的運動狀態,分析中剛度、阻尼都用線性項和非線性項來綜合表達。首先只考慮系統非線性剛度式(4)轉化為

(6)

式中,γ為系統的非線性剛度系數,且與轉子位移的立方成正比。分析中將非線性剛度跟隨轉子位移分解成在x、y方向的分量。

此時系統的運動微分方程為

(7)

當只考慮系統非線性阻尼時,式(4)轉化為

(8)

此時非線性阻尼與轉子速度的n次冪成正比,為體現轉子移動速度的方向則有

(9)

(10)

同時考慮系統內非線性剛度、阻尼時,式(4)轉化為

(11)

運動微分方程為

(12)

1.4 試驗驗證

現有的研究已發現工程實際中轉子內的非線性動態特性對其運動特性有明顯影響但并沒有與激振力耦合研究。為驗證激振力模型及包含非線性動態特性的運動方程的準確性,基于如圖6的轉子-軸承-密封試驗臺進行驗證,其具體參數如表4所示。試驗臺通過調整轉子平臺高度來控制齒頂間隙,使其在0.5～3 mm范圍內可調,并通過多次試驗測量對其進行修正,確保誤差低于5%。通過調控風量來控制轉子轉速,使其穩定增加到316 rad/s,合理布置電渦流鍵相傳感器,采集轉子的位移信息,考察轉子不同轉速的渦動幅值。圖7、8分別為不同轉速的振幅曲線和額定轉速頻譜圖,展示了不考慮非線性動態特性,考慮非線性動態特性與試驗值的對比,其中不考慮的由式(5)模擬獲得,考慮的由式(12)模擬獲得。圖7中可以看出考慮非線性動態特性系數的幅值曲線更接近試驗值,特別是在汽輪機額定轉速附近更為明顯。圖8也可看出額定轉速下考慮非線性動態特性系數的頻譜圖與試驗值更為接近。說明本次研究模型設計合理、運動方程正確,也說明考慮系統內非線性動態特性能夠更準確的反應轉子的運動特性。

圖6 轉子-軸承-密封試驗臺

圖7 轉子振幅

2 計算結果分析

為深度揭示非線性動態特性系數對系統運動特性的影響,把非線性剛度、阻尼及其耦合效果分別進行分析。利用龍格-庫塔法求解,可得不同條件下轉子渦動狀態。計算過程中每完成一步,流場隨之更新,通過不斷迭代計算,考察不同條件下系統的運動特性。

圖8 額定轉速頻譜圖

表4 試驗臺的主要參數

2.1 非線性剛度對系統運動特性的影響

圖9為不同非線性剛度轉子系統分岔圖,直觀來講隨著負荷的增加轉子的運動都經歷了類似混沌二周期-類似混沌一周期-類似混沌二周期-具有分散性的混沌周期四種變化。開始的類似混沌二周期在1%～5%THA(turbine heat acceptance)左右,這是由于系統初帶負荷,突然的載荷使轉子發生擾動,之后轉子位移變小進入類似混沌一周期,直到33%THA左右混沌一周期開始向混沌二周期轉變,經過短暫的過渡,分岔圖主要集中在上下兩部分。負荷達到58%THA后轉子進入復雜的混沌運動,當系統進入VWO(valve whole open)工況后有向混沌三周期轉化趨勢。觀察非線性剛度介入后的變化,非線性剛度對負荷較低的情況影響不大,系統進入混沌二周期時,當γ=3×109N/m3、γ=3×1010N/m3時分岔過程稍加柔和,特別是后者分岔過程柔和且分布集中,而γ=3×1011N/m3時此區域變得零散,混沌特征明顯。對比負荷高于55%THA后的情況,考慮非線性剛度后轉子運動更加集中、位移減小,其中γ=3×1010N/m3時最為明顯。當系統達到VWO工況后非線性剛度會使轉子運動轉向混沌三周期趨勢更加明顯?？傮w來看一定范圍內的γ值可減小轉子位移量。

圖9 不同γ值分岔圖

將轉子的無量綱位移進行快速傅里葉變換,可得到不同非線性剛度轉子系統頻率瀑布圖,如圖10所示。展示了激振力與油膜力對系統升負荷過程中造成的分頻現象?？傮w上看在6%THA之前系統都有幅值較小的1/2工頻,隨著負荷的增加逐漸消失,這是由于系統初帶負荷產生的。28%THA左右1/2工頻又再次出現且幅值較高,60%THA后,1/2工頻逐步分化成1/3、2/3工頻。在分化出現后密頻效果愈發明顯,這是由于負荷越高,密封腔流場越復雜,汽流激振力非線性特征越明顯,多頻渦動中的其它頻率振動開始展現。對比發現耦合非線性剛度后密頻現象更為明顯,尤其在60%～70%THA時。考察四種情況整個分化部分密頻現象的增加情況,其中γ=3×1010N/m3時密頻現象較少。

圖10 不同γ值瀑布圖

圖11為不同非線性剛度的主要頻率無量綱振幅圖,展示了系統的主要頻率分布。本次研究中不僅考察汽流激振力對轉子系統的影響還包含油膜力,隨著系統負荷的增高,激振力逐步加大與油膜力產生耦合作用,引起以1/2工頻為主的振動,隨著激振力進一步增大,1/2工頻開始向1/3、2/3工頻演化,其中1/3工頻主要由激振力作用產生,2/3工頻主要由油膜力作用產生,而且由于高負荷時汽流激振力更大作用更明顯,因此1/3工頻幅值更高。如圖考慮非線性剛度后對工頻基本沒有影響,但對其它主要頻率有一定影響。隨著γ取值的增高1/2工頻位置稍向高負荷方向移動,且區域略有增加,1/3、2/3工頻區域幅值略有降低。這是說明非線性剛度延緩了激振力與油膜力的作用效果。且在主要頻率高幅值集中區都有低幅值出現,這也對應了圖10耦合非線性剛度后密頻現象明顯。說明非線性剛度減小了激振力與油膜力的作用效果,但導致多頻渦動中的其它頻率振動增多。

2.2 非線性阻尼對系統運動特性的影響

圖12為不同非線性阻尼轉子系統分岔圖,非線性阻尼對轉子位移的整體運動趨勢沒有改變,但使系統由混沌一周期進入混沌二周期的分岔過程變得柔和,在混沌二周期階段非線性阻尼介入后此區域范圍明顯變大,并隨著ε值的增加逐步減小。系統進入復雜混沌運動后ε值越高轉子位移越小。當系統達到VWO工況后轉子轉向類似混沌三周期趨勢減弱。

圖13為不同非線性阻尼轉子系統頻率瀑布圖,圖中四種情況在負荷達到30%THA后都出現了相似的密頻現象,且考慮非線性阻尼后更明顯。在系統由1/2工頻向1/3、2/3工頻演化初期集中的密頻現象隨著ε值的升高而減弱,且整個演化過程的密頻振幅也隨ε值的升高減小。

圖14為不同非線性阻尼的主要頻率無量綱振幅圖,由圖14(b)可知,1/2工頻隨著ε值的增加有向高負荷方向的移動,且幅值略有降低,區域有所增加。說明在負荷較低時非線性阻尼的增加使系統產生1/2工頻需要更大的耦合力(激振力與油膜力),且需要更大的激振力開始向1/3、2/3工頻演化。由圖14(b)、(c)可知,隨著ε值的增加1/3、2/3工頻區域明顯減少,且幅值降低。

(a) 工頻

(b) 1/2工頻

(c) 1/3工頻

(d) 2/3工頻

圖12 不同ε值分岔圖

圖13 不同ε值瀑布圖

(a) 工頻

(b) 1/2工頻

(c) 1/3工頻

(d) 2/3工頻

2.3 耦合熱、動載荷后系統的運動特性

為了解非線性剛度、阻尼的耦合作用,使研究更符合轉子的實際工作情況,分析了耦合熱、動載荷前后系統的運動特性。耦合熱、動載荷后系統設計參數有一定變化如密封間隙、轉子直徑等。不僅使激振力發生改變同時改變了系統的非線性剛度、阻尼,如表5所示。研究中耦合熱、動載荷汽流激振力獲取流程如下：設計參數下密封流場計算,獲得流場溫度、壓力分布;分析耦合后密封結構,獲得耦合熱、動載荷后密封腔尺寸變化;計算密封腔變形后的渦動流場,獲得耦合后的激振力。

圖15為耦合熱、動載荷前后的分岔圖。相比之下耦合后系統更早的進入了混沌二周期階段,且此階段轉子位移略有減小。當轉子進入復雜混沌周期后,耦合前系統轉子位移更小尤其在60%～75%THA階段。在高負荷區域耦合熱、動載荷后的系統基本沒有明顯的運動趨勢,但轉子位移有所減小,說明耦合后系統流場更復雜難以形成明顯的周期性運動,且汽流激振力更大,對轉子位移有一定限制。

表5 耦合熱、動載荷參數

圖15 耦合熱、動載荷前后分岔圖

圖16為耦合熱、動載荷前后轉子系統頻率瀑布圖,如圖所示耦合前后系統均有明顯密頻現象,其中耦合后較高的密頻區域主要集中在1/2工頻向1/3、2/3工頻轉化初期與1/3工頻前端區域,其它區域較為緩和,尤其是高負荷區域。這是因為耦合熱、動載荷后密封間隙更小、密封內流場更復雜,汽流激振力較大非線性更強,且系統非線性剛度減小,在激振力作用為主的區域多頻渦動體現明顯。

圖17為耦合前后主要頻率無量綱振幅圖,其中耦合后的1/2工頻區域向低負荷方向有所移動且區域更小,這是由于考慮熱、動載荷后激振力更大且增長更快,與油膜力耦合后使系統更早的出現1/2工頻,更早的向其它工頻演化。1/3、2/3工頻耦合后的幅值有所降低,區域變化不大。

圖16 耦合熱、動載荷前后瀑布圖

(a) 工頻

(b) 1/2工頻

(c) 1/3工頻

(d) 2/3工頻

3 穩定性分析

常用最大Lyapunov指數來衡量轉子動力系統的穩定性,本次研究采用wolf法對其值進行計算。其值大于0代表系統處于混沌的、難以預測的隨機運動;其值小于0則代表系統處于周期運動,屬于可預測的有一定規律可尋的運動狀態?；煦绫硎痉蔷€性系統在某種條件下呈現出難以預測的隨機現象,是有序與無序融為一體的現象。

圖18為不同非線性剛度最大Lyapunov指數曲線與指數平均值。如圖所示雖然指數值很小低于10-3但大部分區域高于0,說明轉子基本處于混沌狀態。比較四種情況,在負荷低于30%THA區域指數變化不大,對比分叉圖此時系統處于類似混沌一周期運動,汽流激振力較小非線性剛度的介入作用不明顯。隨著負荷的增加轉子進入混沌二周期運動,此區域指數值均有所上升這說明混沌二周期運動更為復雜,系統穩定性波動明顯?？紤]非線性剛度后指數值有所下降但小于0的區域沒有了,γ=3×109N/m3與γ=3×1011N/m3在57%THA處還出現了高峰值,說明非線性剛度對此區域轉子的運動狀態影響明顯。當轉子進入復雜混沌運動后,指數值有降低趨勢還出現了小于0的穩定區域,說明此時系統運動的失穩特征減弱。而考慮非線性剛度后此區域指數值出現更多小于0區域,其中γ=3×1010N/m3時比較明顯,說明合理的非線性剛度能夠改善系統高負荷運行的穩定性。觀察指數的平均值,考察非線性剛度后平均值有所上升,γ=3×1010N/m3時最低。

(a) 最大Lyapunov指數

(b) 指數平均值

圖19為不同非線性阻尼的最大Lyapunov指數曲線與指數平均值。由圖19(a)可知,考察非線性阻尼后指數值有所上升且波動變大,但也出現了更多指數小于0的區域,波動也隨著ε值的增加有所改善。當系統負荷達到35%THA時ε=10 000 N·s2/m2、ε=100 000 N·s2/m2指數出現較高峰值,說明這兩種情況轉子在分岔過渡過程中穩定性較差。當達到67%THA時ε=10 000 N·s2/m2、ε=50 000 N·s2/m2出現高峰值。而后指數值有所下降并出現小于0區域。由圖19(b)可知,考慮非線性阻尼后指數平均值明顯上升,隨著其值的升高逐步降低?？傮w來講隨著非線性阻尼值的增加系統的穩定性有所改善。

(a) 最大Lyapunov指數

(b) 指數平均值

圖20考察了耦合熱、動載荷前后的指數曲線與平均值。相比之下耦合前后低負荷區域指數值區別不大,系統進入混沌二周期運動后直至73%THA耦合后的指數值波動更明顯且指數值更高并在50%THA處出現高值,對比分岔圖耦合后的混沌二周期區域更長,說明混沌二周期蘊含的失穩因素更多更復雜。當系統負荷超過75%THA耦合后的指數值明顯更低,且小于0的區域更多,說明耦合后的系統高負荷區域更穩定。這是由于耦合后非線性剛度變小阻尼變大,且密封腔內間隙更小,產生更大的汽流激振力對轉子的渦動幅度產生了限制。

4 結 論

為更準確地反應超超臨界汽輪機轉子-軸承-密封系統的運動特性,本文改善了轉子運動微分方程,并通過試驗對比驗證了方程的準確性。利用龍格-庫塔法對其求解,討論了不同非線性動態特性系數以及耦合熱、動載荷后系統的運動特性與穩定性。具體結論如下：

(1) 通過對比不同條件下系統最大Lyapunov指數,系統均有失穩可能,且考慮非線性動態特性后其均值有所上升;其中一定范圍內的非線性剛度能夠改善系統的穩定性;而非線性阻尼值越高系統穩定性越好;耦合熱、動后的系統高負荷區域更穩定。

(2) 系統中非線性動態特性系數會使轉子不同類型的運動區域發生改變,非線性剛度在一定范圍內可以減小轉子的位移;而非線性阻尼會使混沌二周期區域增加,減小較高負荷區域轉子的渦動幅值。對比耦合熱、動載荷前后,高負荷區域耦合后的轉子位移變小。

(a) 最大Lyapunov指數

(b) 指數平均值

(3) 考察不同條件下系統頻率分布,存在1/2、1/3、2/3工頻,且隨著負荷的增加密頻現象愈發明顯。其中考慮非線性剛度后密頻效果略有增加,但合適的非線性剛度值可以緩解此現象?？紤]非線性阻尼后系統由1/2工頻向1/3、2/3工頻的分岔及演化過程中的高幅值密頻隨著其值的升高而減弱;耦合熱、動載荷后明顯的密頻現象主要集中在激振力作用為主的工頻區域。

(4) 系統中的非線性剛度、阻尼均會使1/2工頻區域向高負荷方向移動且范圍略有增加,1/3、2/3工頻幅值降低,且非線性阻尼會縮小1/3、2/3工頻區域;耦合熱、動載荷后系統的1/2工頻出現更早,更早的向其它工頻演化。