漂浮式風(fēng)力機(jī)在縱搖運(yùn)動(dòng)下的氣動(dòng)性能數(shù)值研究

方 圓，段 磊1,

(1.上海交通大學(xué) 海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200240；2.上海交通大學(xué) 船舶海洋與建筑工程學(xué)院，上海 200240)

0 引 言

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

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

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

本文擬基于計(jì)算流體力學(xué)方法，使用IDDES 模型及重疊網(wǎng)格技術(shù)，對(duì)漂浮式風(fēng)力機(jī)的氣動(dòng)響應(yīng)和周圍流場(chǎng)進(jìn)行數(shù)值模擬，研究其氣動(dòng)性能在縱搖運(yùn)動(dòng)影響下的特性。

1 IDDES 湍流模型

IDDES 是結(jié)合了RANS 和LES 模型的混合模型，其湍動(dòng)能方程為：

式中：k 為湍動(dòng)能；τij為應(yīng)力張量；Sij為平均應(yīng)變率；uj及xj分別為速度及位移分量。LHYBRID為IDDES 長(zhǎng)度尺度：

其中：

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

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

2 計(jì)算模型與網(wǎng)格

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

2.1 幾何模型

研究對(duì)象選用Zhao 等[10]設(shè)計(jì)的1∶50 模型風(fēng)力機(jī)。該模型風(fēng)力機(jī)以美國(guó)國(guó)家可再生能源實(shí)驗(yàn)室5 MW 參考風(fēng)力機(jī)為原型，基于傅汝德數(shù)相似定律，根據(jù)推力相似使用NACA4412 翼型重新設(shè)計(jì)而來(lái)，其幾何模型如圖1 所示。

圖1 模型風(fēng)力機(jī)三維模型Fig.1 The scale wind turbine model

2.2 計(jì)算域與網(wǎng)格

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

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

圖3 葉輪附近網(wǎng)格Fig.3 The mesh around the scale rotor

2.3 網(wǎng)格與時(shí)間無(wú)關(guān)性驗(yàn)證

為保證數(shù)值模擬的可靠性，本節(jié)使用葉輪推力為參數(shù)進(jìn)行網(wǎng)格和時(shí)間無(wú)關(guān)性驗(yàn)證。

在邊界層網(wǎng)格不變的前提下，將其他網(wǎng)格的基本尺寸增大或減小21/3（對(duì)應(yīng)網(wǎng)格體積增大或減小2 倍），形成3 種網(wǎng)格劃分，分別為粗糙網(wǎng)格、中等網(wǎng)格及精細(xì)網(wǎng)格。使用3 種網(wǎng)格，在風(fēng)速1.61 m/s、轉(zhuǎn)速85.41 r/min（尖速比7）的工況下進(jìn)行數(shù)值模擬并與實(shí)驗(yàn)值[10]進(jìn)行比對(duì)，如表1 所示?？芍植诰W(wǎng)格計(jì)算結(jié)果的相對(duì)誤差較大，精細(xì)網(wǎng)格的網(wǎng)格數(shù)量較大。綜合考慮計(jì)算精度和效率，選取中等網(wǎng)格進(jìn)行本文相關(guān)的所有數(shù)值模擬。

表1 不同網(wǎng)格精度對(duì)比Tab.1 Precision comprison of different mesh

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

圖4 不同時(shí)間步長(zhǎng)下推力對(duì)比Fig.4 Comparison of the thrust at different time steps

2.4 風(fēng)力機(jī)及平臺(tái)縱搖運(yùn)動(dòng)參數(shù)設(shè)置

本文相關(guān)的所有數(shù)值模擬均在風(fēng)速1.61 m/s、轉(zhuǎn)速85.41 r/min（尖速比7）的工況下實(shí)施。平臺(tái)縱搖運(yùn)動(dòng)中心為無(wú)運(yùn)動(dòng)下輪轂中心垂直向下1.54 m，運(yùn)動(dòng)形式為簡(jiǎn)諧運(yùn)動(dòng)，其位移和速度分別為：

式中：θ為縱搖運(yùn)動(dòng)振幅，取1.5°，2.25°和3°；T為縱搖運(yùn)動(dòng)周期，取0.35 s，0.7 s 和1.4 s；t為時(shí)間。

3 數(shù)值模型驗(yàn)證

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

4 數(shù)值結(jié)果及分析

4.1 縱搖運(yùn)動(dòng)振幅對(duì)氣動(dòng)性能的影響

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

表2 數(shù)值模型驗(yàn)證工況及結(jié)果Tab.2 Validation of the numerical model at different TSRs

圖5 不同縱搖運(yùn)動(dòng)振幅下的氣動(dòng)響應(yīng)Fig.5 The aerodynamic performance with different pitch amplitudes

上述現(xiàn)象可借由葉輪與風(fēng)之間的相對(duì)速度解釋。漂浮式風(fēng)力機(jī)在縱搖運(yùn)動(dòng)中，推力和扭矩與相對(duì)風(fēng)速的平方成正比，相對(duì)風(fēng)速由葉輪速度和自由風(fēng)速合成，其中葉輪速度的幅值隨縱搖運(yùn)動(dòng)振幅增加而增大，因此推力和扭矩的幅值隨縱搖運(yùn)動(dòng)振幅增加而增大。且推力和扭矩的平衡位置大于零，上述平方關(guān)系將導(dǎo)致推力和扭矩曲線關(guān)于其平衡位置不對(duì)稱。

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

圖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 縱搖運(yùn)動(dòng)周期對(duì)氣動(dòng)性能的影響

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

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

圖7 不同縱搖運(yùn)動(dòng)周期下的氣動(dòng)響應(yīng)Fig.7 The aerodynamic performance with different pitch periods

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

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

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

4.3 縱搖運(yùn)動(dòng)對(duì)葉輪平均功率的影響

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

圖10 葉輪平均功率在各種縱搖運(yùn)動(dòng)和無(wú)運(yùn)動(dòng)工況下的對(duì)比Fig.10 Comparison of averaged rotor power with and without pitch motion

5 結(jié) 語(yǔ)

本文在STAR-CCM+軟件中使用IDDES 模型，對(duì)模型漂浮式風(fēng)力機(jī)在縱搖運(yùn)動(dòng)下的氣動(dòng)響應(yīng)和周圍流場(chǎng)實(shí)施數(shù)值模擬。研究結(jié)果表明，漂浮式風(fēng)力機(jī)氣動(dòng)性能受縱搖運(yùn)動(dòng)影響較大，應(yīng)在設(shè)計(jì)階段予以考慮，具體如下：

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

2）葉輪推力和扭矩的輻值隨縱搖運(yùn)動(dòng)周期增加而減小。葉輪順風(fēng)運(yùn)動(dòng)導(dǎo)致葉輪進(jìn)入尾渦場(chǎng)，可能發(fā)生尾渦干擾現(xiàn)象，期間尾渦將部分能量傳遞給葉輪，對(duì)風(fēng)力機(jī)產(chǎn)生部分有利影響。

3）葉輪平均功率在縱搖運(yùn)動(dòng)下有所提高，隨縱搖運(yùn)動(dòng)振幅增加而增大，隨縱搖運(yùn)動(dòng)周期增加而減小。