方 圓，段 磊1,

(1.上海交通大學 海洋工程國家重點實驗室，上海 200240；2.上海交通大學 船舶海洋與建筑工程學院，上海 200240)

0 引 言

漂浮式風力機的氣動性能與其發(fā)電能力密切相關。但是，與固定式風力機不同，漂浮式風力機的氣動性能受六自由度平臺運動的影響，包括縱蕩、橫蕩、垂蕩、橫搖、縱搖和首搖。其中，縱搖運動尤為重要且特殊：1）改變流體與葉輪之間的相對速度，從而直接影響推力、扭矩及功率[1]；2）導致葉輪處于不同的俯仰狀態(tài)，對其氣動性能也有影響[2]；3）使葉輪持續(xù)進出其尾渦場，進而間接影響其氣動性能[3]。因此，本文擬對漂浮式風力機在平臺縱搖運動下的氣動性能展開研究。

截至目前，已有學者使用不同數值方法就相關內容展開研究，包括葉素-動量理論、動態(tài)尾跡模型和自由渦方法等。Jonkman 等[4]開發(fā)基于葉素-動量理論及動態(tài)尾跡模型的全耦合數值模擬工具FAST。Tran 等[5]使用上述2 種理論研究漂浮式風力機在縱搖運動下的氣動響應。Shen 等和Wen 等[6-7]分別用自由渦方法研究風力機在縱搖運動下的非定常氣動特性。然而，這些方法在本領域均有缺陷：葉素-動量理論結合不同動態(tài)入流模型使用經驗公式及修正系數，無法準確預報瞬態(tài)劇烈變化的氣動響應；動態(tài)尾跡模型基于誘導速度遠小于風速假設，限制其在低風速下的應用；自由渦方法基于勢流理論，無法準確模擬葉輪附近劇烈變化的流動現象。因此，需要更準確的數值方法對本文所述問題進行研究。

隨著計算機的發(fā)展，計算流體力學（CFD）方法越來越多地被應用于復雜流體力學問題的模擬。其中，由于適中的計算量和精度，RANS（Reynolds averaged Navier-Stokes）模型被廣泛運用。Tran 等[5]使用RANS模型研究漂浮式風力機在縱搖運動下的氣動響應。Liu 等[8]使用RANS 模型研究漂浮式風力機在縱蕩、垂蕩及縱搖耦合運動下的氣動特性。然而，RANS 模型時均化抹去重要渦結構，因此無法較好地模擬強非定常流動現象。相反，LES（large eddy simulation）模型雖然能對大尺度渦結構進行直接模擬，但是龐大的計算量制約其在工程上的應用。為兼顧計算效率和精度，本文采用由RANS 模型和LES 模型混合發(fā)展而來的IDDES（improved delayed detached eddy simulation）模型，其在近壁面使用RANS 模型，在遠場使用LES 模型，兼具RANS 模型計算量低和LES 模型計算精度高的優(yōu)點。截至目前，IDDES 已經被成功應用于垂直軸風力機[9]，但尚未被應用于本文所述問題的研究。

本文擬基于計算流體力學方法，使用IDDES 模型及重疊網格技術，對漂浮式風力機的氣動響應和周圍流場進行數值模擬，研究其氣動性能在縱搖運動影響下的特性。

1 IDDES 湍流模型

IDDES 是結合了RANS 和LES 模型的混合模型，其湍動能方程為：

式中：k 為湍動能；τij為應力張量；Sij為平均應變率；uj及xj分別為速度及位移分量。LHYBRID為IDDES 長度尺度：

其中：

式中：γ為常數0.09；CDES為模型常數；ω為比耗散率；fB與 fe用于判斷使用DDES（delayed detached eddy simulation）或WMLES（wall-modeled LES）模型進行計算，當來流有湍流成分時，激活WMLES 對邊界層進行計算。此外，fdt與ΔIDDES分別為：

式中：rdt為壁面區(qū)域標記；Δmin為該網格中心到相鄰網格中心的最小距離；Δ為網格尺寸；d為距壁面長度。

2 計算模型與網格

本文使用基于有限體積原理的STAR-CCM+軟件實施數值模擬。

2.1 幾何模型

研究對象選用Zhao 等[10]設計的1∶50 模型風力機。該模型風力機以美國國家可再生能源實驗室5 MW 參考風力機為原型，基于傅汝德數相似定律，根據推力相似使用NACA4412 翼型重新設計而來，其幾何模型如圖1 所示。

圖1 模型風力機三維模型Fig.1 The scale wind turbine model

2.2 計算域與網格

圖2 計算域及邊界條件Fig.2 The computational field and boundary setting

所有計算域形狀均為圓柱體，以適應風力機旋轉特點，如圖2 所示。最外層計算域的長度和直徑分別為9 倍和6 倍葉輪直徑，入流段和出流段分別為3 倍和6 倍葉輪直徑，其邊界條件設置為速度入口和壓力出口，四周為對稱平面以避免邊界影響。內層計算域包含2 層加密域，即網格尺寸從最外層計算域到最內層加密域以1∶2 進行2 次過渡，最內層加密域網格的基礎尺寸為0.05 m。最內層加密域到旋轉域的界面即為重疊網格界面，其兩邊網格尺寸保持一致，以保證計算精度。為避免葉輪發(fā)生較大變形，對葉片表面、輪轂和葉片隨邊進行面加密，如圖3 所示。邊界層網格的第1 層厚度為3.030 5×10-4m，增長率取1.2，以此確保葉片壁面的Y+值小于5，以符合軟件規(guī)定。

圖3 葉輪附近網格Fig.3 The mesh around the scale rotor

2.3 網格與時間無關性驗證

為保證數值模擬的可靠性，本節(jié)使用葉輪推力為參數進行網格和時間無關性驗證。

在邊界層網格不變的前提下，將其他網格的基本尺寸增大或減小21/3（對應網格體積增大或減小2 倍），形成3 種網格劃分，分別為粗糙網格、中等網格及精細網格。使用3 種網格，在風速1.61 m/s、轉速85.41 r/min（尖速比7）的工況下進行數值模擬并與實驗值[10]進行比對，如表1 所示?？芍?，粗糙網格計算結果的相對誤差較大，精細網格的網格數量較大。綜合考慮計算精度和效率，選取中等網格進行本文相關的所有數值模擬。

表1 不同網格精度對比Tab.1 Precision comprison of different mesh

選取3 個時間步長Δt1，Δt2和Δt3（對 應 每時間步葉輪旋轉1°，2°和4°），在風速1.61 m/s、轉速85.41 r/min、縱搖運動振幅2.25°、周期0.7 s 的工況下進行數值模擬并相互比對。如圖4 所示，不同時間步長下的推力十分接近，僅在峰值處存在差別。其中，Δt2的推力曲線處于Δt1和Δt3的曲線之間。同樣，綜合考慮計算精度和效率，選取Δt2進行本文相關的所有數值模擬。

圖4 不同時間步長下推力對比Fig.4 Comparison of the thrust at different time steps

2.4 風力機及平臺縱搖運動參數設置

本文相關的所有數值模擬均在風速1.61 m/s、轉速85.41 r/min（尖速比7）的工況下實施。平臺縱搖運動中心為無運動下輪轂中心垂直向下1.54 m，運動形式為簡諧運動，其位移和速度分別為：

式中：θ為縱搖運動振幅，取1.5°，2.25°和3°；T為縱搖運動周期，取0.35 s，0.7 s 和1.4 s；t為時間。

3 數值模型驗證

基于實驗數據[10]對數值模型進行驗證。其中，數值計算與實驗中的模型風力機幾何尺寸相同，風速、轉速也相同。為節(jié)省計算時間，采用移動參考系法計算葉輪推力，并與實驗值比對，如表2 所示?？芍?，相對誤差絕對值均小于6%，最小相對誤差絕對值0.51%出現在尖速比7 的工況下，即本文所有數值模擬所使用的工況，因此該數值模型可靠。

4 數值結果及分析

4.1 縱搖運動振幅對氣動性能的影響

圖5 為在相同縱搖周期（T=0.7 s），不同縱搖振幅下的葉輪推力、扭矩時域曲線。由圖可知，存在以下現象：1）推力及扭矩曲線的平衡位置不受縱搖運動振幅的影響；2）推力及扭矩輻值隨縱搖運動振幅增加而增大；3）推力和扭矩曲線關于其平衡位置不對稱，以扭矩曲線較為明顯；4）在推力和扭矩曲線峰值附近發(fā)生波動。

表2 數值模型驗證工況及結果Tab.2 Validation of the numerical model at different TSRs

圖5 不同縱搖運動振幅下的氣動響應Fig.5 The aerodynamic performance with different pitch amplitudes

上述現象可借由葉輪與風之間的相對速度解釋。漂浮式風力機在縱搖運動中，推力和扭矩與相對風速的平方成正比，相對風速由葉輪速度和自由風速合成，其中葉輪速度的幅值隨縱搖運動振幅增加而增大，因此推力和扭矩的幅值隨縱搖運動振幅增加而增大。且推力和扭矩的平衡位置大于零，上述平方關系將導致推力和扭矩曲線關于其平衡位置不對稱。

葉輪與風之間相對速度的增加可能引起失速現象，進而使推力、扭矩曲線在峰值附近發(fā)生波動。為觀察該現象，本文在典型縱搖運動（θ=3°，T=0.7 s）中，監(jiān)測距輪轂中心0.7 倍葉輪半徑處的葉片截面速度云圖，如圖6 所示。其中，圖6（a）～圖6（d）對應圖6（e）中的T1～T4點。因為設置的縱搖運動速度正方向與自由風速正方向相反，所以縱搖運動速度變化趨勢與相對速度變化趨勢一致，即圖6（a）對應的相對速度最大，并開始減小，至圖6（c）對應的相對速度最小，再開始增大。圖6（a）中葉片截面導邊和隨邊均出現流動分離，即發(fā)生動態(tài)失速現象，引發(fā)升力降低、阻力增加，從而可以解釋圖5 推力及扭矩曲線出現在峰值附近的波動。

圖6 典型工況下（θ=3°，T=0.7 s），距輪轂中心0.7 倍半徑處葉片截面速度云圖Fig.6 Velocity magnitude of one blade section at 0.7 radius from the hub center in a typical case (θ=3°，T=0.7 s)

4.2 縱搖運動周期對氣動性能的影響

圖7 為在同一縱搖振幅（θ=2.25°），不同縱搖周期下的葉輪推力、扭矩時域曲線。存在以下現象：1）推力及扭矩曲線的平衡位置不受縱搖運動周期的影響；2）推力及扭矩輻值隨縱搖周期減少而增大；3）推力和扭矩曲線關于其平衡位置不對稱，以扭矩曲線較為明顯。

上述現象也可借由相對速度解釋。如4.1 節(jié)所述，推力和扭矩與相對風速的平方成正比，相對風速由葉輪速度和自由風速合成，其中葉輪速度的幅值隨縱搖運動周期減少而增大，因此推力和扭矩的幅值隨縱搖運動周期減少而增大。且推力和扭矩的平衡位置大于零，上述平方關系將導致推力和扭矩曲線關于其平衡位置不對稱。

圖7 不同縱搖運動周期下的氣動響應Fig.7 The aerodynamic performance with different pitch periods

圖8 典型工況下（θ=3°，T=0.35 s），漂浮式風力機氣動響應時域變化Fig.8 The aerodynamic performance in a typical case (θ=3°，T=0.35 s)

圖9 典型工況下（θ=3°，T=0.35 s），一個縱搖周期內尾渦瞬時變化Fig.9 The instantaneous vorticity within one pitch period in a typical case (θ=3°，T=0.35 s)

此外，葉輪進出尾渦場可能引起尾渦干擾現象。為觀察該現象，本文提取典型縱搖運動（θ=3°，T=0.35 s）中推力和扭矩的時域曲線，如圖8 所示。其中，圖8（a）和圖8（b）中推力和扭矩曲線由上向下穿越平衡位置時發(fā)生曲率變化，圖8（b）中扭矩曲線在谷值附近發(fā)生明顯波動。該現象可以由圖9 監(jiān)測的尾渦解釋。圖9（a）～圖9（d）對應圖9（e）中的T1～T4點。同樣，因為設置的縱搖運動速度正方向與自由風速正方向相反，所以縱搖運動速度變化趨勢與相對速度變化趨勢一致，即圖9（a）對應葉輪處于平衡位置，迎風運動速度最大，葉尖渦間距與固定式風力機相近（對比圖9（f））；圖9（b）對應葉輪處于最遠離尾渦場位置，運動速度為0，即將開始順風運動，葉尖渦間距最大；圖9（c）對應葉輪處于平衡位置，順風運動速度最大，葉尖渦間距與固定式風機相近（對比圖9（g））；圖9（d）對應葉輪處于最深入尾渦場位置，運動速度為零，即將開始迎風運動，葉尖渦間距最小。上述推力、扭矩曲線在平衡位置附近的曲率變化對應圖9（e）中T2點附近，即葉輪處于最遠離尾渦場位置，運動速度為0，即將開始順風運動，同時拍擊尾渦。而扭矩曲線在谷值附近的波動對應圖9（e）的T3點附近，即葉輪處于平衡位置，順風運動速度最大，同時拍擊尾渦最劇烈。從能量角度分析，葉輪拍擊尾渦，尾渦場將部分能量傳遞給葉輪，使其推力及扭矩略增大，從而可以解釋圖8 推力及扭矩曲線出現在平衡位置附近的曲率變化以及扭矩曲線出現在谷值附近的波動。

4.3 縱搖運動對葉輪平均功率的影響

圖10 為在各種工況下的葉輪平均功率。存在以下規(guī)律：1）縱搖運動下的葉輪平均功率大于無平臺運動時的葉輪平均功率；2）相同縱搖振幅下，葉輪平均功率隨著縱搖周期的增加而減??；3）相同縱搖周期下，葉輪平均功率隨著縱搖振幅的增加而增大。上述規(guī)律也可借由相對速度解釋，同4.1 和4.2 節(jié)。

圖10 葉輪平均功率在各種縱搖運動和無運動工況下的對比Fig.10 Comparison of averaged rotor power with and without pitch motion

5 結 語

本文在STAR-CCM+軟件中使用IDDES 模型，對模型漂浮式風力機在縱搖運動下的氣動響應和周圍流場實施數值模擬。研究結果表明，漂浮式風力機氣動性能受縱搖運動影響較大，應在設計階段予以考慮，具體如下：

1）葉輪推力和扭矩的輻值隨縱搖運動振幅增加而增大。葉輪迎風運動導致相對速度增加，可能發(fā)生動態(tài)失速現象，期間升力減小、阻力增大，對風力機產生不利影響。

2）葉輪推力和扭矩的輻值隨縱搖運動周期增加而減小。葉輪順風運動導致葉輪進入尾渦場，可能發(fā)生尾渦干擾現象，期間尾渦將部分能量傳遞給葉輪，對風力機產生部分有利影響。

3）葉輪平均功率在縱搖運動下有所提高，隨縱搖運動振幅增加而增大，隨縱搖運動周期增加而減小。