某型波浪發電裝置3D打印樣機飛逸轉速的風洞測試與數值仿真

鄭貴林，梁燦棉

(武漢大學動力與機械學院，武漢 430072)

0 引 言

憑借存在范圍廣、儲量大、能量形式為機械能等特點，波浪能被視作可再生自然能源開發利用的重要一員，各種形式的波浪能采集裝置也不斷涌現。而海洋具有表面波浪起伏，海面20～30 m以下平穩的特點[1,2]，于是差動垂直波浪能提取發電裝置應運而生。

該發電裝置的特點是利用導流葉片使內部旋轉葉輪在海水相對其從上向下流還是從下向上流都保持一個旋轉方向[3,4]。

目前，對于因受到外界流體的作用而被動轉動的旋轉葉輪例如風力發電機和水輪發電機組等在某一特定負載力矩下的旋轉速度的預測的研究尚少，一般都是通過具體的樣機試驗測定，而葉輪的旋轉速度又是決定變速機構的設計的一個重要參數，如果能夠使用CFD軟件準確地模擬計算出葉輪在不同流體速度和不同負載力矩情況下的旋轉速度，將對變速機構的設計帶來極大的方便，并且方便對葉輪的效率分析以及進一步的優化設計工作。近年來，國內外有部分學者使用CFD軟件來進行水輪機飛逸狀態(負載力矩為零的狀態)的三維數值模擬并取得可觀的成效。周大慶等[5]對軸流轉槳式水輪機在不同槳葉啟閉規律下的飛逸狀態進行計算，得出飛逸轉速在不同槳葉啟閉規律下的變化情況以及尾水管負壓值、壓力脈動等的變化。李金偉等[6]對某型混流式水輪機在不同開度下的飛逸狀態進行三維非定常計算，得出飛逸狀態下單位轉速、單位流量以及水力矩的變化情況。羅興锜等[7]對貫流式水輪機的飛逸狀態進行數值計算，獲得飛逸轉速、流量、力矩、軸向力等外特性參數以及水輪機內部流場的動態特性。

如今對于利用海洋能量來進行發電的被動型旋轉葉輪的飛逸狀態的三維數值模擬的研究并不多見，差動垂直波浪能提取發電裝置其中的能量獲取機構便是一種旋轉葉輪，對該旋轉葉輪的轉速預測可以為變速機構的設計、效率分析以及優化設計提供參考指導作用。由于對大尺寸的發電裝置進行海試實驗的成本高、操作難度大、改進的研發周期長等，而Fluent等一些計算流體力學軟件進行三維數值仿真具有成本低、試驗方便、改進研發周期短的特點，因此可以使用Fluent為差動垂直波浪能提取發電裝置的優化改進提供一定的輔助作用。另外隨著3D打印技術日漸成熟，能夠打印出具有較為復雜表面的實體，因而便借助3D打印來制作小尺寸的發電裝置樣機來進行實物實驗來驗證裝置模型設計的正確性和數值計算方法的可行性等[8]。而要進行實物實驗除了需要發電裝置樣機還需要實驗環境，由于發電裝置的工作過程與風力發電機的工作過程比較類似，即是利用流體與裝置之間的相對運動的能量，另外波浪的垂直方向速度曲線近似梯形波時存在速度比較平穩的階段，因而決定使用相對較簡單的風洞實驗來取替水槽實驗來模擬起伏運動中速度比較平穩的階段并近似為單向恒速狀態。為此，本文通過Fluent軟件建立差動垂直波浪能提取發電裝置中旋轉葉輪在單向恒定流速空氣流作用下的飛逸轉速的三維湍流數值計算方法，并與樣機的風洞測試實驗進行對比，以驗證數值計算方法中被動型動網格、湍流模型等參數設置的正確性，為后續雙向變流速的水流作用下以及帶負載作用下的仿真和優化提供參考指導。

1 差動垂直波浪能提取發電裝置的結構、工作原理以及3D打印的樣機

1.1 結構組成及工作原理

差動垂直波浪能提取發電裝置的總體結構如圖1所示，其中包括海面漂浮平臺、漂浮平臺上安置的發電機、傳動轉軸和能量轉換波輪機，而波輪機由導流圓管、導流葉片和轉子等組成。

差動垂直波浪能提取發電裝置的運行原理為：如圖1所示，漂浮平臺位于海洋表面波浪區，波輪機置于穩定區。漂浮平臺隨波浪上下運動并通過傳動轉軸帶動波輪機在垂直方向上下運動，因此波輪機與其周圍穩定區中的水流存在垂直方向上的相對往復運動，把參考系放在波輪機上，便可以看作穩定區中的水流在垂直方向上下往復地流經波輪機。水流在上下往復的過程中，經由上部和下部導流葉片的引導始終單方向作用于轉子葉片從而使轉子單方向地旋轉，再經由傳動轉軸帶動發電機旋轉發電。

圖1 差動垂直波浪能提取發電裝置示意圖Fig.1 The differential vertical wave energy converter

1.2 3D打印的樣機

由于使用風洞實驗取替水槽實驗，直接有空氣流體與波輪機作相對運動，不需要漂浮平臺來帶動波輪機做上下起伏的運動，故只使用3D打印機打印差動垂直波浪能提取發電裝置中的波輪機部分。打印的波輪機主要參數如下：導流圓管最大直徑：D1=100 mm；導流圓管最小直徑：D2=80 mm；轉子輪轂直徑：D=20 mm；導流葉片數量：Z1=12或16；轉子葉片數量：Z2=12或16；轉子葉片截面形狀：圓弧形和角形。打印完成的轉子零件以及進行組裝所得的波輪機樣機如圖2，圖3所示。

圖2 圓弧形葉片轉子以及角形葉片轉子Fig.2 Rotor with circular arc blades and rotor with folded angle blades

圖3 完整組裝的波輪機樣機Fig.3 The whole propotype

2 數值計算方法與模型

2.1 運動方程及轉速計算方程

在波輪機轉子飛逸過渡過程中，其外負載為0，轉子被水流沖擊而轉速迅速上升，其運動方程為：

(1)

式中：Mt為主動力矩，N·m；Mg為負載力矩，N·m；J為機組總轉動慣量，kg·m2；ω為轉子旋轉角速度，rad/s。

每一步的轉速ωi+1可以通過下式來求得：

(2)

式中：Δt是計算時間步長，s。

2.2 湍流模型與邊界條件

本文采用商業CFD軟件Fluent16.0來模擬波輪機在風洞測試儀中的飛逸過程，計算中采用非定常(即瞬態)的標準k-ε湍流模型。流體介質選擇為空氣，密度為1.225 kg/m3，動力黏度為1.789 4×10-5kg/s，進口采用速度邊界條件，流體速度設置為5 m/s，方向垂直于進口截面。出口采用自由出流邊界條件。由于轉子轉速是根據式(2)來確定，屬于被動轉動，因此采用6DOF被動型動網格模型來設置轉子所在的流體域，其他部分的流體域則為靜流體域，不同的流體域之間的重合面設置為interface交界面。通過編寫用戶自定義函數(UDF)來輸入機組的總轉動慣量、質量、運動約束條件等參數。計算過程中以靜止狀態為初始值，時間步長設為0.001 s。

由于本文的目的是通過對波輪機樣機進行單向恒速空氣流風洞測試與仿真計算來對比以驗證被動型動網格、湍流模型等參數設置的正確性，因此對于實際環境中的雙向變流速水流的邊界條件未作考慮。對于垂直方向速度曲線近似為梯形波的波浪，由于存在速度比較平穩的階段，此階段內的力矩特性可以用單向恒速流體作用下的力矩特性來近似，假若速度曲線梯形波的斜線段時間很短即速度曲線接近矩形波的時候，整個周期內的力矩特性都能夠用單向恒速流體作用下的力矩特性來近似。工作介質為空氣時，計算結果所得的力矩值小于0.001 N·m，現有的測量設備難以測量如此小的力矩值，因此風洞測試的力矩值沒有進行測量。計算工作介質改為水時的飛逸轉速有約10%的提高，力矩值提升約三個數量級，與水的密度和空氣的密度比值相接近。

2.3 流體場網格

按照風洞測試儀的尺寸來建立直徑為25 cm，長度為53 cm的圓柱體風洞流體場，波輪機軸線與該圓柱體軸線重合，波輪機前端距離風洞流體場入口截面16 cm，后端距離風洞流體場出口截面30 cm，對于圓柱體風洞流體場可以通過切割分塊并采用六面體網格劃分，如下圖4所示。由于轉子葉片和導流葉片的幾何外形比較復雜，故這些葉片附近的流體域采用非結構化四面體網格來劃分，而距離葉片較遠的區域則切割分塊并采用六面體網格劃分以減少網格數量，如圖5所示。

圖4 風洞流體場網格劃分示意圖Fig.4 Flow field grid division

圖5 波輪機內部流體域網格劃分示意圖Fig.5 Internal flow field grid division of the converter

2.4 網格無關性驗證

為減小網格數量對仿真計算結果的影響，分別選取網格數為442萬(具體是4 423 871)、513萬(具體是5 134 559)、550萬(具體是5 503 763)、606萬(具體是6 061 361) 這4種網格模型對進行了計算，以波輪機轉子飛逸轉速為目標變量來進行網格無關性驗證，驗證結果如圖4所示。由圖6可知，網格數量達到513萬以上波輪機轉子飛逸轉速幾乎無變化，考慮計算精度和計算經濟性，確定模型網格數量為513萬。

圖6 網格無關性測試結果Fig.6 Gird independence test results

3 數值計算結果與風洞測試結果

3.1 飛逸轉速對比

對定子葉片數量與轉子葉片數量均為12和定子葉片數量與轉子葉片數量均為16，轉子葉片形狀分別為圓弧形和折角形這四種樣機在風速5 m/s下的飛逸轉速的數值計算結果與風洞測試結果如表1所示。從表1中可以看出數值計算所得的飛逸轉速與風洞測試所得的飛逸轉速相對誤差在10%以內，說明數值計算飛逸轉速的可行性，以及圓弧形葉片的轉子的飛逸轉速比折角形葉片的轉子的飛逸轉速更高。

表1 兩種形狀葉片轉子的數值計算飛逸轉速與風洞測試飛逸轉速的對比Tab.1 Comparison of the runway speeds of the numerical simulations and the wind tunnel tests(rotors with circular arc blades and folded angle blades)

對于數值計算所得的飛逸轉速與風洞測試所得的飛逸轉速的誤差分析，作者認為有主要兩個部分：第一部分，由于用于風洞測試的樣機內部安裝有軸承，雖然是測試飛逸狀態也就是不帶外負載力矩，但是軸承的滾動摩擦力矩依然會對樣機的飛逸轉速產生輕微的影響，而在數值計算中沒有考慮對軸承滾動摩擦力矩的設置；第二部分，3D打印所得的樣機模型與在三維繪圖軟件中的原始設計模型仍然存在一些外形和尺寸誤差，而在數值計算中使用的是三維繪圖軟件中的原始設計模型。

3.2 轉動慣量與飛逸轉速關系分析

對定子葉片數量與轉子葉片數量均為12和定子葉片數量與轉子葉片數量均為16，轉子葉片形狀為圓弧形，但轉子轉動慣量不一樣的4種樣機在風速5 m/s下的飛逸轉速的數值計算結果與風洞測試結果如表2所示。

從表2中可以看出，轉子轉動慣量的改變對飛逸轉速的影響幾乎沒有(在不考慮因為轉動慣量變化帶來機械摩擦變化的情況下)。經過仿真得到該裝置的力矩-轉速曲線如下圖7所示。結合力矩-轉速曲線以及式(2)顯而易見，當Mt=Mg=0以后，ωi+1=ωi，即轉子最終穩定的轉速總是達到力矩為0時對應的轉速值也就是飛逸轉速值，而轉動慣量主要影響了轉子從初始狀態到飛逸狀態之間變化的時間，該時間長度T=iΔt，i為式(2)達到飛逸狀態迭代計算的所使用的步數。

表2 不同轉動慣量轉子的數值計算飛逸轉速與風洞測試飛逸轉速的對比Tab.2 Comparison of the runway speeds of the numerical simulations and the wind tunnel tests(rotors with different moments of inertia but the same circular arc blades)

圖7 仿真計算的力矩-轉速曲線Fig.7 Curves of torque against rotate speed

4 結 語

通過對差動垂直波浪能提取發電裝置中的能量轉換機構波輪機樣機的飛逸轉速進行單向恒速風洞測試以及使用Fluent軟件進行三維數值仿真，得出以下結論。

(1)樣機在風洞測試中的飛逸轉速值與Fluent仿真計算所得的飛逸轉速值接近，誤差在10%以內，說明數值計算的可行性。

(2)在其他條件相同的情況下，圓弧形葉片的轉子與折角形葉片的轉子相比具有更高的飛逸轉速。

(3)保持其他條件不變的情況下，轉子轉動慣量的變化對轉子的飛逸轉速的影響幾乎沒有。