潮流能水平軸水輪機葉片雙向流固耦合特性分析

袁 鵬， 徐孝輝， 譚俊哲， 張金輝， 于鵬鵬， 王樹杰

（1.中國海洋大學 工程學院， 山東 青島 266100； 2.青島市海洋可再生能源重點實驗室， 山東 青島 266100）

0 引言

潮流能水輪機運行于復雜的海洋環境中，在水流的影響下，水輪機的葉片結構會發生水彈性變形， 而葉片結構的變化又反作用于水流的流動，對流場產生影響，干擾水輪機的正常運轉，進而影響發電效率[1]。 如果這種耦合作用越來越大，葉片會喪失穩定性而產生裂紋甚至斷裂，這要求葉片要有較好的強度和剛度，從而對葉片結構的設計提出了較高要求[2]，[3]。 因此，進行流固耦合條件下的水輪機葉片的性能分析和結構優化，對水輪機的安全穩定運行具有重要意義。

Y Bazilevs 采用風機真實工況下的風速和轉速對海上風力機進行了全尺寸的雙向耦合模擬，得出了風機在此工況下的動力性能以及結構性能[4]。Jong Won Lee 使用流固耦合方法對風力機葉片的氣動彈性進行了研究，結果表明，隨著風速和載荷的增加，葉片的彈性變形逐漸增大，氣動外形也會發生明顯改變，從而影響葉片的性能[5]。荊豐梅采用ANSYS 和CFX 軟件進行了垂直軸潮流能水輪機的單向流固耦合瞬態計算，并完成了水輪機工作狀態主軸應變計算，且計算結果與實驗結果吻合較好[6]。 胡丹梅采用ANSYS 軟件對5 MW 離岸型風力機葉片進行了流固耦合計算，分析了額定風速下葉片的變形、推力和扭矩等參數的變化和耦合條件下葉片周圍流場的變化[7]。 程友良建立了海上風力機塔架葉片的耦合模型，結合Davenport 脈動風速譜， 利用雙向流固耦合方法， 對塔架葉片的耦合結構進行了動力學參數變化分析， 發現塔架葉片耦合后的葉尖位移偏移量呈脈動性變化，塔架的慣性矩較小[8]。

本文建立了20 kW 潮流能水輪機葉片的實尺度模型，采用ANSYS Fluent 模塊進行了三維流場計算， 使用Transient structure 模塊進行了葉片結構響應計算， 在每個時間步上采用System coupling 模塊進行了流固耦合數據交換， 得到了葉片在實際工況下的水動力學性能和結構性能，通過與不考慮流固耦合條件的計算結果進行對比分析， 得到了流固耦合條件下的葉片性能的變化規律，并將仿真結果與實驗結果進行對比，從而驗證計算結果的準確性， 為潮流能水輪機葉片的設計提供參考。

1 流固耦合的計算理論

流固耦合問題按其耦合機理可以分為兩大類， 一種是流固耦合作用僅發生在流體和固體交界面處，另一種是發生在流體和固體的重疊部分，并且兩相域重疊在一起難以明顯分開[9]。本文研究的內容屬于第一類問題。

1.1 流固耦合方程

基于流體力學三大方程和理想流體及小變形假設， 運用Galerkin 法對流體運動微分方程離散化，可以得到離散化后的流體運動方程式[10]：

式中：H，A，E 和B 均為系數矩陣；ρ 為流體密度；p 為壓力矢量；q0為輸入激勵矢量；r 為位移矢量；C 為流體的壓縮波速度；N 為形狀函數矢量；Ω 為流體域體積；Sr 為無限遠邊界處的表面積；SF為自由表面處的表面積；g 為重力加速度。

與流體接觸的結構運動方程式為

1.2 ANSYS workbench 耦合計算流程

ANSYS 提供的System coupling 模塊可以實現Fluent 和Transient structure 的雙向流固耦合計算。 單向流固耦合分析在耦合交界面處的數據傳遞是單向的，而雙向流固耦合計算既有CFD 分析計算的結果（如壓力和溫度）傳遞給固體結構分析，又有固體結構分析的結果（如位移和速度）反向傳遞給流體分析。 由于考慮了結構變形對流場的擾動，雙向流固耦合計算可以更真實地模擬水輪機葉片的運轉情況。

雙向流固耦合的計算采用分離解法[11]，即在每一個時間步長內，按照順序首先在Fluent 求解器中求解流場，通過流固耦合面傳遞壓力值給結構場，然后求解結構場，在流固耦合面處返回位移信息給流場，最后檢查流場、結構場和流固耦合面處是否達到各自的收斂要求， 若不收斂則返回第一步重新計算，直到達到收斂標準或最大迭代數，若收斂則進行下一時刻的計算。

2 計算模型的建立

2.1 葉片幾何模型計算模型

選定齋堂島海域作為20 kW 潮流能水平軸水輪機機組的安裝區域， 該島西北部與陸地形成了一條水道，該水道深8～14 m。 水輪機的額定功率為20 kW，額定流速為1.5 m/s，設計使用年限為20 a，葉片翼型為NACA63-8 系列。

水輪機的整體參數如表1 所示。

表1 水輪機整體設計參數Table 1 Overall parameters of the turbine

2.2 流場模型與結構模型的建立

圖1 為葉片旋轉域與流體域劃分圖。 流體域的高度為15 m，長度為30 m，入口邊界條件為速度入口， 流速為1.5 m/s， 水輪機葉片距離入口8 m，出口邊界條件為自由出口。 葉片旋轉域的半徑為2.8 m，長度為1 m，流體域與旋轉域的交界面采用滑移網格，葉片表面設定為流固耦合交界面，采用適合水輪機的k-ω SST 模型進行計算求解[12]。

圖1 葉片旋轉域與流體域劃分圖Fig.1 Blade rotation domain and fluid domain

在葉片流固耦合問題中， 由于固體結構存在變形運動， 所以不僅需要耦合求解結構動力學和流體動力學方程， 還要采用動網格技術不斷調整計算網格來準確描述不斷變化的流體和固體幾何空間。 在葉片旋轉域采用的彈簧光順方法是一種實現非結構網格隨體變形的有效方法， 計算效率較高[13]。

圖2 葉片結構模型Fig.2 The model of blade structure

在Solid Works 中建立葉片結構模型（圖2），葉片內部采用雙腹板加橫向隔板支撐。 其中葉片蒙皮采用45°環氧E-玻璃纖維復合材料鋪層，葉片支撐結構為結構鋼。

2.3 網格無關性驗證

網格劃分的質量決定了計算結果的精度，也會影響計算機速度。 葉片表面采用三角形網格劃分，流體域和旋轉域采用四面體網格劃分。 在穩態條件下，當網格總數為125 萬～456 萬時，輸出功率與設計功率的相對誤差如表2 所示。 由表2可知，隨著網格數量的增加，相對誤差逐漸減小。綜合考慮計算量和計算精度， 選擇網格總數為386 萬，其中流體域網格數為313 萬左右，旋轉域網格數為73 萬左右。

表2 網格無關性驗證Table 2 Mesh independence verification

3 結果與分析

3.1 葉片水動力學性能分析

在相同的工況下，對考慮流固耦合條件以及不考慮流固耦合條件的葉片水動力學性能進行分析。

3.1.1 葉片獲能效率

在尖速比不同的情況下，考慮流固耦合條件和不考慮流固耦合條件的葉片的獲能效率如圖3所示。 從圖3 可以看出，考慮流固耦合條件時，葉片的獲能效率低于不考慮流固耦合條件時。 這是因為在流固耦合條件下， 葉片會發生水彈性變形，同時由于前后緣壓差的作用，葉片形狀和功角會發生改變，葉片會發生扭轉變形，導致葉片獲能效率降低。 由此可見，在流固耦合條件下，葉片獲能效率的計算值更接近于真實情況，且獲能效率低于設計值。

圖3 不同尖速比下，流固耦合和無流固耦合葉片的獲能效率對比Fig.3 Comparison of energy efficiency between coupled and uncoupled blades with different tip speed ratios

3.1.2 葉片截面表面的壓力分布

當尖速比為5 時， 在流固耦合與無流固耦合條件下，沿葉片展向26%，39%，53%，66%，79%和92%處的葉片截面表面的壓力分布如圖4 所示。

圖4 葉片截面表面的壓力分布Fig.4 The pressure of blade surface

從圖4 可以看出，沿葉片展向，流固耦合條件下的葉片上、 下表面壓差與無流固耦合條件下的葉片上、 下表面壓差的差值呈現出先增大后減小的變化趨勢。 沿葉片展向39%和53%處，兩種條件下的葉片上、下表面壓差的差值最大，這是因為這兩處的流固耦合作用最顯著， 在流固耦合條件下，葉片發生水彈性變形，葉片中部以及靠近葉根處應變值較大，導致葉片形狀發生一定改變，影響葉片表面流場， 從而導致迎流面壓力降低， 葉片上、下表面壓差減??；而在無流固耦合條件下，葉片默認為剛體，不會產生變形，葉片表面流場不會受到結構變形的影響。 在葉片展向39%以下和53%以上，兩種條件下的葉片上、下表面壓差逐漸減小，這說明在此長度以下葉片形狀變化不大，流固耦合作用不明顯，流固耦合條件下的葉片上、下表面壓差與不考慮流固耦合條件時相差不大。

3.2 葉片結構力學性能分析

圖5 為葉片的總體變形圖。從圖5 可以看出，葉片變形位移較大的區域在葉片中上部， 而且越靠近葉尖，累計變形位移越大。葉片最大變形位移可達44.845 mm，因此，在設計之初就應該考慮到葉片在流固耦合條件下的變形。

圖5 葉片總體變形圖Fig.5 Overall deformation of the blade

葉片在流固耦合條件下的應力及應變如圖6所示。

圖6 葉片應力和應變云圖Fig.6 Blade stress and strain diagram

從圖6 可以看出， 應變和應力主要集中在葉片中部和翼型葉根過渡區， 即葉片中部和翼型葉根過渡區的應變較大， 這會對葉片表面的流場造成影響，也進一步驗證了圖4 中數據的合理性。葉片內部最大應力為14.157 MPa，未超出玻璃纖維材料的彈性極限，葉片強度滿足要求。

圖7 葉片最大應變的變化曲線Fig.7 Maximum strain curve of the blade

耦合計算期間葉片最大應變隨時間的變化曲線如圖7 所示。剛開始計算時，葉片突然承受水流施加的載荷而產生了較大的應力應變， 同時會出現較大位移，之后在持續水流的沖擊下，初始擾動逐漸被消耗， 葉片的最大應變表現出圍繞平衡位置的低幅高頻振動。

3.3 驗證實驗

為驗證流固耦合計算對20 kW 潮流能水平軸水輪機水動力學性能的數值模擬結果的正確性，使用1∶5 的模型葉片進行水槽實驗，根據相似原理將實驗結果和數值模擬結果進行對比， 以驗證數值模擬結果的正確性。 本次實驗在中國海洋大學海洋工程實驗室進行，實驗水槽（圖8）長30 m，寬1 m，高1.2 m。 實驗須要使用流速儀監測流速，使用扭矩儀測量轉速和扭矩，以便計算獲能效率。

圖8 實驗水槽Fig.8 Test sink

通過實驗測得了水輪機的轉速和扭矩， 根據獲能效率計算公式和相似原理得到了水輪機的獲能效率（圖9）。 從圖9 可以看出：考慮流固耦合條件時的計算結果更加接近于實驗結果， 驗證了流固耦合計算方法的正確性； 實驗結果略低于仿真結果，這是由于數值仿真是在理想條件下進行的，而水輪機在實驗過程中會受到擾動， 實驗結果也會受到測量裝置誤差的影響。

圖9 不同尖速比下，獲能效率實驗結果與仿真結果的對比Fig.9 Comparison of energy efficiency between test and simulation with different tip speed ratios

4 結論

本文建立了20 kW 潮流能水輪機實尺度流固耦合分析模型，研究了葉片水彈性變形對葉片水動力學性能和結構性能的影響， 得到以下結論。

①葉片在流固耦合條件下的獲能效率計算值低于無耦合剛體葉片的計算值，且流固耦合條件下的計算值更加接近于真實情況。

②沿葉片展向， 流固耦合條件下的葉片上、下表面壓差與無耦合狀態下的葉片上、下表面壓差的差值呈現出先增大后減小的變化趨勢，在葉片展向39%和53%處，這一差值最大。

③葉片在流固耦合條件下發生水彈性變形，在葉片中部以及中部靠近葉根處應力應變較大，這是造成流固耦合條件下葉片中下部上、下表面壓差與無耦合狀態相差較大的主要原因。