基于OpenFOAM的層結性對風力機尾流大氣參數(shù)的影響

吳正人，張智博，劉 梅，杜思源，張 雷

（1.華北電力大學 動力工程系，河北 保定 071003；2.華北電力大學經(jīng)濟管理系，河北 保定 071003）

0 引言

隨著全球經(jīng)濟的快速發(fā)展，優(yōu)化能源結構已成為可持續(xù)發(fā)展的重要組成部分。風能為優(yōu)化全球整體能源結構和滿足世界整體能源供應做出了巨大貢獻。風力機作為開發(fā)風能的發(fā)電裝置，其運行必定會破壞大氣層的物質(zhì)能量循環(huán)，因此研究風力機運行對大氣參數(shù)的影響具有重要的意義。袁仁育[1]采用動態(tài)入流條件對風力機運行進行數(shù)值模擬，開發(fā)多尺度下的風力機模型，發(fā)現(xiàn)風電場會對風力機尾流10 km內(nèi)的大氣邊界層產(chǎn)生影響。剛蕾[2]研究了大氣穩(wěn)定性對風力機尾流的影響，其提出的Park-Gauss模型對尾流區(qū)風速模擬取得了較好的效果。Xie S[3]通過大渦模擬分析不同溫度層結下的風力機尾流，發(fā)現(xiàn)尾流的風速、濕度和溫度變化受溫度層結影響較大。Sharma V[4]利用大渦模擬研究發(fā)現(xiàn)，區(qū)域大氣層在風電場影響下升高了200 m。徐榮會[5]通過研究蘇尼特右旗某一風電場區(qū)域的氣候變化情況，發(fā)現(xiàn)風速和陸地地表溫度會隨著風向變化逐漸下降，風力機尾流區(qū)域的近地層會產(chǎn)生降溫增濕的現(xiàn)象。通過對蘇格蘭Black Law風電場的氣象數(shù)據(jù)進行分析，Armstrong A[6]發(fā)現(xiàn)風力機周圍環(huán)境溫度升高了0.18℃、濕度升高了0.03 g/m3。Lu H[7]通過大渦模擬與風力機致動線技術結合的方法，發(fā)現(xiàn)科里奧利效應和葉片運動引起的旋轉效應會導致邊界層高度增加，其表面動量通量減小了30%，表面浮力通量減小了15%。Emanuel K[8]通過全球氣候模型設想在全球各地部署足夠的風力機來開發(fā)整個風能資源的極端情況，發(fā)現(xiàn)平均表面溫度可能變化幾攝氏度，兩極可能降低多達10℃。張亞光[9]通過模擬3種不同的大氣穩(wěn)定度，發(fā)現(xiàn)隨著穩(wěn)定度的增加，尾流湍流強度和剪切應力下降。楊祥生[10]通過新修正的Park-Gauss和Park-Polynomial模型來探究大氣邊界層穩(wěn)定性對風電場微觀選址的影響，發(fā)現(xiàn)大氣邊界層越不穩(wěn)定，越能促進風力機尾流的恢復。

本文采用開源CFD平臺OpenFoam，實現(xiàn)對PimpleDyMFoam求解器的二次開發(fā)，并結合大氣能量守恒方程編譯求解溫度場的分布，采用codedFixedValue邊界條件建立大氣邊界層的風廓線模型和大氣溫度層結模型，從而建立不同溫度層結下的風力機模型，并在此基礎上研究風力機運行對大氣邊界層參數(shù)的影響。

1 模型建立及數(shù)值模擬方法

1.1 控制方程組

不可壓縮流動大渦模擬的連續(xù)方程和N-S方程可表示為

式中：υθ為熱擴散系數(shù)；Rj為xj方向的凈輻射分量；cp為濕空氣的定壓比熱；E為單位時間相變的水量。

方程左邊第一項表示能量的貯存，第二項表示平流作用。方程右邊第一項表示分子熱傳導作用，第二項表示輻射散度，描述xj方向凈輻射通量密度的梯度引起的空氣層溫度的變化，第三項表示因相變水汽凝結或蒸發(fā)時所產(chǎn)生的能量變化。

本文將大氣作為不可壓縮均質(zhì)流體處理，同時不考慮氣壓的變化對大氣溫度的影響，所以采用溫度T代替位溫θ，根據(jù)大氣能量守恒方程簡化的溫度方程為

1.2 幾何模型及網(wǎng)格劃分

以1.2 MW水平軸風力機為研究對象，應用坐標變換理論得到葉片截面的坐標數(shù)據(jù)，從而建立風力機葉片實體模型，進而完成風力機流場的建模工作。風力機參數(shù)如表1所示[12]。

表1 風輪及風況參數(shù)Table 1 Wind wheel and wind condition parameters

風力機運行時承載負荷的主要區(qū)域在輪轂和槳葉連接處，所以選用NACA634翼型，而葉尖為了提升葉片的氣動特性選用FX66S196翼型。風力機d為60 m，槳葉安裝角和風輪錐角均為0°，輪轂直徑為2.8 m，轉軸傾角為4°。通過Gambit軟件建立的葉片實體模型如圖1所示。

圖1 葉片實體模型圖Fig.1 Solid model diagram of blade

流場的計算域尺寸設置如圖2所示。

圖2 整體流場區(qū)域尺寸Fig.2 Overall flow field area size

流向沿Z軸，風輪中心高度為60 m。為了使入口風得到充分發(fā)展，把風力機設置在流場入口后3d的位置處。葉片區(qū)域網(wǎng)格劃分如圖3所示。

圖3 葉片區(qū)域網(wǎng)格劃分Fig.3 Grid division of blade area

通過比較模型中心線上湍流動能變化曲線來驗證網(wǎng)格無關性（圖4）。對比網(wǎng)格數(shù)分別為595萬、652萬和724萬的模型，發(fā)現(xiàn)數(shù)目較大的兩個模型的計算結果接近，而數(shù)目較少的網(wǎng)格模型計算結果有所偏差，所以選用652萬的網(wǎng)格模型進行計算。

圖4 網(wǎng)格無關性驗證Fig.4 Grid independence verification

1.3 模擬條件設置

為了更好地適應模擬的實際工況，本文采用codedFixedValue邊界條件對入口風速的垂直分布以及入口的溫度進行編譯。入口壓力設置為zeroGradient邊界條件；對于流場區(qū)域的出口，壓力采用toalPressure邊界條件，速度采用pressureInletOutletVelocity邊界條件；在風力機所處的旋流區(qū)域，因為與流場的其他區(qū)域進行物質(zhì)交換，采用cyclicAMI邊界條件；對于流場區(qū)域的頂部和兩側采用symmetry邊界條件，而底部需采用無滑移邊界條件。

風力機的來流風速由風廓線方程定義。

式中：u（zr）為參考風速，u（zr）=10.0；zr為參考高度，zr=10.0；α為指數(shù)冪，α=0.12。

入口風溫為305 K，為了加快收斂速度防止發(fā)散，采用adjustableRunTime，即變時間步長的控制方法，時間步長設為10-4s，所有參數(shù)的殘差設置為10-4。在LES下的S-A湍流模型基礎上，對時間項采用歐拉離散，梯度項、散度項和拉普拉斯項均采用高斯方法，線性插值為二階離散。風力機周圍小區(qū)域內(nèi)流體隨著葉片及輪轂共同以19.27 r/min的轉速勻速旋轉，其余外部區(qū)域的流場靜止。

1.4 模擬方法驗證

將風力機尾流z=7d以及z=10d處的速度剖面和瑞典航空研究院（FFA）風洞實驗數(shù)據(jù)結果[13]進行對比（圖5），該風洞采用L1回流式風洞，長、寬、高分別為6，2.4 m和1.4 m，風洞自由來流的湍流水平為0.5%，風速為8 m/s；風力機模型直徑D為0.25 m，輪轂高度為1D。

由圖5可知：數(shù)值模擬的數(shù)據(jù)與風洞實驗的數(shù)據(jù)在分布規(guī)律上具有一定的相似性，風力機的近尾流區(qū)速度虧損明顯；風洞實驗數(shù)據(jù)較數(shù)值模擬數(shù)據(jù)偏小。

圖5 模擬結果與風洞實驗結果尾流速度剖面比較Fig.5 Comparison of wake velocity profile between simulation results and wind tunnel test results

2 計算結果與分析

2.1 速度場的對比分析

對風力機尾流影響最大的因素是當?shù)氐奶鞖鈼l件和風電場入口自由風的狀態(tài)，自由風被風力機吸入后動能減小，使得尾流區(qū)原本正常的湍流輸運和交換循環(huán)發(fā)生改變。風力機輪轂高度處的速度分布如圖6所示。

圖6 輪轂高度處的速度分布Fig.6 Velocity distribution at hub height

由圖6可知，兩種溫度層結狀態(tài)下的速度分布規(guī)律基本一致，但穩(wěn)定層結下尾流對區(qū)域大氣邊界層的影響范圍明顯比不穩(wěn)定層結遠。這是因為在不穩(wěn)定層結下，高處空氣的溫度低、密度大，從而有空氣向下流動及流速增大的現(xiàn)象。另外尾流會與不穩(wěn)定邊界層中的大尺度旋渦產(chǎn)生劇烈的相互作用，使得尾流剪切層的形狀被快速破壞，且空間分布更加均勻，速度的補充也較快。

在不穩(wěn)定層結下，尾流在3d后呈現(xiàn)蜿蜒特性，而在穩(wěn)定層結下，尾流的形狀相對固定，這反映了穩(wěn)定層結對剪切速度的抑制作用。

圖7為尾流區(qū)不同距離處風速的垂直分布曲線。圖中向下凸起的范圍即為尾流的高剪切層，也是尾渦結構的主要部分。

圖7 尾流區(qū)不同距離風速的垂直分布Fig.7 Vertical distribution of wind speed at different distances in the wake region

由圖7可知：在穩(wěn)定層結下，大氣的速度損失最大為7 m/s，而不穩(wěn)定層結下為5 m/s；穩(wěn)定層結和不穩(wěn)定層結下風速的變化規(guī)律基本相似。不穩(wěn)定層結下的速度損失較小是由于自由風的風速較大，動能交換劇烈。觀察不穩(wěn)定層結z=1d處可以發(fā)現(xiàn)，葉片尖端處的風速比輪轂處要小，這正好描述了尾流中高剪切層是呈圓柱形分布的，而穩(wěn)定層結中該現(xiàn)象并不明顯，體現(xiàn)了穩(wěn)定層結的自由風分布較穩(wěn)定，并對剪切速度有明顯的抑制作用。其它曲線由于逐漸遠離風力機，損失的動能逐漸從周圍空氣中得到補充，所以不斷貼近入口風速的分布曲線。

2.2 溫度場的對比分析

溫度作為研究大氣邊界層的關鍵參數(shù)，其影響程度雖不及風速，但仍值得考慮。對比風廓線方程和溫廓線方程可以看出，風廓線的變化率要遠大于溫度，使得局部區(qū)域的溫度分布差異很小。同時伴隨著尾渦結構的破裂，對流換熱的作用在遠離風力機的位置也逐漸減弱，但并不影響從中探討溫度參數(shù)的變化。圖8為輪轂高度處的溫度分布圖。

圖8 輪轂高度處的溫度分布Fig.8 Temperature distribution at hub height

由圖8可知：在遠尾流區(qū)，不穩(wěn)定層結的溫度較來流溫度的變化比穩(wěn)定層結明顯；在不穩(wěn)定層結下，遠尾流區(qū)域一側溫度較來流溫度下降約為0.7 K，這是因為由于風輪的轉動，高處的低溫空氣與低處的高溫空氣進行混合，出現(xiàn)了風力機一側溫度較低的現(xiàn)象；而穩(wěn)定層結下的溫度變化規(guī)律與不穩(wěn)定層結相反，其溫度升高約為0.15 K；在遠尾流區(qū)域，風輪輪轂高度處的溫度較來流溫度的變化程度均會隨著距離的增大而增大。

圖9為尾流區(qū)不同距離溫度的垂直分布曲線。

圖9 尾流區(qū)不同距離溫度的垂直分布Fig.9 Vertical distribution of temperature at different distances in the wake region

由圖9可知：在高度為1.8d附近，溫度有明顯的界限點，高于此值的尾流區(qū)溫度變大，而且作用的范圍超過1d以上，并且隨著尾流的發(fā)展，其溫度升高的幅度減小；在1.8d高度以下的尾流區(qū)，溫度較低，這說明風力機的運行會使其處于降溫的環(huán)境中；在高度為0.5d以下時，尾流區(qū)各溫度曲線的變化相似，說明在此區(qū)域的溫度分布不受尾流發(fā)展的影響，在高度為3d以上時也有類似規(guī)律。

在穩(wěn)定層結下，高度為0.5d以下區(qū)域的溫度升高，而高度為0.5d以上區(qū)域的溫度則下降，并且溫度升高的幅度小于溫度降低的幅度，這反映了風輪的旋轉會使高溫空氣向下擴散，并與下墊面的低溫空氣混合，使得下墊面出現(xiàn)溫度升高的現(xiàn)象。在尾流1d處降幅達到0.15 K，5d處降幅達到0.075 K，在高度為3.5d以上的區(qū)域溫度已不發(fā)生明顯變化。在兩種層結狀態(tài)下，隨著尾流的發(fā)展，曲線彎曲程度變小，說明尾渦結構的不斷分解使其作用的范圍也在變小。

上述分析印證了文獻[6]中得到的風電場下游區(qū)域的溫度比上游區(qū)域的溫度明顯升高的現(xiàn)象，并由此得出夜晚風電場會顯著影響尾流區(qū)下墊面空氣溫度的結論。

3 結論

本文采用OpenFoam對1.2 MW風力機進行了模擬分析，得到以下結論。

①通過開發(fā)編譯特定環(huán)境下的OpenFoam求解器，以實現(xiàn)不同層結狀態(tài)下的風力機模型，與固定求解器相比，其模擬方式有較好的靈活性和針對性。

②在穩(wěn)定層結下，風力機尾流對大氣邊界層的影響范圍比不穩(wěn)定層結遠，不穩(wěn)定層結下的風速損失比穩(wěn)定層結小。但兩種層結下風力機尾流的風速變化趨勢基本相似，隨著尾流距離的增加其變化的趨勢變小，且速度逐漸向入口風速恢復。

③在不穩(wěn)定層結下，遠尾流區(qū)溫度較來流溫度的變化比穩(wěn)定層結明顯。在不穩(wěn)定層結下，輪轂高度處遠尾流區(qū)的溫度會下降，而穩(wěn)定層結下的溫度則會上升，并且這種變化趨勢會隨著距離的增加不斷加大。

④隨著高度的增加，不穩(wěn)定層結下的風力機尾流溫度呈先下降再上升后下降趨勢，而穩(wěn)定層結下的尾流溫度變化趨勢與其相反；在不穩(wěn)定層結下，風輪高度以下區(qū)域的尾流溫度分布不受尾流發(fā)展特性變化的影響，且高度為3d以上的區(qū)域有相似的規(guī)律。