水泵水輪機啟動過程轉輪動應力分析

桂中華 夏 翔 王 薇 周凌九 劉殿海 王正偉

(1.國網新源控股有限公司抽水蓄能技術經濟研究院， 北京 100761; 2.中國農業大學水利與土木工程學院， 北京 100083;3.清華大學能源與動力工程系， 北京 100084)

0 引言

在運行過程(特別是過渡過程)中，水輪機轉輪承受劇烈的水壓力脈動，可能產生嚴重振動和高水平動應力，甚至可能產生共振和疲勞[1-3]。與常規水輪機組相比，水泵水輪機長期處于啟動、停機和負荷變化等過渡過程中，面臨的穩定性問題更加突出。準確預測水泵水輪機轉輪在過渡過程(如啟動過程)中的動應力水平具有重要意義。

試驗測量是一種有效分析轉輪動應力的手段[3-4]。但是，試驗方法的成本和難度很高，并且無法應用在設計階段。隨著計算機技術的進步，數值模擬為轉輪動應力的研究提供了便利。水輪機轉輪在水壓力作用下產生動力響應是典型的流固耦合(Fluid-structure interaction, FSI)問題。而流固耦合問題的求解方法可分為強耦合法和弱耦合法。其中弱耦合(即分離式求解)方法的實施相對容易，因此得到了廣泛的應用。一些學者采用分離式求解的雙向流固耦合方法對水輪機和水泵的內部流動及結構響應進行了研究[5-9]。上述工作都對研究對象作了一定簡化，因為對于復雜結構而言，雙向耦合方法需要消耗大量的計算資源，在實際工程中并不適用。事實上，在水力機械正常運行時，轉輪的彈性變形量遠小于流體流動的特征長度，忽略結構振動對流體流動的影響可滿足精度要求。因此，也有學者采用分離式求解的單向流固耦合方法進行研究[10-14]。然而，上述單向流固耦合瞬態計算方法都沒有考慮水體的附加質量，理論上對于動應力幅值特別是共振工況動應力的預測存在一定誤差。文獻[15]基于聲固耦合提出了改進的單向流固耦合(Acoustic-structure based one-way FSI, ASOW-FSI)方法，并以此預測了模型水泵水輪機轉輪在不同轉速工況下的動應力，結果表明計算得到的共振曲線與試驗數據高度吻合。文獻[16]借助圓柱渦激振動研究了ASOW-FSI方法在流固耦合問題中的適用性，結果表明這一方法不能預測頻率鎖定等自激振動現象，但對于受迫振動計算具有較高的精度。更多關于流固耦合問題的研究現狀可參見文獻[17-19]。

對于復雜的水力機械轉輪，動應力的獲取需要消耗大量的計算資源。而一個完整的過渡過程通常長達幾十秒，通過數值計算準確預測其動應力特性十分困難。因此，目前關于瞬態過程轉輪動應力數值計算的報道較少。本文采用基于聲固耦合構建的單向流固耦合方法研究某原型混流式水泵水輪機轉輪在啟動過程中的振動響應，重點分析水力激勵力的瞬態特性對動應力的影響。

1 數值模型

1.1 動態響應分析方法

線性系統的結構動力學方程可以表示為

(1)

式中Ms——結構系統的質量矩陣

Cs——結構系統的阻尼矩陣

Ks——結構系統的剛度矩陣

u(t)——節點位移向量

F(t)——外部激勵力

在流固耦合問題中，F(t)即為水壓力，可通過計算流體動力學(CFD)分析得到。

在進行雙向流固耦合計算時，流體流動產生的非恒定水壓力使得結構產生振動。同時，結構的運動也會反作用于流場。而在單向耦合中，忽略了結構運動對流場帶來的影響。事實上這種簡化不僅引起流場計算不準確，還會導致對結構動力特性的誤判。當結構在流體中變速運動時，將受到額外的與加速度方向相反的流體力。這種現象被稱為附加質量效應，現已被廣泛研究[20-23]。其中，文獻[20]通過試驗研究了靜水中淹沒平板的附加質量效應，結果表明附加質量是結構模態振型的函數，附加質量在低階模態尤為顯著。文獻[21]利用試驗和數值方法研究了水體附加質量對水泵水輪機轉輪動力特性的影響，并從振型分析出發討論了這種現象產生的機理。文獻[22-23]研究了壁面對附加質量的影響，結果表明壁面的接近程度可以極大地改變附加質量，結構離壁面越近，附加質量效應越明顯。其中，文獻[22]采用聲固耦合方法計算了淹沒圓盤的模態，并與試驗數據進行了比較。結果表明，聲固耦合方法可以準確預測考慮壁面影響的水下結構固有模態。

ASOW-FSI方法在動態響應分析中采用聲學單元對結構周邊的水體進行建模，可以較充分地考慮水體對結構動力特性的影響。此時，系統的結構動力學控制方程可以寫為

(2)

式中Ma——附加質量矩陣

Ca——附加阻尼矩陣

Ka——附加剛度矩陣

文獻[15]運用這種方法預測了模型水泵水輪機轉輪在不同相對轉速(即轉速與額定轉速之比N/Nrated)下的動應力，結果如圖1所示。可以看出，ASOW-FSI方法對共振轉速和共振幅值的預測都與試驗數據吻合較好。與常規的單向耦合(OW-FSI)方法相比，ASOW-FSI方法在共振點附近的計算精度得到了很大的提升。

1.2 流動分析模型

本文研究的混流式水泵水輪機組的基本參數如下：轉輪直徑Dr=5.1 m，葉片數Zb=9，額定轉速Nrated=300 r/min，固定導葉及活動導葉數Zs=Zg=20，導葉最大開度Gmax=36°，額定水頭Hrated=295 m，額定流量Qrated=118 m3/s。

為了分析啟動過程中水力激勵力的瞬態特性對動應力的影響，本文針對某個特定啟動過程選取不同轉速工況點分別進行穩態和局部瞬態動應力計算。其中，穩態動應力是指假定機組在某工況點穩定運行時的動應力水平；而局部瞬態動應力是指在過渡過程中選取某一時間段進行獨立求解得到的結果。機組啟動過程中轉速和導葉開度的變化規律如圖2(圖中G表示導葉開度)所示。在轉速上升段(圖中虛線框所示時間段)選取了4個工況點，對應的運行參數展示在表1(表中H表示運行水頭)中。

為了給動應力計算提供相對準確的邊界條件，建立了全流道三維非定常流動數值分析模型，其計算域包括蝸殼、固定導葉、活動導葉、轉輪流體域、間隙、均壓管及尾水管等部分。使用六面體和四面體網格對計算域進行離散，其中，六面體網格主要分布于間隙、均壓管、尾水管及活動導葉等區域。計算域網格如圖3所示，蝸殼及固定導葉區域的網格單元數8.6×105，活動導葉2.78×105，轉輪流體域1.499×106，間隙4.43×105，均壓管1.896×106，尾水管8.24×105。

表1 選定工況點的運行參數Tab.1 Operating parameters of selected working conditions

在蝸殼進口設置總壓邊界，尾水管出口為開放邊界，進出口壓力均由實測數據換算得到。設置轉輪流體域和上、下間隙為轉動域，其中靠近轉輪室一側的壁面保持靜止，而其它區域均設為靜止域，在各計算子域之間設置交界面(包括動-靜和靜-靜交界面)。為了保證交界面上數據傳遞的準確性，采用CFX提供的GGI網格連接方式進行插值。

選用SSTk-ω湍流模型求解URANS方程。采用有限體積法對控制方程在空間上進行離散，而時間上的離散采用二階全隱式格式。非定常計算的時間步長Δt設為轉動周期的1/200，待測點壓力呈現出明顯周期性后，至少再計算5個轉動周期。

針對該水泵水輪機組在額定工況下穩定運行時的水壓力進行了現場測試，重點監測了球閥前后、蝸殼以及尾水管進口處的壓力脈動(采樣頻率為100 Hz)。其中，蝸殼處測點的位置、現場實測得到的壓力脈動以及本文CFD計算的結果如圖4所示。

可以看出，本文通過CFD計算得到的壓力脈動的均值和幅值都與實測數據吻合較好。其中，均值誤差為0.08%，幅值誤差為7.55%，測點位置的偏離可能是造成壓力脈動幅值產生誤差的原因之一。本文所述“幅值”均指某一變量時域信號中的峰峰值。由于現場測試時的采樣頻率較低，其結果無法真實反映壓力脈動的頻率，因而不能直接用來判定數值模擬對于脈動頻率預測的準確性。對數值計算的結果重新進行了采樣(采樣頻率同樣為100 Hz)，結果也展示在圖4中。結果與實測數據基本吻合，證明了該數值分析模型對壓力脈動頻率的預測是準確的。

1.3 有限元模型

采用基于聲固耦合的單向流固耦合方法計算轉輪在運行過程中的動態響應，求解過程在ANSYS APDL中進行。計算域包括轉輪結構以及由聲學單元模擬的間隙流體和轉輪流體域。其中流體域沿用了流動分析中的網格，并且流固耦合交界面兩側的節點一一對應。重點針對葉片及其倒圓處進行了網格加密。有限元模型如圖5所示。

作用在轉輪結構上的荷載包括非定常水壓力、重力及離心力。其中，水壓力由流體動力學計算得到，加載的時間步長為轉動周期的1/200。在主軸螺栓中心線上施加固定約束，將轉輪進出口和間隙進出口的流體域邊界設為全吸收邊界，而頂蓋和座環壁面設為全反射邊界。轉輪材料為結構鋼，密度為7 850 kg/m3，泊松比0.3，彈性模量200 GPa。而水體的密度為1 000 kg/m3，水中聲速取1 482 m/s。

2 結果與討論

2.1 轉輪模態

為了解轉輪的動力特性，采用聲固耦合法計算了轉輪在流道中的固有模態。模態分析中的網格、邊界條件和材料屬性均與動應力分析保持一致。本文主要關注轉輪的節徑模態(節徑是指結構在自由振動過程中位移保持為零的直線)。圖6展示了前4階節徑模態的振型，其中1ND、2ND和3ND模態以上冠和下環的振動為主，而4ND模態的最大位移在葉片上。前4階節徑模態的固有頻率分別為24.84、57.76、119.04、161.90 Hz。

2.2 動靜干涉

由于高水頭水泵水輪機轉輪的進水邊(發電工況)設計較矮，且活動導葉的出水邊較厚，在運行過程中轉輪葉片和導葉間的干涉作用較強。文獻[24]研究表明，高水頭水輪機轉輪在最佳效率運行點附近時，動靜干涉貢獻的動應力約占總應力幅值的80%。因此，開展動靜干涉理論計算有助于理解和分析非定常水壓力以及結構的動態響應。

由轉輪葉片和活動導葉聯合作用產生的最終壓力場可以表示為旋轉壓力場和靜止壓力場的累加，經調制，最終壓力場在靜止參考系中可表示為[25]

(3)

其中

ω=2πN

(4)

式中ω——轉輪轉動角速度

N——轉輪轉速

m——轉動系統的諧波階數

n——固定系統的諧波階數

Amn——聯合幅值

θs——靜止坐標系中的角位置

φm——轉動系統對應階次諧波的相位

φn——固定系統對應階次諧波的相位

這個壓力場具有兩種徑向壓力模態，節徑數k分別為

k1=mZb-nZg

(5)

k2=mZb+nZg

(6)

一般情況下，諧波階次越高，振幅越小，因此可以只考慮k=k1的情況。其在旋轉坐標系和靜止坐標系下角速度和激勵頻率分別為[1]

ωr=nZgω/k

(7)

ωs=mZbω/k

(8)

fr=nZgfN

(9)

fs=mZbfN

(10)

式中ωr、ωs——旋轉坐標系、靜止坐標系下角速度

fr、fs——旋轉坐標系、靜止坐標系下激勵頻率

fN——轉動頻率

對于本文研究的對象，動靜干涉模式(部分結果)如表2所示。以2節徑振型為例：在靜止坐標系下，機組各部件受到的動靜干涉激勵頻率為18fN，角速度為-9ω；針對轉動部件進行分析時，需將其置于旋轉坐標系下來考慮，其2ND水力激振力的頻率為20fN，角速度為-10ω。當轉輪結構的2ND模態頻率等于20fN時，可能引發共振。

2.3 穩態動應力

首先針對前文選定的4個工況點進行穩態計算。在轉輪振動和動應力研究過程中，需分析其周邊的水壓力特性，因此設置了如圖7所示的7個壓力脈動測點。以OP1工況為例，不同測點處的壓力脈動時程曲線如圖8所示。由于p2和p3處的壓力脈動特性基本一致，在分析中略去p3的結果。

表2 動靜干涉模式(部分結果)Tab.2 Modes of rotor-stator interaction

從圖8可以看出，OP1工況下無葉區和葉片進水邊處的壓力脈動最顯著，其幅值分別為0.12 MPa和0.10 MPa。此外，葉片出水邊處也出現了明顯的壓力脈動，約為0.07 MPa。而間隙內的壓力脈動相對較弱，幅值約為0.04 MPa。受篇幅限制，本文沒有將OP2～OP4工況測點處的壓力-時間曲線列出。但是圖9展示了不同工況下各測點處的壓力脈動頻域信號。各工況下壓力脈動的頻率特征基本一致：無葉區的壓力脈動頻率為葉片通過頻率9fN及其倍頻；轉輪葉片進水邊及間隙處為導葉通過頻率20fN；出水邊的頻率成分相對復雜，有20fN、18fN以及由漩渦運動引起的壓力脈動頻率(如OP1中的2.33fN和OP4中的0.667fN)。其中，0.667fN對應尾水管渦帶的轉動頻率0.333fN，由于壓力測點處于旋轉坐標系下，監測到的脈動頻率與渦帶轉頻相差fN。

OP1至OP4工況下的轉頻fN分別為2.87、3.2、3.8、4.5 Hz，對應的20fN分別為57.4、64、76、90 Hz。其中，OP1工況的20fN約等于轉輪結構在流道中的2ND模態頻率57.76 Hz，結合動靜干涉理論分析可知，此時轉輪可能會發生共振。

完成CFD計算后，將轉輪表面的水壓力轉換為節點荷載，即可進行結構響應分析。圖10為OP1工況下某一時刻轉輪的等效應力云圖。可以看出，葉片與上冠、下環相接處出現了明顯的應力集中，其中葉片進水邊靠近上冠處最為明顯，最大等效應力為317 MPa。在后續分析過程中，分別在葉片進水邊靠近上冠處和靠近下環處設置動應力監測點pc和pb。

用測點處每一時刻的等效應力減去時均值，即可得到動應力時程曲線。不同工況下pc和pb處動應力如圖11所示。可以看到，在OP1工況，轉輪產生了明顯的共振，pc和pb處的動應力幅值分別達到了47.9 MPa和47.1 MPa，并且還有進一步增大的趨勢。此時，轉輪結構的2ND模態被激發，振動頻率為20fN。而在轉速較大的OP2、OP3、OP4工況，轉輪的振動較弱，pc處的動應力幅值分別為7.5、6.2、6.6 MPa。其中OP3和OP4工況下動應力的頻率成分相對簡單，主頻為20fN。而在OP2工況，轉輪的振動形式較為復雜，上冠和下環處的動應力主頻分別為18fN和20fN。

事實上，OP2工況的轉頻fN為3.2 Hz，而靜止坐標系下機組的動靜干涉頻率18fN=57.6 Hz，恰好對應轉輪結構的2ND模態頻率，但是其激勵振型與結構模態振型存在一定偏差。因此，轉輪的固有頻率在結構響應中表現突出，但不會引發共振。根據動靜干涉分析結果可知，在旋轉坐標系下，節徑數為2的水力激振力頻率為20fN，轉速為-10N，每隔時間Δt=1/(20fN)激勵一次，并且在此過程中激振力旋轉了半個周期，而轉輪結構在旋轉坐標系下保持靜止，因此恰好可以激發出結構的2ND振型；同理，在靜止坐標系下，節徑數為2的激振力頻率為18fN，轉速為-9N，每經過時間Δt=1/(18fN)，激振力旋轉半個周期，但是由于轉輪本身的轉動，每一次激勵都會產生20°的偏差。從另一個角度來說，即18fN對應的激振力在旋轉坐標系下并不是嚴格對應2ND振型。上述分析可以用圖12進行描述。

2.4 瞬態動應力

實際上在啟動過程中，轉輪所受的水力激勵力是非周期性的。總體上，隨著轉速逐漸增大，激勵力的脈動頻率越來越高。針對上述4個工況點進行局部瞬態計算，分析激勵力的瞬變特性對轉輪動應力的影響。首先，轉速在OP1工況(即共振點)附近時危險點pc處的瞬態應力計算結果如圖13所示。從應力水平(即時均應力)的角度來看，隨著轉速不斷增大，等效應力水平逐漸上升。從動應力的角度看，當轉速增至OP1工況時，并沒有出現明顯異常的應力波動。但是隨著時間推移，應力的波動逐漸增強，達到峰值后逐漸減弱。最大應力波動出現的時間較共振點滯后約0.25 s。這一現象是由轉輪結構的慣性引起的。另外，應力波動的幅值約為24 MPa，比穩態計算的結果(47.9 MPa，見圖11a)小得多。這是因為結構的振動需要一定時間才能激勵出來，而在機組啟動過程中激振力的頻率一直在增大，與轉輪固有頻率重合的時間較短，所以并未產生非常強烈的振動。由此可以推斷，啟動過程中轉速變化得越慢，在共振點附近能夠激勵出的振動越強烈。上述結果表明，在共振點附近采用穩態方法估算過渡過程中的動應力水平是不可行的。

接下來，以OP4工況為例，分析非共振區域的瞬態計算結果(如圖14a所示)。同樣，可以看出隨著轉速不斷增大，應力水平逐漸上升。實際上，這是由于轉輪離心力和流場中的壓力梯度逐漸增大所致。為了剔除應力水平的變化對動應力幅值評估的干擾，假定應力水平與轉速在局部呈線性關系，然后采用最小二乘法進行線性回歸分析，結果如圖中的虛線所示。由此，可以得到轉速在OP4工況附近時pc處的動應力，如圖14b所示。結果表明，動應力幅值為6.7 MPa，與穩態計算的結果(6.6 MPa，見圖11d)非常接近。說明在非共振區域可以采用穩態方法估算過渡過程中的動應力水平。

根據上述分析可知，可以采取穩態和局部瞬態相結合的方法來研究過渡過程中結構的動應力。具體的實施過程如下：首先，進行非穩態流動分析；第2步，進行流固耦合結構模態計算；第3步，找到共振工況點以及由于導葉開度、轉速等運行參數瞬變引起高水平壓力脈動的工況點；第4步，針對上述工況點進行局部瞬態動應力計算，其它區域可選擇性地取點進行穩態分析。

3 結論

(1)采用基于聲固耦合構建的單向流固耦合方法研究了某原型混流式水泵水輪機組在水輪機模式啟動過程中轉輪的振動響應和動應力。結果表明，在非共振區域，激振力的瞬變特性對轉輪的振動幾乎沒有影響，具體表現為穩態計算和瞬態計算得到的動應力水平基本一致；而在共振區域，激振力的瞬變特性對轉輪的振動有顯著影響，一方面會使啟動過程中的動應力幅值遠比采取穩態方法估算的值低，另一方面，還會使得高水平振動滯后于共振點。上述結果可為過渡過程轉輪動應力的評估提供參考。

(2)借助模態分析和動靜干涉理論對轉輪的共振特性進行了分析。研究表明，僅當某一水力激勵力與轉輪固有模態的頻率和振型均對應時才會產生共振；當頻率相等而振型之間存在相對轉動時，這種激振力會在一定程度上主導轉輪的振動，但振動幅值只是較常規情況略微增大。