， ， ，，

(1.中海油研究總院，北京 100027;

2.中國石油大學博士后流動站，北京 102200;

3.中國石油大學機械與儲運工程學院，北京 102200;

4. 哈爾濱工程大學海洋可再生能源研究所，哈爾濱 150001)

0 引 言

潮流能作為一種新興的清潔、可再生能源正在為世界各種爭相開發(fā)，而垂直軸水輪機作為潮流能轉換的核心裝置之一，已成為潮流能開發(fā)的主要途徑。水輪機技術的研究主要來自實驗驗證與計算機的數值模擬。垂直軸水輪機和垂直軸風力機基本原理是相同的，在數值模擬方面可以相互借鑒。2005年，Horiuchi等人[1]使用STAR-CD 模擬二維假設下的垂直軸風力機速度場，并與試驗值進行比較。2007年，Guerri等人[2]使用STAR-CD模擬二維垂直軸風力機葉片的受力。2008年孫科提出應用滑移網格理論研究H型葉輪的動邊界CFD數值模擬問題[3-4]。2009年，Hwang等人[5]應用STAR-CD，采用k-ε湍流模型對擺線式垂直水輪機進行了研究。

CFD數值模擬的方法在水輪機技術研究過程中得到了廣泛的應用，但是模擬結果的精度及可信度一直是所有研究者共同關注的焦點。盡管垂直軸水輪機或風力機CFD數值預報方面取得了很大的發(fā)展，但是有限的精度仍舊迫切需要得到改進。這其中如物理模型、網格劃分、邊界條件、湍流模型、CFD求解器等因素引起的計算結果中的誤差，都是不能被忽略的。因此，本文針對CFD方法在垂直軸水輪機數值模擬中的有效性進行了詳細的研究，包括湍流模型、網格尺度，以及時間步長對計算精度的影響。

1 CFD理論基礎

計算流體力學(Computational Fluid Dynamics，簡稱CFD)是流體力學研究領域中的重要技術之一，其基本思路為使用數值方法在計算機中對流體力學的控制方程進行求解，從而達到模擬預測流場的目的。下面介紹本文中涉及的湍流模型及其特點，公式的推導過程可以參考Ansys CFX用戶使用手冊[6]。

1.1 k-ε模型

(1)標準的k-ε模型：

標準k-ε模型：是從實驗現象總結出來的半經驗公式，也是最簡單的兩方程模型，主要基于湍流動能k和湍流耗散率ε。k方程是個精確方程，ε方程是個由經驗公式導出的方程。該模型假定流場完全是湍流，分子之間的粘性可以忽略。因此，只對完全是湍流的流場有效。

(2)RNGk-ε模型：

1.2 k-ω模型

(1)標準的k-ω模型：

兩方程模型，包含了低雷諾數影響、可壓縮性影響和剪切流擴散，因此適用于尾跡流動計算、混合層計算、射流計算以及受壁面限制的流動計算和自由剪切流計算。

(2)SSTk-ω模型：

它綜合了k-ω模型在近壁區(qū)計算的優(yōu)點和k-ε模型在遠場計算的優(yōu)點，將k-ω模型和標準k-ε模型都乘以一個混合函數后再相加就得到這個模型。與標準k-ω相比，SSTk-ω模型中增加了橫向耗散導數項，同時在湍流粘性定義中考慮了湍流剪切應力的輸運過程，模型中使用的湍流常數也有所不同。這些特點使得 SSTk-ω模型的適用范圍更廣，比如可以用于帶逆壓梯度的流動計算、翼型計算、跨音速激波流動計算等等。

1.3 雷諾應力模型(RSM)

雷諾應力模型中沒有渦粘度的各向同性假設，因此從理論上來說湍流模式理論要更精確。從理論上說，RSM 應該比一方程和兩方程模型的計算精度更高，但實際上雷諾應力模型的精度受限于模型的封閉形式，因此雷諾應力明顯具有各向異性的特點時才必須使用 RSM，比如龍卷風、燃燒室內流動等帶有強烈旋轉的流動問題。

2 垂直軸水輪機的CFD有效性驗證

2.1 計算算例

1983年Strikland等人[5]等人做的一系列研究垂直軸風力機動態(tài)失速試驗中，測量和記錄的數據包括葉片的瞬時切向和法向受力、葉片表面壓力分布以及尾流流場等，試驗結果作為一個可靠的參考依據。 Deglaire等人[7]、Hwang等人、Tchon和Paraschivoiu[8]和汪魯兵[9]的模擬研究中都是與其比較得出結論的。本文研究同樣采用Strikland等人的試驗模型為驗證算例。模型的具體參數如表1所示。

表1 Strickland試驗輪機的主要參數

2.2 網格劃分

網格劃分是CFD模擬過程中較為耗時的環(huán)節(jié)，也是直接影響模擬精度和效率的關鍵因素之一。對于所有的CFD數值模擬，結構網格模擬精度和速度要明顯優(yōu)于非結構網格情況，因此本文系采用結構化網格對流場進行劃分。根據垂直軸水輪機的結構特點，葉片在展長方向上具有相同的橫截面，為了節(jié)省計算時間，本文簡化為二維模型。

水蛭交配一般發(fā)生在早上或夜間，交配后一個月左右，雌體生殖器分泌出稀薄的黏液，中包被卵帶，形如“蠶繭”，排出體外，在濕泥中孵化。產繭一般在土壤（含30%水分）中進行，產繭時間一般為4月中下旬溫度在20-25℃以上。

劃分網格之前，還需要對模型中的流場進行基本的預估。過疏或過密的網格都會對計算有極大影響。過疏的網格往往會得到不精確甚至完全錯誤的結果；過密的網格會使計算量增大，使計算難以收斂。文中的計算網格如圖1和圖2所示。

圖1 Turbine 網格分布圖

圖2 Turbine 葉片周圍網格分布圖

2.3 邊界條件的設置

流體域內是以大氣壓為參考氣壓的；在進口邊界處設置為速度進口條件，給定均勻來流速度V及湍流參數；出口邊界采用開敞邊界；流體計算域的側面采用滑動墻面(slip wall)設定，這種墻面設定可以有效地提高計算速度，葉片表面為不可滑移壁面；旋轉域給定旋轉角速度，旋轉域和靜止域通過交界面連接；使用入口處流動參數來初始化流場；采用瞬態(tài)求解器。

圖3 邊界條件設定

3 結果對比

3.1 不同湍流模型計算結果對比

圖4給出了標準、RNG 、SST和Omega雷諾應力四種湍流模型計算得到的水輪機在速比為5.0時，葉片法向力系數和切向力系數隨位置角的變化歷程與試驗測量結果的比較。從總體來看，不同湍流模型計算結果和試驗結果趨勢基本一致，法向力系數要比切向力系數吻合的好。但是在數值上，不同湍流模型結果還是存在一些差異。

從圖4(a)法向力系數曲線可以看出，結果的主要差別出現在谷值和峰值上，四種湍流模型模擬結果谷值均大于試驗結果，雷諾應力模型結果偏差更大一些，其他三種模型結果相當；模擬結果峰值均小于試驗值，SST模型模擬結果與試驗值更為接近，其次是雷諾應力模型，標準模型和RNG 模型結果相當。從4(b)切向力系數曲線可以看出，四種湍流模型模擬結果在上游盤面()與試驗結果符合程度好于下游盤面時的情況，這主要是因為下游盤面經過葉片擾動，流場變得更加復雜，使模擬精度下降。在上游盤面谷值處四種湍流模型模擬結果很接近，均略大于試驗結果，同時與試驗結果相位也有偏差；上游盤面峰值處，RNG 模型結果高于試驗結果，標準模型結果低于試驗結果，SST模型和雷諾應力模型結果與試驗結果都比較接近；下游盤面谷值處四種湍流模型模擬結果很接近，但是都小于試驗值，同樣也存在相位偏差；下游盤面峰值情況，四種湍流模型模擬結果均出現了平峰，這一點與試驗值差別較大，SST模型模型結果更接近試驗值。

總的來說，采用SST模型所計算的結果與試驗結果吻合更好，這也進一步體現了SST模型對逆向壓力梯度與流動分離的流場模擬的優(yōu)勢，雖然采用雷諾應力模型計算結果與SST模型計算的準確性相差不大，但是雷諾應力模型比SST模型計算耗時更長，因此在本文以后的計算中選用SST模型。

3.2 網格尺度對計算結果的影響

網格的疏密程度與計算迭代的收斂速度以及最終的計算精度緊密相關。網格節(jié)點數量大則意味著計算結果精度較高，但這是以增加計算時間為代價的。通常情況下，網格的生成過程中人的主觀因素對網格的質量將產生很大的影響。 CFD 用戶在生成網格時僅僅憑個人經驗，在平衡計算時間和精度的前提下確定網格分布的拓撲結構、密度等。文中需要對垂直軸水輪機非定常繞流流場做大量計算，希望網格在保證一定精度的前提下，網格數盡可能少。因此在開展數值計算前，有必要對網格進行測試，測試具有不同疏密程度的網格系統(tǒng)對垂直軸水輪機計算結果的影響。

垂直軸水輪機數值模擬中，實質是葉片的繞流問題，因此葉片周圍網格最為重要，為了進行網格測試，葉片周圍生成4套密度不同的“⊙”形結構網格，周向節(jié)點布置相同，僅變化徑向的網格密度，相關信息見表2。

表2 不同網格結構相關信息

仍然以水輪機在速比5.0時的試驗結果與數值結果進行對比，以此判斷所選網格密度是否合理。圖5給出了四種不同網格密度(Mesh 1～4)計算得到葉片法向力系數和切向力系數隨位置角的變化歷程與試驗測量結果的比較。從圖5(a)法向力系數的對比曲線看出，網格疏密對法向力的影響不大，四種不同網格密度的計算結果十分接近。而網格疏密對切向力的影響比較大，從圖5(b)切向力系數的對比曲線看出，隨著網格密度的增加而呈越來越接近試驗結果的趨勢，但是Mesh2和Mesh1的結果差別不大，基本滿足網格無關性的要求。

圖5 葉片法向力和切向力系數隨位置角的變化(λ=5.0)

這也和文獻推薦值很吻合，當時可以充分利用低雷諾數模型預測分離流的優(yōu)勢。鑒于速比5.0是一個比較高的速比，分離流動還不是很明顯，因此文中后面計算全部采用網格密度類型Mesh2。

3.3 時間步長對計算結果的影響

在非穩(wěn)態(tài)問題求解時，時間步長影響模擬的收斂性及CPU計算時間。正如存在網格無關解一樣，非穩(wěn)態(tài)數值模擬應當存在時間步長無關解：數值解在工程允許的偏差范圍內不再隨時間步長的減小而發(fā)生變化。

但是，時間步長的減小導致了計算時間的迅速增加。由于描述計算對象的微分方程的剛性較強，較小的時間步長有助于得到較為準確的結果，且能提高時間步內的收斂速度；但是若計算的總時間一定，這并不足以抵消時間步數增加而引起的計算量增加。時間步長往往需要在計算精度和收斂速度之間進行折衷。下面采用簡單的增加時間步長的方法來分析時間步長對垂直軸水輪機數值模擬結果的影響，時間步長主要根據葉輪每步轉過角度來設定，采用上節(jié)中介紹的Mesh2網格模型，時間步長設置見表3。

表3 Strickland實驗輪機的主要參數

圖6給出了四個不同時間步長計算得到的水輪機在速比為5.0時，葉片法向力系數和切向力系數隨位置角的變化歷程與試驗測量結果的比較。從圖中可以看出，時間步長對計算結果的影響和網格影響很相似，都是對法向力影響比較小，對切向力影響比較大，湍流模型的影響也有類似結論，究其原因：一是因為葉片切向力系數的大小比法向力系數小一個量級；二是切向力系數的變化規(guī)律比法向力系數復雜，在一個周期內，切向力系數在上游盤面和下游盤面各出現一次峰值和谷值，即其變化頻率比法向力系數大。從圖6(b)葉片切向力系數曲線中可以看出，時間步長越大，計算結果與試驗結果的偏差越大，但是Tstep2與Tstep1的計算結果相差無幾，可以說當葉輪每個時間步轉過角度小于3°時，計算結果基本達到了時間步長無關性。

圖6 葉片法向力和切向力系數隨位置角的變化(λ=5.0)

4 結 論

通過對比以上計算結果可以得到以下結論:

(1)通過與試驗結果對比，驗證了CFD方法在垂直軸水輪機性能預報方面應用的可靠性；

(2)湍流模型選擇方面，和其他湍流模型相比，SST湍流模型計算精度更高；

(3)當滿足時基本滿足網格無關性要求；

(4)時間步長設定時，當葉輪每個時間步轉過角度小于3°時，計算結果基本達到了時間步長無關性要求。

[1] Horiuchi, K., Ushiyama, I., and Seki, K., Straight Wing Vertical Axis Wind Turbine: A Flow Analysis, Wind Engineering, Vol. 29, No.3, pp.243-252, May 2005.

[2] Guerri, O., Sakout, A. and Bouhadef, K., Simulations of the Fluid Flow around a rotating Vertical Axis Wind Turbine, Wind Engineering Vol.31, no. 3, pp 149-163, 2007.

[3] 孫 科, 豎軸H型葉輪及導流罩流體動力性能數值模擬[J].哈爾濱工程大學,2008.

[4] 姜子豪，李浩，闞江明.減振降噪在機械設計中的應用[J].森林工程，2013，29(5)：55-59.

[5] Hwang, S., Lee, Y. H., and Kim, S. J., Optimization of Cycloidal Water Turbine and the Performance Improvement by Individual Blade Control, Applied Energy, Vol. 86, Issue 9, 2009.

[6] Ansys Inc: “ANSYS CFX User Manual”, 2005

[7] Deglaire, P., Engbom, S., Agren, O., and Bernhoff, H., Analytical Solutions for a Single Blade in Vertical Axis Turbine Motion in Two-Dimensions, European J. of Mechanics B/Fluids, Volume 28, Issue 4, pp. 506-520, July-August 2009.

[8] Tchon, Ko-Foa and Paraschivoiu, I., Navier-Stokes Simulation of the Flow Around an Airfoil in Darrieus Motion, Journal of Fluids Engineering, Vol. 116, pp. 870-876., ASME, 1994.

[9] 汪魯兵.豎軸潮流水輪機水動力性能理論與實驗研究[J]. 哈爾濱工程大學,2006.