王 建， 宋志強， 劉昱杰

(西安理工大學 水利水電學院， 陜西 西安 710048)

1 研究背景

水電在國家經(jīng)濟、能源發(fā)展及環(huán)境安全等方面有著巨大的影響，隨著水輪發(fā)電機容量和機組尺寸的增大，水電站機組及廠房的振動問題成為制約機組高效運行，甚至威脅其安全的關鍵。引發(fā)水電站機組及廠房振動的因素除了機組和廠房結構本身的非線性特性外，水力、機械、電磁振源的復雜性及其相互作用也是一個重要的方面。由于過流部件構造復雜，水輪機運行工況多變使得獲知、表征和施加水力振源成為一大難題。轉(zhuǎn)輪受到水流作用產(chǎn)生一種近似隨機的變值荷載，在以往的軸系統(tǒng)分析中，只能將其忽略或簡化為只考慮一種優(yōu)勢頻率的簡諧荷載，該荷載的無法準確模擬給軸系統(tǒng)振動分析造成了很大影響[1-3]。

關于水輪發(fā)電機組水力振源的研究目前已經(jīng)取得了較大進展，文獻[4]應用Newmark算法計算了不同工況條件下水輪機轉(zhuǎn)輪部位應力，考慮了周期性的水力振源作用下結構的振動分析，但并未考慮脈動壓力的非周期性部分；文獻[5]運用轉(zhuǎn)子動力學軟件ARMD給予了機組軸系橫向振動的評價，并準確預估了不同工況下水力激勵力對機組瞬態(tài)響應的振動分析，提出了水力激勵力對軸系穩(wěn)定性研究的重要意義；文獻[6]建立ANSYS主軸模型進行了軸系結構自振特性分析，分析了多種外激勵共同作用下時振動響應。隨著計算流體動力學(CFD)的應用及發(fā)展，合理獲知機組過流部件表面脈動壓力成為了可能[7]，在流體計算結果越來越準確的情況下，將CFD獲知的脈動壓力導入軸系結構進行水力振源的模擬，有了更為可靠的理論依據(jù)和試驗參考，很多學者已經(jīng)在CFD計算的基礎上對流道內(nèi)的復雜流體擾動進行了探究并結合模型試驗驗證了CFD流體計算的可行性[8-12]，在此基礎上，Khare等[13]得出了模型試驗所不能觀測到的細節(jié)部位的流體運行情況，并分別計算了不同導葉開度和轉(zhuǎn)輪轉(zhuǎn)速條件下的流體運行特性；尤建鋒等[14]建立了與試驗模型一致的完整CFD流體計算模型并進行了多工況的數(shù)值模擬，對比分析了壓力脈動特性；錢忠東等[15]通過改變導葉開度對比分析了不同工況條件下流道內(nèi)各部位脈動壓力分布規(guī)律。上述研究雖然準確獲知了流道脈動壓力，但更多是研究機組水力效率、轉(zhuǎn)輪葉片疲勞破壞及空化等，合理施加轉(zhuǎn)輪脈動壓力進行機組整個軸系統(tǒng)結構振動分析的研究目前尚不多見。

機組在受到轉(zhuǎn)輪處脈動壓力作用的同時也受到定轉(zhuǎn)子氣隙不均勻產(chǎn)生的不平衡電磁力和轉(zhuǎn)子-轉(zhuǎn)輪偏心產(chǎn)生的離心力影響，本文通過CFD計算提取轉(zhuǎn)輪表面脈動壓力時程，轉(zhuǎn)化為荷載文件，實現(xiàn)轉(zhuǎn)輪脈動壓力更為合理的模擬和施加，在此基礎上綜合電磁、機械振源進行了多振源耦合情況下水電機組三維轉(zhuǎn)子-軸承系統(tǒng)結構的振動分析，并對轉(zhuǎn)子-轉(zhuǎn)輪不平衡力相位差等進行了敏感性分析。

2 脈動壓力的CFD計算及提取

2.1 全流道流體域計算模型

本文以某電站混流式水輪機為研究對象，其轉(zhuǎn)輪直徑2.165 m，轉(zhuǎn)輪葉片13個，固定導葉12個，活動導葉24個，導葉開度角23°，模擬工況條件為最高水頭、最優(yōu)負荷(H=77 m，P=102 MW)和額定水頭、額定負荷(H=61 m，P=88 MW)下水輪機的內(nèi)部流動特性，采用Unigraphics NX 8.0軟件建立包括蝸殼、固定導葉、活動導葉、轉(zhuǎn)輪和尾水管在內(nèi)的全流道計算域模型。由于水輪機部件結構復雜，各部件之間相互關聯(lián)，為確保模型的完整性，減小分析誤差，采用ANSYS ICEM-CFD軟件進行高精度的六面體網(wǎng)格劃分，流道模型共有4 842 301個單元，4 880 919個節(jié)點，流道內(nèi)各部件網(wǎng)格劃分如圖1所示。

圖1 流體計算域網(wǎng)格劃分示意圖

2.2 控制方程

對三維不可壓縮非穩(wěn)態(tài)Navier-Stokes方程采用時間平均法，進行Reynolds時均法得到連續(xù)方程和動量方程分別為式(1)和(2)：

(1)

(2)

湍流模型采用剪切應力輸運SST(Shear Stress Transport)模型，用以模擬旋轉(zhuǎn)部件與靜止部件的動靜干涉，固體壁面采用標準壁面函數(shù)進行處理，控制方程為式(3)：

(3)

(4)

式(1)～(4)中：ρ為流體的密度，kg/m3；μ為動力黏性系數(shù)；u為表征速度，m/s；τ為雷諾應力張量，MPa；Pk為湍流脈動動能k的生成項，k由k-ε模型求得；Pω為湍流脈動頻率ω的生成項，由k-ω求得；Γk，Γω為k與ω的有效擴散系數(shù)；Yk，Yω為k與ω的耗散項；Dω為正交擴散項。

2.3 邊界條件

計算域的進口為金屬蝸殼的進口，流體垂直進口流入，進口邊界條件設置為質(zhì)量-流量邊界；出口邊界條件為壓力出口，指定平均靜壓為0；固體壁面為無滑移、絕熱壁面；交界面方式采用動靜轉(zhuǎn)子交界。

2.4 控制方程的離散和求解設置

模型采用有限體積法進行離散，建立多種運行條件下的原型全流道模型，利用大型有限元軟件ANSYS中的CFX模塊，采用Realizable 模型和SIMPLEC算法對水輪機過流部件的內(nèi)部流動進行計算，從振源特性和水流特性分別進行流體計算。轉(zhuǎn)輪轉(zhuǎn)速為136.4 rad/min，轉(zhuǎn)頻為2.273 Hz，旋轉(zhuǎn)周期為0.44 s，轉(zhuǎn)輪瞬態(tài)分析中，將周期分為120個時間步進行脈動壓力提取，時間步長取為Δt=0.00367 s，其中，轉(zhuǎn)子-轉(zhuǎn)輪軸系結構旋轉(zhuǎn)角速度為ω=2πn=14.2838 rad/s，定常計算結果提取轉(zhuǎn)輪旋轉(zhuǎn)兩個周期的脈動壓力時程值。

2.5 脈動壓力結果分析

為了直觀地分析引水部件和轉(zhuǎn)輪內(nèi)部的壓力脈動情況，分別在水輪機蝸殼進口、轉(zhuǎn)輪進口、葉道間，轉(zhuǎn)輪出口布置了4個壓力觀測點如圖2所示。通過完成蝸殼進口脈動壓力的提取并進行傅里葉變換后與模型試驗實測數(shù)據(jù)對比如圖3所示?？梢奀FD脈動壓力計算結果可信，并在此基礎上進行轉(zhuǎn)輪部位的脈動壓力提取。

水輪機轉(zhuǎn)輪是整個水輪發(fā)電機組的核心部件，轉(zhuǎn)輪中流態(tài)的好壞直接影響了機組的性能。對兩種工況轉(zhuǎn)輪處脈動壓力提取，處理流體計算結果可知最優(yōu)工況和額定工況下固定導葉和活動導葉通道內(nèi)均無脫流，流線比較光順，固定導葉出口角和活動導葉進口角配合很好，提取的脈動壓力結果可以作為本文研究數(shù)據(jù)參考?，F(xiàn)以最高水頭H=77 m、最優(yōu)工況條件下的轉(zhuǎn)輪脈動壓力提取結果為例，就轉(zhuǎn)輪脈動壓力的時域特性和頻域特性進行簡要分析。

脈動壓力的近似隨機性通過提取轉(zhuǎn)輪進口、葉道間、轉(zhuǎn)輪出口3個監(jiān)測點脈動壓力時域曲線進行表征，并就傅里葉變換后的頻域曲線進行轉(zhuǎn)輪部位脈動壓力分析如圖4所示。

由圖4中壓力脈動時域圖(圖4(a)、4(c)、4(e))和頻域圖(圖4(b)、4(d)、4(f))可以看出，轉(zhuǎn)輪進口脈動壓力幅值較大，轉(zhuǎn)輪進口壓力幅值為0.423 MPa，葉道間脈動壓力幅值為0.125 MPa，轉(zhuǎn)輪出口脈動壓力幅值為0.079 MPa；轉(zhuǎn)輪進口和葉道間主頻均為2.23 Hz的低頻，該頻率是轉(zhuǎn)輪轉(zhuǎn)頻2.273 Hz的約1倍，在整個水輪機流道內(nèi)傳播，此外，轉(zhuǎn)輪進口處還存在53.4 Hz的高頻，是轉(zhuǎn)輪葉片數(shù)(Zr=13)的4倍，該高頻脈動是由于轉(zhuǎn)輪和導葉的動靜干涉作用產(chǎn)生的，轉(zhuǎn)輪出口主頻為接近1倍轉(zhuǎn)頻的低頻，由各部位頻域圖對比可以看到葉道間振幅增加明顯。

3 軸系統(tǒng)模型建立及荷載施加

3.1 軸系統(tǒng)三維有限元模型建立

利用有限元軟件Ansys建立軸系統(tǒng)三維模型，模型由上導軸承、發(fā)電機轉(zhuǎn)子、下導軸承、水導軸承、法蘭、主軸及轉(zhuǎn)輪組成。轉(zhuǎn)輪采用solid187單元離散，主軸結構采用solid45單元離散。用combin14彈簧單元模擬上導、下導、水導軸承，軸系統(tǒng)單元總數(shù)為64 936個，節(jié)點數(shù)為113 904個，實際單元劃分及結構模型如圖5所示。

圖2 流體計算壓力觀測點位置示意圖 圖3 蝸殼進口處脈動壓力頻域圖結果對比

圖4 轉(zhuǎn)輪特征部位脈動壓力時域圖和頻域圖

導軸承的剛度和阻尼隨著軸頸運動及導軸承間隙變化而非線性變化，而且剛度與阻尼之間存在一定的相互關系，即當導軸承的剛度發(fā)生變化時阻尼也隨之相應變化[1]，如圖6所示，其中λk=Kx/μnl(R/C)3，λc=ωCx/μnl(R/C)3,μ為軸承內(nèi)油的黏度，N·s/cm2；n為機組額定轉(zhuǎn)速，rad/s；l為軸瓦的軸向長度，m；R為導軸承軸瓦半徑，m；C為導軸承間隙，m；C′為軸瓦在支承點處的間隙，m；ω為大軸旋轉(zhuǎn)圓頻率，rad/s；Kx為導軸承剛度，N/m。

以某實際水電站機組為例，取導軸承軸瓦半徑R=0.77 m，轉(zhuǎn)速n=136.4 r/min，ω=14.2838 rad/s，導軸承間隙C=0.02 cm，黏度μ=0.33×10-5N·s/cm2。代入上述數(shù)據(jù)后，Kx=1.709×106λk，Cx=0.12×106λc。

擬取導軸承剛度為Kx=2×108N/m，λk=1.171，由圖查得C′/C=0.79，得λc=2.25，導軸承阻尼系數(shù)為Cx=2.69×107N/m，為使計算簡便，本文計算模型擬取3個導軸承剛度系數(shù)和阻尼系數(shù)均為Kx和Cx。

3.2 荷載計算及施加

3.2.1 水力振源 作用在轉(zhuǎn)輪上的水力分為周期性和非周期性兩種，分別對應水力不平衡力和脈動壓力，周期性力即水力不平衡力，水力不平衡力主要是由轉(zhuǎn)輪部分不圓或機組大軸擺度使迷宮間隙發(fā)生周期性變化產(chǎn)生的迷宮壓力脈動，取值參照式(5)計算[16]。

F1=K水ε

(5)

式中：K水為水的不平衡系數(shù)；ε為轉(zhuǎn)輪處擺度，m；算例取F1=39.12 kN。

圖6 導軸承剛度與阻尼系數(shù)的關系

近似隨機的非周期性水力振源一般難以用傳統(tǒng)方法考慮。非周期性脈動壓力即上述CFD計算所得轉(zhuǎn)輪表面實時脈動壓力。模擬過程是在Ansys軟件中建立以坐標原點為基準點的轉(zhuǎn)輪模型，進行網(wǎng)格劃分，提取轉(zhuǎn)輪區(qū)域各節(jié)點編號及其坐標，寫入節(jié)點文件，同時建立以坐標原點為基準點的全流道計算模型，保證結構體與流道模型的位置對應，在此基礎上應用CFD-ICEM進行流體域網(wǎng)格的劃分，完成CFD流體非定常計算后，進入CFX-Post后處理中編寫節(jié)點導入命令，提取各節(jié)點部位的脈動壓力時程后，利用APDL語言將脈動壓力施加到結構中，至此完成脈動壓力的合理模擬及施加。

3.2.2 機械不平衡力 假設只存在動偏心，各部位軸心在初始位置都是原始位置零點，系統(tǒng)偏心是由于系統(tǒng)轉(zhuǎn)動所引起的旋轉(zhuǎn)部件質(zhì)量分布不均勻引起的，所產(chǎn)生的離心力為：

F2=me0ω2

(6)

式中：m為轉(zhuǎn)動部件的質(zhì)量，kg；ω為轉(zhuǎn)動角頻率，rad/s。

離心力主要產(chǎn)生于轉(zhuǎn)子和轉(zhuǎn)輪處，根據(jù)已知參數(shù)轉(zhuǎn)子質(zhì)量為mr=6.5×106kg，轉(zhuǎn)輪質(zhì)量為mt=3.5×106kg，求得F2r=856.913 kN、F2t=583.51 kN ，不平衡力作為外力施加在相應的結構上。偏心距參考水輪機制造標準：轉(zhuǎn)子偏心er=6.3/ω，轉(zhuǎn)輪偏心et=40/ω，取er=0.44 mm，et=1.2 mm。

3.2.3 不平衡電磁拉力 由于制造、安裝偏差及運行過程中的磨損、變形等會造成定子和轉(zhuǎn)子間產(chǎn)生不均勻間隙，進而產(chǎn)生不平衡電磁拉力，當轉(zhuǎn)子偏心為e時，不平衡磁拉力為[17]：

(7)

式中:L、D分別為轉(zhuǎn)子高度和直徑，m；取L=1 134 cm；D=275 cm;B為磁通量密度，取2393×10-4T；β為系數(shù)，一般取值0.3左右；δ為正常間隙，取值為1.6 cm，e=0.44 cm。本文不平衡電磁拉力取F3=152.7 kN。

為了保證數(shù)值模擬更接近機組實際運行條件，不平衡電磁拉力和離心力均由互相垂直的正弦荷載和余弦荷載來模擬。

4 基于水力振源合理施加的軸系振動分析

4.1 水力振源對軸系振動的影響

4.1.1 脈動壓力對軸系結構的影響分析 將CFD提取的脈動壓力結果進行處理后單獨施加到軸系結構中進行動力響應分析，并提取轉(zhuǎn)輪、轉(zhuǎn)子周期內(nèi)的軸心偏移情況如表1，并根據(jù)XY向(水平向)位移進行偏移合成繪制轉(zhuǎn)輪、轉(zhuǎn)子軸心偏移圖如圖7。

結合表1和圖7可以看出，在兩種工況條件下，導入脈動壓力，轉(zhuǎn)輪、轉(zhuǎn)子軸心運行均受到影響。當水頭H=61 m時，轉(zhuǎn)輪軸心偏移為0.005～0.238 mm，轉(zhuǎn)子軸心偏移為0.001～0.01 mm；當水頭H=77 m時，轉(zhuǎn)輪軸心偏移為0.009～0.233 mm，轉(zhuǎn)子軸心偏移為0.001～0.007 mm。水頭H=61 m時轉(zhuǎn)輪、轉(zhuǎn)子位移偏移較大，原因是由于此工況條件下轉(zhuǎn)輪部位脈動壓力振幅最大，轉(zhuǎn)輪處于部分負荷狀態(tài)，振動更加明顯，轉(zhuǎn)輪部位脈動壓力較大，對整個軸系結構有明顯的影響。

表1 兩種水頭工況下轉(zhuǎn)輪、轉(zhuǎn)子偏移情況

4.1.2 振源耦合條件下軸系統(tǒng)振動分析 通過在轉(zhuǎn)子處施加不平衡電磁拉力模擬電磁振源分布、加載離心力的方式模擬轉(zhuǎn)子和轉(zhuǎn)輪處的機械振源以及導入實時脈動壓力時程模擬水力振源，研究整個轉(zhuǎn)輪-轉(zhuǎn)子軸系結構在機械振源、電磁振源、水力振源3種振源耦合情況下的振動響應。圖8為軸系結構在只導入機械振源和電磁振源作用下的轉(zhuǎn)輪、轉(zhuǎn)子軸心軌跡圖，圖9為兩種運行工況下，軸系結構在3種振源耦合情況下的振動分析，并根據(jù)響應結果繪制的軸系統(tǒng)關鍵部位軸心軌跡圖。

由圖8、9轉(zhuǎn)輪和轉(zhuǎn)子軸心軌跡圖對比可知，在脈動壓力的非線性影響下，結構各部位軸心運行軌跡呈現(xiàn)不規(guī)則變化，從圖9軸心軌跡圖可以看出轉(zhuǎn)輪、下導和水導受脈動壓力影響較為劇烈，當水頭H=61 m時，導入脈動壓力后轉(zhuǎn)子擺動幅值為0.102～0.107 mm，轉(zhuǎn)輪擺動幅值為0.586～0.725 mm，轉(zhuǎn)子擺動幅值增加了4.91%，轉(zhuǎn)輪擺動幅值增加了23.7%，上導軸承、下導軸承和水導軸承的擺動幅值分別為0.072、0.059和0.063 mm；當水頭H=77 m時，轉(zhuǎn)子擺動幅值為0.102～0.103 mm，轉(zhuǎn)輪擺動幅值為0.586～0.664 mm，轉(zhuǎn)子擺動幅值增加了0.98%，轉(zhuǎn)輪擺動幅值增加了13.2%，上導軸承、下導軸承和水導軸承的擺動幅值分別為0.071、0.058和0.062 mm。軸系結構各部位偏移情況對比可知H=61 m水頭作用下軸系結構振動較H=77m水頭更為劇烈。

圖7 兩種水頭工況下轉(zhuǎn)輪、轉(zhuǎn)子受脈動壓力影響的軸心偏移圖

圖8 電磁振源和機械振源耦合作用下轉(zhuǎn)輪、轉(zhuǎn)子軸心軌跡圖

圖9 3種振源耦合作用下軸系結構關鍵部位軸心軌跡圖

綜上可得，脈動壓力對整個結構的影響主要體現(xiàn)在轉(zhuǎn)輪上，轉(zhuǎn)輪軸心運行受脈動壓力非線性影響較為劇烈，轉(zhuǎn)子由于距離脈動壓力施加部位較遠，受到的影響較??；在額定水頭，額定工況條件下，脈動壓力對轉(zhuǎn)輪振動影響高于最大水頭，對于最優(yōu)工況條件下的轉(zhuǎn)輪振動，其脈動壓力對轉(zhuǎn)輪結構的振動有著重要影響，為使得分析更為符合實際，脈動壓力的合理模擬和施加尤為重要。

4.2 轉(zhuǎn)子與轉(zhuǎn)輪荷載相位差對系統(tǒng)振動的影響

圖10 不同相位差條件下各部位軸心偏移對比圖

由圖10(a)、(b)可以看出不考慮轉(zhuǎn)輪脈動壓力影響下，轉(zhuǎn)輪、轉(zhuǎn)子在電磁振源和機械振源耦合作用下振動平穩(wěn)，轉(zhuǎn)輪、轉(zhuǎn)子位移幅值均出現(xiàn)在θ=180°時，分別為0.707和0.116 mm，可見當θ=180°時為荷載的最不利組合；由圖10(c)可以看出，在考慮脈動壓力的振源耦合作用下，在最不利荷載θ=180°時，轉(zhuǎn)輪偏移幅值為0.757 mm，且當結構在振源作用0.07 s之后，轉(zhuǎn)輪開始劇烈振動，當θ=90°時，轉(zhuǎn)輪軸心運行軌跡變動幅度最大，振動幅值為0.751 mm，略小于θ=180°時的位移幅值，θ=0°和θ=45°時轉(zhuǎn)輪運行軌跡變動范圍較小，且θ=0°時，轉(zhuǎn)輪運行軌跡范圍最小，但仍在0.41～0.66 mm之間擺動；由圖10(d)可以看出，轉(zhuǎn)子受相位角θ影響較小，但θ=90°時，轉(zhuǎn)子軸心運行軌跡變幅大于其他3種相位差，運行范圍在0.031～0.117 mm之間，可見相位角θ大小對轉(zhuǎn)輪運行穩(wěn)定有著重要影響，對轉(zhuǎn)子影響相對較小。 此外由各圖對比可知，在θ=180°的時候，轉(zhuǎn)輪位移幅值由0.707 mm變?yōu)?.757 mm，增長幅值為7.08%，轉(zhuǎn)子位移幅值由0.116 mm變?yōu)?.117 mm，變化相對較小，原因是因為轉(zhuǎn)輪受到脈動壓力直接作用變化較明顯，而轉(zhuǎn)子距離轉(zhuǎn)輪位置較遠，受脈動壓力影響較小。

5 結 論

(1) 建立水輪機全流道湍流模型，基于CFD數(shù)值模擬獲取轉(zhuǎn)輪葉片表面脈動壓力時程，通過Ansys命令流將脈動壓力導入軸系統(tǒng)結構的振動分析中，可實現(xiàn)水力脈動振源的合理模擬和施加，在此基礎上，綜合電磁和機械振源進行了多振源耦合作用下軸系統(tǒng)的振動分析。

(2) 軸系統(tǒng)考慮非周期水力即脈動壓力作用后，轉(zhuǎn)子振動變化較小，轉(zhuǎn)輪軸心偏移幅值較大。由于各點脈動壓力的方向是垂直該點所在表面的，轉(zhuǎn)輪葉片沿圓周循環(huán)對稱且本身具有的空間扭曲性，所以非周期性脈動壓力對轉(zhuǎn)輪某一個方向上的振幅影響不會很大，但使轉(zhuǎn)輪運行軌跡趨于復雜，振動趨于劇烈。軸系統(tǒng)在非周期性脈動壓力作用下的非線性變化明顯，合理模擬和施加轉(zhuǎn)輪葉片表面脈動壓力是必要的。

(3)轉(zhuǎn)子和轉(zhuǎn)輪荷載相位差為90°時，轉(zhuǎn)子受非周期性脈動壓力影響較大，相位差為180°時，轉(zhuǎn)輪受非周期性脈動壓力影響較大。導軸承支撐條件和荷載相位差對軸系統(tǒng)振動有著重要影響，在軸系統(tǒng)動力特性分析及評價中值得重視。