風(fēng)機尾流的CFD數(shù)值模擬

孫 毅， 吳鳳波， 程 旭， 李 藝， 李天成

(西南交通大學(xué)土木工程學(xué)院，四川成都 610031)

風(fēng)能是一種重要的清潔能源，隨著人們保護(hù)環(huán)境意識的逐步提高，對于這種可再生能源的開發(fā)更加迫切，我國風(fēng)力發(fā)電量在能源供給中占了重要地位[1]。隨著風(fēng)力發(fā)電場在全球范圍的快速發(fā)展，風(fēng)機尾流和大氣湍流的相互作用對于風(fēng)力發(fā)電的效率有很大影響，并且這種相互作用對于尾流區(qū)的風(fēng)機產(chǎn)生巨大的動力荷載，對于風(fēng)機的使用壽命也有很大的影響[2-5]；同時，風(fēng)在通過風(fēng)機的時候會產(chǎn)生動能的損失，這種損失的動能轉(zhuǎn)化為電能，這就造成風(fēng)機后方的風(fēng)的速度減小，為了保證后方風(fēng)機的發(fā)電效率，對于風(fēng)機尾流的研究就顯得至關(guān)重要[6]。

隨著計算機技術(shù)的快速發(fā)展，計算流體力學(xué)也成為一大熱門[10]，相比傳統(tǒng)的風(fēng)洞試驗，其具有準(zhǔn)確、方便、迅速等特點，因此，數(shù)值模擬在風(fēng)機的尾流的研究這方面具有巨大的潛力[7-8]。也有些學(xué)者采用雷諾應(yīng)力模型(Reynolds Stress Model ,RSM)對風(fēng)機尾流進(jìn)行了數(shù)值模擬[9]，本文使用另外一種湍流模型(LES)對風(fēng)機進(jìn)行數(shù)值模擬。

1 數(shù)值計算

1.1 控制方程

在大渦數(shù)值模擬當(dāng)中，大尺度渦可以被直接計算，當(dāng)渦的影響小于設(shè)置的網(wǎng)格間距時，采用Boussinesy假設(shè)和動力亞格子模型(Smagorinsky-Lill model)來計算亞格子尺度的雷諾應(yīng)力(The subgrid-scale Renolds Stress)。

控制方程被過濾成不可壓縮的N-S方程，Smagorinsky-Lilly模型通常被用于SGS湍流粘度的計算：

(1)

(2)

μt代表SGS湍流粘度，δij是克羅內(nèi)則符號，在式中LS代表亞格子模型的混合尺度，并且可以使用下面的公式對LS進(jìn)行計算：

LS=min(κδ,CsV1/3)

(3)

κ代表卡門常數(shù)，取值為0.42，δ是到近壁面的距離，V是計算單元的體積。在本文中Smagrinsky常數(shù)Cs的取值為0.032。

(4)

(5)

式中：μ是垂直于壁面的過濾速度，y是網(wǎng)格中心到壁面的距離，μτ代表摩擦速度，常數(shù)E的取值為9.793。

本文在ANSYS FLUENT的平臺上采用有限體積法(Finite Volume Method)對風(fēng)機尾流進(jìn)行數(shù)值模擬，二階中心差分格式被用于離散對流項和粘度項，二階隱式格式被用于計算非穩(wěn)定項。使用SIMPLE(Semi-implicit pressure linked equations)算法求解差分方程。本文中的數(shù)值模擬一共計算了80 s，最后50 s用于采樣以求得變量的時間平均。

1.2 葉素動量理論

假設(shè)風(fēng)流過風(fēng)機的過程相當(dāng)于空氣流過流管的過程，并對此過程作以下假設(shè)：(1)氣流均勻、不可壓縮；(2)氣流無摩擦；(3)具有無限多的葉片；(4)推力均勻作用在轉(zhuǎn)子葉輪上；(5)尾流無旋轉(zhuǎn)。在此過程中考慮空氣的質(zhì)量守恒、動量守恒和能量守恒(圖1)。

ρA0U0=ρADUD=ρAWUW

(6)

(7)

(8)

(a) 葉片微元

(b) 每個葉片單元的受力和速度圖1 葉片構(gòu)造

下標(biāo)為0表示上游變量，下標(biāo)為D表示風(fēng)機位置變量，下標(biāo)為W表示下游變量，P代表壓力。并引入軸向誘導(dǎo)因子a：

(9)

因為對于風(fēng)機前的軸向速度UD=U0(1-a)，風(fēng)機后的軸向速度UW=U0(1-2a)即在通過風(fēng)機前減少aU0，通過風(fēng)機后減少aU0。氣流在通過風(fēng)機之前，不帶有角速度，通過風(fēng)機后會帶有角速度，所以風(fēng)機對于氣流的角速度起到增加的作用，類似于軸向誘導(dǎo)因子a速度損失的貢獻(xiàn)，那么切向誘導(dǎo)因子a′對于氣流的角速度增量:

(10)

根據(jù)一維動量理論可以得出作用于每一環(huán)空氣微元上的推力和轉(zhuǎn)矩：

(11)

(12)

ADM能很好地模擬風(fēng)機在風(fēng)洞中的情況，并且根據(jù)這種方法能很好求解風(fēng)機的軸向誘導(dǎo)因子和切向誘導(dǎo)因子，并運用葉素理論求解作用于風(fēng)機葉片上的升力和阻力。葉素理論中，圖1所示，每個葉片被分成N段，每段考慮其受到空氣運動所產(chǎn)生的升力和阻力的影響，并且每個葉素之間沒有相互影響。作用于每個葉素上風(fēng)機上升力和阻力可以表示為：

(13)

(14)

c表示葉素弦長，CD(α)表示阻力系數(shù)，CL(α)表示升力系數(shù)，dr表示葉素長度，Urel表示相對于每個葉素的來流速度，ρ表示空氣密度。對于每個葉素上的軸向力和切向力可以表示為：

dFx=dFLcosφ+dFDsinφ

(15)

dFθ=dFLsinφ-dFDcosφ

(16)

現(xiàn)在作用于每個葉素上的力轉(zhuǎn)化成體積微元Δx·dA上的體積力，Δx代表體積微元的厚度，微元面積dA=2πr·dr，體積力可以表示為：

(17)

(18)

其中B是葉片數(shù)量。

1.3 數(shù)值計算

采用1∶100的縮尺模型對風(fēng)機流場進(jìn)行數(shù)值模擬，計算域采用長寬高分別為13 m、1.5 m、1.8 m的長方體，在距離入口5.5 m中心處放置風(fēng)機，風(fēng)機直徑D為0.4 m，輪轂高度H為0.7 m。在入口采用速度入口的邊界條件，在出口采用壓力出口的邊界條件，兩側(cè)和頂部采用對稱邊界條件，底面采用無滑移壁面邊界條件。將風(fēng)機的旋轉(zhuǎn)中心，發(fā)動機艙和支撐塔架使用尺寸為0.01 m的四面體非結(jié)構(gòu)化網(wǎng)格進(jìn)行劃分。在風(fēng)機網(wǎng)格和其它網(wǎng)格區(qū)域之間采用0.02 m的四面體非結(jié)構(gòu)化網(wǎng)格進(jìn)行過渡。在計算域的兩頭采用結(jié)構(gòu)化網(wǎng)格(圖2)。

(a) 風(fēng)機網(wǎng)格

(b)計算域

1.4 入口來流

速度入口的選擇見圖3。

數(shù)值計算中的大氣邊界層采用域前模擬的方式生成，結(jié)果證明這種方法能很好的模擬速度入口的來流風(fēng)剖面的平均速度和湍流度。風(fēng)機輪轂高度H處的平均風(fēng)速Uh為10 m/s，本文所有的平均速度風(fēng)剖面都通過這個速度進(jìn)行無量綱化。湍流度被定義為：

(19)

式中：σu表示x方向的速度標(biāo)準(zhǔn)差。

(a) 速度風(fēng)剖面

(b)湍流度圖3 來流速度風(fēng)剖面和湍流度

1.5 數(shù)值模擬結(jié)果

數(shù)值模擬的結(jié)果見圖4，圖4是由通過風(fēng)機中心x-z平面的平均速度云圖和通過風(fēng)機中心距離原地分別為2D、4D、6D、8D位置的尾流區(qū)域速度風(fēng)剖面，展示了尾流區(qū)域不同位置的平均速度大小。實線是大渦模擬的結(jié)果，圓圈是風(fēng)洞實驗結(jié)果，圖中計算曲線很好的同實驗結(jié)果相吻合，證明大渦模擬的有效性。圖中所有的速度都是被Uh無量綱化的量。

(a) 尾流速度

(b)速度云圖圖4 尾流速度剖面和云圖

2 結(jié)論

本文通過大渦模擬實現(xiàn)了對于風(fēng)機尾流的數(shù)值模擬，通過數(shù)值模擬對風(fēng)機尾流的平均速度做出了進(jìn)一步的研究，數(shù)值模擬結(jié)果能很好地同實驗數(shù)據(jù)吻合，從結(jié)果可以得出以下結(jié)論：

(1)基于葉素動量理論的制動盤模型的大渦模擬能夠很準(zhǔn)確地實現(xiàn)風(fēng)機尾流的數(shù)值模擬，同實驗結(jié)果相比，數(shù)據(jù)基本吻合。在一定的條件下甚至可以取代風(fēng)洞試驗，節(jié)約成本，提高經(jīng)濟(jì)效益。

(2)隨著尾流區(qū)域到風(fēng)機中心的距離增大，風(fēng)機尾流的平均速度越來越大，速度受影響的范圍越來越寬，影響峰值愈來愈小。

(3)在風(fēng)機輪轂高度的水平面上，速度最小點產(chǎn)生在風(fēng)機中心，隨著到風(fēng)機中心的展向距離的增大，速度愈來愈大。