垂蕩板對浮式風機水動—氣動耦合性能影響研究

鄭建才，趙偉文，萬德成, 2

(1. 上海交通大學 船海計算水動力學研究中心(CMHL) 船舶海洋與建筑工程學院，上海 200240; 2. 浙江大學 海洋學院，浙江 舟山 316021)

海上風場相比陸上風場而言，具有平均風速高、湍流強度低等特點，這些豐富的風資源條件為浮式風場提供了良好的經濟效益基礎。此外，浮式風電產業還能與海上石油和天然氣開發等部門形成配套設備，這使得漂浮式風力發電機具有十分廣泛的應用前景[1]。然而，由于浮式風機所處海洋環境載荷十分復雜，同時平臺受力產生大幅運動響應會增加結構系統的疲勞載荷[2]，且浮式風機作業時水動力和氣動力之間耦合效應顯著，因此如何抑制浮式風機水動力響應，提升風機風能利用率便成為了浮式風機領域的研究重點。垂蕩板作為一種安裝在浮式平臺上的被動抑制裝置，具有能增加浮式支撐平臺附加質量，進而增大固有頻率的優點，因此被廣泛使用在單立柱式Spar平臺[3]、張力腿平臺[4]以及半潛式平臺[5]。

垂蕩板不僅能夠減小浮式支撐平臺的大幅運動，還能影響浮式風機水動—氣動耦合后的氣動性能，很多國內外的相關學者都對該問題進行了研究。Tao和Dray[6]對帶垂蕩板的單柱式支撐平臺進行了水池試驗，試驗結果表明垂蕩板能夠減小原始平臺的垂蕩和縱搖運動響應幅值，這說明垂蕩板能夠有效抑制浮式平臺的運動響應。而水池試驗與數值模擬相比，代價昂貴且受限于試驗設備，不能進行垂蕩板結構參數研究，如垂蕩板安裝位置、形狀、以及開孔大小等。黃志謙等[7]采用Fluent商業軟件對不同垂蕩板模型進行數值模擬，對比分析了不同分形階數以及規則孔垂蕩板的水動力特性，研究結果表明附加質量系數隨分形階數增加而減小，阻尼力系數值隨分形階數增加而先增大后減小，且當透空率相同時，分形孔阻尼力系數較規則孔有較大提升。Moreno等[8]采用CFD方法研究了垂蕩板對浮式平臺運動響應結果的影響，并將數值模擬結果與試驗數據進行比較，研究結果表明數值模擬結果與試驗結果擬合良好，但是該數值模擬并未考慮風機氣動性能的影響。Ding等[9]采用等效力模型進行了垂蕩板抑制作用的研究，研究結果表明在等效氣動載荷的作用下，垂蕩板仍能有效抑制平臺的運動響應結果。但是等效氣動載荷的方法是一種理想化模型，且未能考慮風輪旋轉效應以及平臺與風機之間耦合效應的影響，有必要進行浮式風機氣動—水動耦合影響下垂蕩板的作用分析。Tran等[10]基于葉素動量理論結合三維勢流理論，使用FAST開源軟件進行了帶垂蕩板浮式風機耦合性能的分析，研究發現垂蕩板一方面能夠減小平臺的運動幅值，另一方面則增加了風機的輸出功率，由于勢流理論未考慮水的黏性效應，葉素動量方法也不能進行風機后方的尾流場計算，因此有必要進一步采用CFD方法進行帶垂蕩板浮式風機氣動—水動—系泊耦合性能數值模擬。

目前已經有不少浮式風機采用垂蕩板裝置的實際案例[11-12]，這表明垂蕩板對浮式風機系統具有實際的工程意義。基于開源工具箱OpenFOAM，將非穩態致動線模型嵌入兩相流CFD求解器naoeFOAM-SJTU，形成浮式風機水動—氣動耦合求解器FOWT-UALM-SJTU，并基于該求解器進行平臺水動力、風機氣動載荷和錨鏈系泊載荷三者的耦合計算，并在此基礎上分析了垂蕩板的形狀對浮式風機耦合響應的影響。

1 數值方法

1.1 非穩態致動線模型

致動線模型(actuator line model, 簡稱ALM)[13]是使用具有等效體積力的線代替風機葉片的方法，在進行流場計算時，不需要求解葉片表面邊界層，從而提升了計算效率。在致動線模型中假定風機葉片沿徑向分割為若干葉素，并使用致動線上的點代替葉片在該位置處的葉素，若已知每個葉素的弦長c和寬度dr，則該致動點處的受力可以根據該處攻角α對應的升力系數Cl、阻力系數Cd和相對于葉片的流速Urel求解得出。

在初始致動線模型中，旋轉葉片相對于葉片的流速Urel矢量關系如圖1所示。其中在進行浮式風機水動力和氣動力的耦合傳遞過程中，需要在此基礎上附加由于平臺運動引起的非穩態速度UM，形成非穩態致動線模型(unsteady actuator line model, 簡稱UALM)。

圖1中，Ω是旋轉角速度，UZ和Uθ分別是軸向和切向速度，UM,Z和UM,θ分別是平臺誘導速度在軸向和切向的分量。則該葉素處的相對速度大小Urel為：

(1)

圖1 葉片截面上速度矢量圖Fig. 1 Velocity triangle on a local blade section

若已知致動點處的α，Cl，Cd，并基于上式求得該點的相對速度Urel，可得到每個致動點上產生的力f。但是每個致動點處產生的力f為一系列離散的點力，這些點力并不能直接作用于流場中進行計算，以避免數值發散，需采用高斯權函數將反作用于流場體積力進行光順。

(2)

(3)

其中，(xi,yi,zi)是第i個致動點，di是點(x,y,z)與點(xi,yi,zi)之間的距離,ε稱為高斯光順參數。為了保證數值穩定性，文中取ε≈2Δx(Δx為葉片附近網格單元的長度)。Troldborg等[14]指出，2Δx是防止在應用空間中心差分格式時速度場出現數值振蕩的最小值。

1.2 氣動—水動耦合模型

將上述非穩態致動線模型嵌入氣液兩相流求解器naoeFOAM-SJTU，形成浮式風機氣動—水動耦合求解器FOWT-UALM-SJTU，同時考慮到大渦模擬方法(large eddy simulation, 簡稱LES)在捕捉非平衡、非穩態結構時的大尺度效應的優越性，因此在該數值模擬中采用大渦模擬方法進行浮式風機氣動—水動耦合計算，其控制方程：

(4)

式中：帶“～”符號意為空間過濾，fσ為表面張力項，fs為消波區源項，fε為體積力項。其中Smagorinsky模型被用于求解亞格子應力τSGS：

(5)

在該求解器中對于系泊力的求解使用的是準靜態方法(piecewise extrapolating method, 簡稱PEM)，與懸鏈線方程相比考慮了系泊預張力和流體力的影響。對于錨泊作用下的浮式平臺的六自由度運動響應通過naoeFOAM-SJTU求解得到，在水動力響應求解時分別采用隨平臺運動的隨體坐標系以及大地坐標系，兩者之間的關系如圖2所示。通過時間步進行PISO[15]循環求解N-S方程，更新六自由度運動引起的致動點的速度分量UM，并將求解得到的速度分量傳遞給非穩態致動線模型。與此同時受平臺運動影響后的風機氣動載荷也反饋到系泊系統作用下的平臺運動方程中，進而實現浮式風機氣動—水動耦合求解。綜上所述，使用FOWT-UALM-SJTU進行浮式風機氣動—水動耦合求解流程如圖3所示。

圖2 六自由度運動響應坐標轉換示意Fig. 2 Two coordinate systems in solving motions equation of floating platform

圖3 浮式風機氣動—水動耦合求解流程圖Fig. 3 Flow chart of solving strategy of coupled aero-hydro simulation

2 計算設定

采用的浮式風機模型由美國國家能源實驗室的NREL-5 MW[16]風力機和OC3-Hywind Spar浮式平臺組成。為了研究垂蕩板對浮式風機耦合性能的作用效果，需要確保相同的波浪載荷，同時還需要考慮風剪切效應對風機氣動性能的影響。因此在進行數值模擬時，參考Jonkman等[17]的工作，設置規則波波高為6 m，波浪周期為10 s；輪轂高度H處風速為U0=11.4 m/s，且不同高度位置Z處風速UZ服從指數規律：

(6)

分別設置三個算例進行不同形狀垂蕩板對浮式風機氣動—水動—系泊耦合性能的影響分析，其中帶圓形和正方形垂蕩板的浮式平臺模型如圖4(b)～(c)所示，三個算例采用的數值計算模型參數如表1所示，其中垂蕩板的直徑、厚度參考文獻[18]，垂蕩板的吃水位置在-130 m處。

圖4 數值計算模型Fig. 4 The model of numerical simulation

表1 算例設置

設置以上算例的計算域如圖5所示。其中風輪旋轉平面的中心距離入口處1.25D，風機后方沿流向長度為3D，計算域的寬度為3D，輪轂高度以上為1.5D，為了避免水深對平臺運動造成的干擾，沿水深方向設置為1.7d(如圖5所示，D為葉片直徑，d為平臺吃水)。同時為了進行致動線附近網格的插值計算并較好地捕捉尾流中的渦結構，在風機后方區域使用兩級網格加密，一級網格加密區域的寬為1D，高度從自由面處至輪轂高度上方1D；二級網格加密區域分別在寬度和高度方向上減小0.2D；考慮到風機后方的螺旋狀尾渦在黏性作用下的延伸及膨脹現象，加密區域一直延伸至出口邊界。在數值計算中，基于開源平臺OpenFOAM，使用自定義求解器FOWT-UALM-SJTU進行浮式風機氣動—水動全耦合數值模擬，設置入口為三維數值造波邊界條件，物面為移動壁面邊界條件，三個算例的網格總量分別為308萬、368萬和367萬。

圖5 計算域布置Fig. 5 Sketch of computational area

在求解器正確性驗證方面，Li等[19]通過使用FOWT-UALM-SJTU求解器，在相同條件下對不同俯仰運動幅值下的瞬態氣動功率和推力計算結果進行數值模擬，并與Tran等[10]采用star-CCM+商業軟件中的非定常葉素動量理論進行的多自由度數值模擬結果進行對比，結果表明該求解器能夠滿足浮式風機氣動—水動耦合問題的求解精度。

3 結果與分析

3.1 氣動載荷

觀察風力機的氣動功率時歷曲線(圖6)可知，在耦合運動影響下，功率波動周期與平臺發生的縱搖波頻運動周期相同，但帶垂蕩板浮式風機系統氣動功率達到穩定平衡的時間相對提前。這是由于垂蕩板的阻尼作用，使得平臺的縱搖運動響應穩定時間提前，而風力機的平臺運動和氣動載荷相互耦合后，風力機的輸出功率主要受縱搖運動影響。

圖6 氣動功率性能Fig. 6 The performance of aerodynamic power

此外由于垂蕩板的作用，氣動功率在初始時刻和穩定時刻幅值相應增大。這是由于在初始時刻，帶垂蕩板的浮式支撐平臺會提供一個較大的水動力阻尼力以抵消風機產生的氣動推力，這提升了浮式風機系統在深遠海抵抗湍流風的性能。而當功率趨于穩定時，此時縱蕩和垂蕩運動響應幅值降低，浮式風機處于更加穩定的動態平衡狀態，因此氣動功率會相應增加。將氣動功率的結果進行時均化處理可知，帶圓形垂蕩板的浮式風機系統功率提升約0.844%，而帶正方形垂蕩板卻使得平均功率減小1.492%。由此可見對氣動功率的提升而言，圓形垂蕩板的效果優于正方形。

為了進一步區分圓形垂蕩板與正方形垂蕩板作用效果的差異，將兩者的功率結果進行快速傅里葉變換并取模，并將橫坐標軸無量綱化為斯特勞哈爾數得到功率頻譜如圖7所示。從圖7中可以看到，圓形垂蕩板會增大氣動功率在St=0.221(fr=0.02) 外的幅值，即氣動功率波動幅度增大。

圖7 浮式風機機械功率頻譜Fig. 7 Spectrum of power output of FOWT

在非穩態致動線模型中，氣動載荷是根據葉片攻角對致動點處的升阻力系數進行插值求解，從而得到反作用于流場的體積力。在垂蕩板的作用下，浮式風機系統的波頻運動隨之改變，進而會影響風機葉片翼型攻角的變化。圖8分別給出了采用不同計算模型時風機葉尖處翼型攻角的時空云圖(考慮攻角周期性變化，圖中給出40 s內攻角隨時間和空間變化的云圖)。由圖8可知垂蕩板能夠明顯地增大0.45r～0.95r之間的翼型攻角，這有利于風能的吸收，有助于提升風機氣動功率。此外，四個波浪周期內帶圓形垂蕩板浮式風機葉尖處攻角的平均值比帶正方形浮式風機大0.002°，即在相同直徑時圓形垂蕩板增加翼型攻角的作用優于正方形，但總體相差并不大。

圖8 葉片攻角時空云圖Fig. 8 Contour of attack angle of blades

此外在進行風機的氣動載荷分析時，均方根值(root mean square，RMS)能夠代表氣動結果的有效性，表示載荷值的大小，而標準差(standard deviation，STD)能夠指出氣動結果隨時間波動的劇烈程度，值越大表示疲勞載荷越大。

對風機的軸向彎矩和切向彎矩分別進行分析如圖9所示。其中帶圓形垂蕩板浮式風機系統的軸向彎矩的均方根值相對增大0.404%，而切向彎矩的均方根值相對減小0.595%，兩者的標準差減小都在4.5%以內。而帶正方形垂蕩板軸向彎矩和切向彎矩均方根值分別減小0.595%和3.075%，兩者的標準差減小都大于30%。這說明附加圓形垂蕩板的浮式風機能在保持疲勞受力的條件下，增大浮式風機捕獲風能的效率。相對而言，雖然帶正方形垂蕩板的浮式風機系統氣動載荷波動較小，但其氣動載荷的值也降低較大，并不利于對氣動功率的提升。

圖9 風機旋轉過程中受力及誤差分析Fig. 9 Moment of turbine during rotation

3.2 水動力響應

在進行風浪聯合作用下浮式風機耦合性能研究時，運動響應與平臺固有周期、波浪頻率以及風輪旋轉周期等多種因素相關，這一方面會影響平臺的水動力性能，另一方面還會影響水動—氣動耦合后的氣動力性能。圖10給出了浮式風機耦合狀態下平臺運動響應隨時間的變化曲線。對整個波浪周期內的運動響應結果統計可知，帶圓形和正方形垂蕩板平臺的縱蕩運動響應幅值分別相應減小23.339%和6.255%。這是由于垂蕩板增大了平臺的水動阻尼力，使得系泊系統作用下的縱蕩運動幅值相應減小。帶圓形垂蕩板和正方形垂蕩板平臺的垂蕩運動響應幅值分別減小25.637%和30.524%，這主要是垂蕩板增大了浮式風機系統的質量力和阻尼力的影響，進而使得平臺的垂蕩板運動響應幅值減小。

圖10 平臺運動響應Fig. 10 Motion responses of platform

為了分析兩種形狀的垂蕩板對垂蕩運動響應的影響，對其結果進行快速傅里葉變換如圖11所示，由圖可知圓形垂蕩板能明顯減小平臺在垂蕩固有頻率處的運動響應幅值，但在波浪誘導頻率f=0.1 Hz處的幅值卻較正方形垂蕩板增大。這是由于當考慮風機和平臺的耦合作用時，風力機受到的氣動載荷會傳遞到平臺上，因此平臺不僅會受到波浪載荷的作用，還需要承受到由風力機傳遞的氣動載荷，這是兩者產生耦合后共同作用的結果。此外由于風剪切效應和偏航力矩的存在，Spar平臺存在較大的艏搖運動響應，由于垂蕩板的存在，使得帶圓形垂蕩板和正方形垂蕩板平臺的艏搖運動響應平均值分別減小2.006%和1.955%。而垂蕩板對縱搖運動響應結果并無明顯作用，在系統處于穩定平衡狀態時，帶圓形垂蕩板系統的縱搖運動響應平均值反而增大了0.439%，而帶正方形垂蕩板縱搖運動平均值減小0.072%。

圖11 垂蕩運動響應傅里葉變換結果Fig. 11 Result of heave motion with FFT transform

在進行浮式風機氣動—水動耦合計算時，浮式風機的動態平衡基于風機的氣動載荷、平臺的水動力載荷以及系泊系統的系泊載荷三者的共同作用。采用分段外推法進行系泊纜的受力分析，系泊纜張力隨時間的變化規律如圖12所示。由圖12可知，系泊纜繩的總張力變化趨勢和平臺的縱蕩運動趨勢相似，這是由于張力變化主要取決于平臺縱蕩運動響應結果。同時，由于縱蕩和縱搖運動響應之間的耦合效應，系泊張力呈現小幅度周期性波動。對比可知，初始時刻帶垂蕩板浮式風機受到的系泊張力與不帶垂蕩板浮式風機系泊張力相同，這是為了保證系泊纜在初始時刻具有相同的預張力，并在整個計算中盡量保證系泊系統的靜回復力剛度一致。此后由于浮式風機受到風浪載荷的聯合作用，開始偏離平衡位置，此時由于垂蕩板增加了浮式平臺的阻尼力，因此系泊纜產生一個較小的系泊張力即可維持系統的動態平衡位置。此外由于帶垂蕩板浮式風機系統系泊載荷波動周期隨之減小，此時應注意系泊纜疲勞載荷，避免因系泊纜疲勞載荷引起的系統失效。

圖12 系泊張力時歷曲線Fig. 12 Time history curves of mooring tension

3.3 尾流場

沿流向中縱平面的尾流速度場和渦結構組合如圖13所示。由圖可知，風浪聯合作用下入流風通過風輪后，尾流速度明顯降低，且速度虧損主要集中在風輪后方葉片中間對應的位置。同時由于風剪切效應的存在，風輪后方的不對稱性相應增加。此外隨著流向距離增大，風機后方的尾流寬度逐漸增大。

圖13 沿流向中縱平面的軸向尾流速度場與尾渦結構組合Fig. 13 Axial direction wind velocity combined with the wake vortex structure contours in the longitudinal section in center plane

在相同的風浪條件下，帶垂蕩板浮式風機的尾流速度場略微小于不帶垂蕩板的浮式風機，但三者之間的尾流速度場和渦結構總體差異不大。這說明在垂蕩板的作用下，雖然支撐平臺的水動力性能得到了改善，但風機系統的氣動功率并不一定增加。

4 結 語

使用浮式風機氣動—水動耦合求解器FOWT-UALM-SJTU，基于大渦模擬方法，針對帶垂蕩板的OC3-Hywind Spar平臺與NREL 5 MW風機組成的浮式風機系統進行了耦合性能計算，并對其數值結果及垂蕩板作用下氣動—水動耦合性能變化的原因進行了分析。此外，對比分析了相等直徑時，圓形和正方形垂蕩板的作用效果。研究結果表明：相等直徑和吃水位置時，圓形垂蕩板能使該浮式風機的平均氣動功率增大約0.844%，且圓形垂蕩板對浮式風機系統的作用效果優于正方形；帶圓形垂蕩板的浮式風機系統軸向彎矩載荷增大0.404%，切向彎矩載荷減小0.595%，且兩者的標準差增大都保持在4.5%內，這說明帶圓形垂蕩板浮式風機能在滿足疲勞載荷增加不大的條件下，增大風輪的捕獲效率；圓形垂蕩板使得縱蕩和垂蕩運動響應幅值分別減小23.339%和25.637%，正方形垂蕩板使得縱蕩和垂蕩運動響應幅值分別減小6.255%和30.524%；在穩定平衡階段，帶圓形垂蕩板浮式風機縱搖平均值增大了0.439%，而帶正方形垂蕩板浮式風機縱搖運動響應平均值減小了0.036%。

從數值計算結果可知，基于CFD技術和非穩態致動線模型進行帶垂蕩板浮式風機氣動—水動—系泊耦合性能的研究，能夠給出精細的浮式風機氣動載荷、水動力載荷和風機后方尾流場的數值模擬結果。這說明圓形垂蕩板是一個增強浮式風機系統穩定性的可行方案，且為浮式風場的數值模擬及布置優化提供了支撐。