黃 揚，趙偉文，萬德成, 2

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

相較于陸上風電場，海上風電場不占用寶貴的土地資源，受到的空間限制較小，且能避免噪音污染和視覺干擾對人類活動的影響。同時，海上風能的質(zhì)量更高，風速更加平穩(wěn)，發(fā)電效率更高，因此海上風電場逐漸成為風能開發(fā)的主要趨勢[1]。近年來，隨著風電產(chǎn)業(yè)的快速發(fā)展，風力發(fā)電逐漸從近海走向遠海。隨著風機工作水深的增加，固定式海上風力機的成本急劇增大。相較而言，采用漂浮式基礎的海上風力機的性價比明顯提高，這使得浮式風機逐漸從概念設計走向工程應用[2]。為了降低風電場的建設和運行成本，海上風力機逐漸向大型化方向發(fā)展，風力機的裝機容量不斷增大，對應的風力機葉片尺寸也顯著增加[3]。大型浮式風機成為未來深遠海優(yōu)質(zhì)風能開發(fā)的重要方式和研究熱點。

浮式風機是由風力機、浮式支撐平臺以及系泊系統(tǒng)組成的復雜海洋結(jié)構(gòu)物系統(tǒng)。在風—浪—流載荷的作用下，風力機的氣動載荷和浮式支撐平臺的水動力響應之間存在顯著的相互干擾作用，使得浮式風機的氣動—水動耦合響應呈現(xiàn)出明顯的非穩(wěn)態(tài)特性[4]。此外，風力機大型化使得葉片細、長、薄的特點愈發(fā)突出。隨著葉片尺寸的增加，葉片柔性變形越來越顯著，進而會對浮式風機的耦合性能產(chǎn)生不利影響[5]。這些因素使得浮式風機推廣應用還存在較大的挑戰(zhàn)。為了加深對大型浮式風機在風浪聯(lián)合作用下耦合響應特性的認識，需要針對浮式風機的氣動—水動—氣彈性耦合響應特性開展相關(guān)研究，并進一步為浮式風機的優(yōu)化設計提供參考建議。

模型試驗是研究浮式風機耦合響應特性的可靠手段，但受限于試驗場地和設備，試驗結(jié)果不可避免地會受到尺度效應以及低雷諾數(shù)的影響，進而會影響試驗結(jié)果的適用范圍[6]。為了獲得更加可靠的浮式風機的運行數(shù)據(jù)，海外風電巨頭開展了相關(guān)的浮式風機實測項目，如挪威國家石油公司的Hywind項目[7]、Principal Power在葡萄牙的WindFloat半潛式浮式風機實測項目[8]等。盡管實測數(shù)據(jù)的價值很高，但浮式風機實測試驗的成本和風險也十分巨大。相較而言，數(shù)值方法通過對浮式風機系統(tǒng)的各個部分進行數(shù)值建模分析，可以獲得實尺度浮式風機的耦合動力響應特性，具有時間短、成本低、精度高的特點，因此被廣泛地應用于浮式風機耦合性能的研究。浮式風機氣動—水動—氣彈性耦合響應數(shù)值模擬涉及到風力機的氣動性能、浮式平臺及系泊系統(tǒng)的水動力性能以及葉片結(jié)構(gòu)變形響應之間的復雜相互作用。為了簡化計算，一些學者提出采用葉素動量理論(blade element momentum，簡稱BEM)或者渦尾跡方法(wake vortex method，簡稱WVM)來計算風力機的氣動性能，并結(jié)合Morison公式或者勢流理論預報浮式支撐平臺的水動力響應。同時，基于等效梁理論求解葉片結(jié)構(gòu)變形響應，從而建立了浮式風機氣動—水動—氣彈性耦合響應計算模型，并開發(fā)出一系列浮式風機耦合性能分析軟件，如FAST[9]、HAWC2[10]等。盡管這些耦合分析軟件具有計算快速的優(yōu)勢，但其對復雜入流工況條件下浮式風機的耦合響應問題的適用性較差。隨著的計算機技術(shù)和計算能力的快速發(fā)展，計算流體力學(computational fluid dynamics，簡稱CFD)方法因其適用范圍廣、計算精度高以及數(shù)據(jù)豐富的特點而被廣泛應用于浮式風機在復雜入流工況條件下耦合動力響應特性的數(shù)值模擬研究。Tran和Kim等[11]基于STAR-CCM+軟件和自主開發(fā)程序?qū)Π霛撌礁∈斤L機在風浪聯(lián)合作用下的耦合響應特性進行了CFD模擬；Liu等[12]等基于開源軟件平臺OpenFOAM開發(fā)了一套針對浮式風機氣動—水動耦合響應的CFD分析工具，并結(jié)合MBD(multi-body dynamics)方法求解葉片結(jié)構(gòu)變形。值得注意的是，CFD方法存在對計算資源需求大、計算速度相對較慢的特點。為此，需要進一步地對浮式風機耦合計算模型繼續(xù)完善，提高計算效率。

為了減少對計算資源的消耗，采用彈性致動線模型來模擬風力機，并基于等效梁理論和一維有限元方法計算葉片結(jié)構(gòu)變形響應，結(jié)合兩相流CFD求解器naoe-FOAM-SJTU，建立了浮式風機氣動—水動—氣彈性耦合響應計算模型，實現(xiàn)了風浪聯(lián)合作用下浮式風機耦合響應的數(shù)值模擬。通過對規(guī)則波和剪切風作用下Spar型浮式風機的氣動—水動—氣彈性耦合響應進行數(shù)值模擬，探究風力機氣動載荷、浮式平臺運動響應以及葉片結(jié)構(gòu)變形響應之間的耦合作用，并對浮式風機的尾流場特性進行了分析。

1 浮式風機氣動—水動—氣彈性耦合計算方法

1.1 風力機氣動性能計算模型

致動線模型(actuator line model，簡稱ALM)[13]是一種模擬風力機氣動性能的簡化方法。在ALM中，風力機葉片被等間距離散為一系列葉素單元，每個葉素單元被虛擬的致動點代替，然后基于葉素理論計算作用在葉素單元上的氣動載荷，并通過高斯核函數(shù)將致動點處的氣動力投影到流場，模擬風力機的尾流。ALM不需要求解葉片表面的邊界層，因此可以極大地減少了計算量。在ALM基礎上，提出了一種彈性致動線模型(elastic actuator line model，簡稱EALM)[14]，如圖1所示。該模型考慮了平臺運動以及葉片變形對風力機氣動性能的影響，實現(xiàn)了對大型浮式風機非穩(wěn)態(tài)氣動響應特性的模擬。

圖1 彈性致動線模型原理示意

為了體現(xiàn)平臺運動以及葉片變形對風力機氣動性能的影響，EALM在風力機葉片截面處的速度矢量分析中引入了由于平臺運動導致的附加速度UM和葉片變形引起的附加速度US。因此，葉片截面處的速度矢量Urel可以通過下式進行計算：

Urel=Uin+UB+UM+US

(1)

其中，Uin表示葉片截面處的入流風速，UB=Ω·r表示葉片的旋轉(zhuǎn)速度，Ω是風輪轉(zhuǎn)子的旋轉(zhuǎn)角速度，r代表葉片截面距離葉根的距離。根據(jù)葉片截面處的矢量關(guān)系，可以得到入流角φ以及葉片截面的局部攻角α：

φ=arctan(Ua/Ut)

(2)

α=φ-θtotal

(3)

θtotal=θpitch+θtwist+θtorsion

(4)

其中，Ua和Ut分別是Urel在葉片隨體坐標系下的投影，θpitch代表葉片截面所在位置的槳距角，θtwist為局部扭轉(zhuǎn)角，θtorsion為葉片扭轉(zhuǎn)變形角。根據(jù)局部攻角，結(jié)合二維翼型的氣動性能數(shù)據(jù)進行插值，可以得到葉片截面處的升力系數(shù)CL和阻力系數(shù)CD。進一步地，風力機葉片上的氣動載荷f可以通過下式進行計算：

(5)

其中，L和D分別表示作用在翼型截面處氣動升力和氣動阻力，ρ表示空氣的密度，c代表翼型截面的弦長，Nb表示風力機葉片的個數(shù)，eL和eD分別表示葉片隨體坐標下不同方向的單位矢量。

進一步地，計算得到的風力機氣動載荷需要通過高斯核函數(shù)光順地投影到葉片周圍的流場。投影到流場中的體積力可以通過下式進行計算：

(6)

(7)

其中，ηε代表高斯核函數(shù)，d表示網(wǎng)格節(jié)點與致動點的距離，ε代表高斯寬度系數(shù)，一般取值為風輪葉片附近最小網(wǎng)格尺寸的2倍或者葉片弦長的1/4。

1.2 葉片結(jié)構(gòu)變形響應計算模型

大型浮式風機的葉片具有很大的展弦比，葉片的主要變形是彎曲變形和扭轉(zhuǎn)變形，剪切變形程度較小。因此，可以采用Euler-Bernoulli梁模型對風力機葉片進行簡化處理，將之視為一端固定一段自由的懸臂梁，并選用離散參考系的多自由度模型進行計算[15]。

如圖2所示，采用一維有限元方法對風力機葉片進行離散，將連續(xù)梁結(jié)構(gòu)劃分成用有限節(jié)點相連接的離散單元體系。每個梁單元包含2個節(jié)點，每個節(jié)點處包含3個自由度，分別表示葉片的揮舞變形δx，擺振變形δy和扭轉(zhuǎn)變形δθ。根據(jù)達朗貝爾原理，可以構(gòu)建葉片結(jié)構(gòu)變形響應的動力平衡方程：

圖2 風力機葉片結(jié)構(gòu)計算模型示意

(8)

其中，M、C和K分別代表風機葉片結(jié)構(gòu)離散系統(tǒng)的整體質(zhì)量矩陣、整體阻尼矩陣和整體剛度矩陣，F(xiàn)表示作用在葉片結(jié)構(gòu)單元上的外載荷向量，y代表節(jié)點位置的全自由度向量：

(9)

其中，y1、y2和yθ分別代表全部梁單元節(jié)點沿揮舞方向的線位移向量、沿揮舞方向的線位移向量和扭轉(zhuǎn)變形的角位移向量。整體阻尼矩陣C可以通過M和K的線性組會得到：

C=a0M+a1K

(10)

其中，a0和a1是比例系數(shù)，由葉片振動的自然頻率和阻尼比進行確定：

(11)

其中，ξ代表阻尼比，fn1和fn2分別表示葉片一階振動自然頻率和二階振動自然頻率。

1.3 浮式風機氣動—水動—氣彈性耦合響應求解

采用兩相流CFD求解器naoe-FOAM-SJTU對浮式風機的水動力響應特性進行預報，并采用分段外推法求解系泊系統(tǒng)[16]。進一步地，通過結(jié)合EALM和葉片結(jié)構(gòu)變形計算模型，建立了浮式風機氣動—水動—氣彈性耦合響應數(shù)值計算模型，如圖3所示。氣動力模塊、水動力模塊以及結(jié)構(gòu)求解模塊之間通過速度、位移、力和力矩信息的傳遞，實現(xiàn)了風力機氣動載荷、平臺運動響應、系泊張力以及葉片結(jié)構(gòu)變形響應之間的耦合作用。需要說明的是，各計算模塊的可靠性已經(jīng)在之前的工作進行了驗證[5，17]，這里不再贅述。

圖3 浮式風機氣動—水動—氣彈性耦合計算模型

浮式風機氣動—水動—氣彈性耦合像響應的計算流程如圖4所示。為了體現(xiàn)各模塊之間的耦合作用，在平臺運動響應的計算中，風力機的氣動載荷和系泊張力會作為外載荷作用在平臺上，而平臺運動響應會反過來影響系泊點的位置以及風力機葉片的位置和相對風速，實現(xiàn)風力機—平臺—系泊系統(tǒng)的耦合作用。此外，風力機的氣動載荷作為外載荷會加劇葉片的結(jié)構(gòu)變形響應，而葉片變形帶來的位置和相對風速的變化也會影響到風力機的氣動載荷，從而實現(xiàn)風力機氣動響應和葉片結(jié)構(gòu)變形之間的耦合，并最終實現(xiàn)浮式風機氣動—水動—氣彈性耦合響應模擬。

圖4 浮式風機氣動—水動—氣彈性耦合響應計算流程

采用不可壓縮N-S方程描述浮式風機耦合響應數(shù)值模擬中涉及的不可壓兩相流的非定常流動，流場的控制方程：

(12)

(13)

其中，U表示速度場；ρ代表兩相流體的混合密度；Ug代表網(wǎng)格節(jié)點的速度；pd代表流場中的動壓力；g代表重力加速度；μeff代表流體的有效動力黏度系數(shù)；fσ代表表面張力的源項，僅在氣液交界面產(chǎn)生作用；fs代表消波源項，僅在消波區(qū)內(nèi)起作用；fε代表由于風力機氣動載荷投影到流場的體積力源項。此外，采用SST k-ω兩方程模型對流體控制方程進行封閉求解。

2 計算設置

2.1 浮式風機模型

選用OC3(offshore code comparison collaboration，簡稱OC3)項目中的Spar型浮式風機作為研究對象[18]。該浮式風機由NREL 5-MW風力機和Hywind Spar平臺組成，并采用懸鏈線式系泊系統(tǒng)，三根系泊纜繩均勻地分布在浮式支撐平臺的周圍。浮式風機的主要參數(shù)如表1所示[19-20]。

表1 OC3 Hywind Spar型浮式風機主要參數(shù)

2.2 計算域及網(wǎng)格布置

風浪聯(lián)合作用下浮式風機耦合響應數(shù)值模擬采用的計算域如圖5所示。計算域的長為784 m，寬為384 m，空氣的高度設置為280 m，約為2.2D(D=126 m表示的風力機的直徑)。為了減小計算量，水深設置為224 m，約為0.7d(d=320 m表示真實水深)，此時水深對平臺運動響應的影響可以忽略。風力機位于計算域的中間位置，距離入口邊界約為1.1λ(λ=146.9 m表示入射波的波長)，距離出口邊界約為5D。為了避免波浪反射對平臺運動產(chǎn)生影響，在出口邊界前200 m的區(qū)域內(nèi)設置消波區(qū)。

圖5 計算域布置情況

計算域中的網(wǎng)格布置情況如圖6所示，包括背景網(wǎng)格和加密網(wǎng)格。其中，背景網(wǎng)格沿x軸和y軸方向的尺寸均為8 m。為了減小計算域中的網(wǎng)格數(shù)量，背景網(wǎng)格的尺寸從自由液面處沿z軸正負方向線性增大，自由液面附近的背景網(wǎng)格尺寸為2 m，而計算域頂部和底部邊界處的背景網(wǎng)格尺寸增大到20 m。進一步地，為了準確捕捉到浮式風機的尾流發(fā)展情況，并減小波浪的耗散，對風力機的尾流區(qū)以及自由液面附近的背景網(wǎng)格進行加密，加密后的背景網(wǎng)格尺寸為2 m×2 m×0.5 m(自由液面)和2 m×2 m×2 m(尾流區(qū)域)。同時，對平臺附近的背景網(wǎng)格進行加密處理。計算域中的網(wǎng)格總量為559萬。

圖6 計算域中網(wǎng)格分布

2.3 風浪耦合入流條件

采用一階Stokes規(guī)則波作為入射波，參考相關(guān)文獻[18]中5級海況條件的設定，波高設置為3.66 m，波浪周期為9.7 s。此外，入流風速采用均勻入流形式，風速大小設置為11.4 m/s，對應的風輪轉(zhuǎn)速為12.1 r/min。在數(shù)值模擬過工程中，不考慮風力機的轉(zhuǎn)矩控制和變槳控制，因此風力機的轉(zhuǎn)速始終保持不變，而葉片的槳距角始終設置為零。

計算域的入口邊界采用Dirichlet邊界條件，根據(jù)入流風速和波浪類型計算不同位置處的速度，壓力梯度設置為零；計算域的出口邊界和頂部邊界均采用零壓力梯度條件；計算域的底部邊界采用滑移邊界條件；計算域的左右兩側(cè)邊界采用OpenFOAM中定義的對稱邊界條件，即垂直于邊界的方向?qū)?shù)為零。

3 算例結(jié)果分析

3.1 網(wǎng)格收斂性驗證

在開展浮式風機氣動—水動—氣彈性耦合響應數(shù)值模擬之前，首先對計算網(wǎng)格進行收斂性驗證。采用三套不同密度的網(wǎng)格開展收斂性驗證計算，具體網(wǎng)格信息如表2所示。計算風速同樣設置為11.4 m/s，采用均勻入流條件，不考慮風力機的控制系統(tǒng)。同時為了減少計算穩(wěn)定所需的時間，不放開浮式平臺的自由度。計算采用的時間步長為0.01 s，小于文獻[21]中對計算收斂所需的最小時間步，同時也滿足CFL條件。

表2 網(wǎng)格收斂性測試中的網(wǎng)格參數(shù)信息

圖7顯示了不同網(wǎng)格密度條件下風力機的氣動載荷系數(shù)隨時間變化的情況。從圖中可以看出，風力機的氣動載荷在20 s后逐漸穩(wěn)定下來。取30～50 s內(nèi)的氣動載荷數(shù)據(jù)進行統(tǒng)計分析，可以發(fā)現(xiàn)，采用中等網(wǎng)格計算得到的風力機氣動載荷系數(shù)與采用精細網(wǎng)格的計算結(jié)果之間的誤差僅為1.5%(Cp)和0.9%(Ct)，而粗糙網(wǎng)格與精細網(wǎng)格之間的誤差達到了8.8%(Cp)和4.5%(Ct)。進一步地，對采用不同網(wǎng)格密度計算得到的風力機葉尖處的結(jié)構(gòu)變形響應進行統(tǒng)計，如表3所示?？梢园l(fā)現(xiàn)，中等網(wǎng)格和精細網(wǎng)格計算得到的葉尖處變形響應的誤差在1%以內(nèi)。綜上，為了平衡計算效率和計算時間，下面將采用中等密度網(wǎng)格開展浮式風機耦合響應數(shù)值模擬。

圖7 網(wǎng)格收斂性測試中氣動載荷系數(shù)的時歷曲線

表3 網(wǎng)格收斂性測試中風力機葉尖結(jié)構(gòu)變形響應

3.2 葉片結(jié)構(gòu)變形響應

浮式風機在正常運行過程中，其葉片會在重力、離心力以及氣動載荷的作用下產(chǎn)生明顯變形，并會進一步對浮式風機的耦合性能產(chǎn)生影響。圖8顯示了浮式風機的葉片(Blade #1)在不同方向上的結(jié)構(gòu)變形響應，包括不同時刻瞬時變形情況以及結(jié)構(gòu)變形響應的時空分布云圖。

圖8 風力機葉片(Blade #1)結(jié)構(gòu)變形響應

對于葉片結(jié)構(gòu)瞬時變形圖，黑色實線表示葉片在某一時刻的結(jié)構(gòu)變形位移，虛線代表葉片瞬時變形的包絡線；對于葉片結(jié)構(gòu)變形時空分布云圖，橫坐標表示時間，縱坐標表示葉片截面沿徑向的相對位置，并采用葉片瞬時變形幅值對云圖進行染色。從圖中可以看出，浮式風機的葉片在不同方向均產(chǎn)生了明顯結(jié)構(gòu)變形。在葉片揮舞方向上，由于氣動推力的作用，葉片彎曲變形的平均值較大，由于是均勻入流條件，氣動推力的變化幅值不大，因此葉片揮舞變形的波動幅值相對較小。沿葉片擺振方向，葉片受到周期性變化的重力的影響，其擺振變形的波動幅值較大。相較于揮舞變形，由于葉片擺振方向的抗彎剛度更大，且受到的氣動載荷較小，因此葉片擺振變形的平均值明顯小于揮舞變形。此外，葉片沿軸向的扭轉(zhuǎn)變形主要受到氣動扭矩的影響，與氣動推力類似，葉片受到的氣動扭矩在均勻入流條件下的波動幅值較小，因此葉片扭轉(zhuǎn)變形的波動程度較低。值得注意的是，葉片扭轉(zhuǎn)變形會直接影響到葉片的局部攻角，進而會明顯改變風力機的氣動性能。

從葉片結(jié)構(gòu)變形響應的時空分布云圖可以看出，葉片的結(jié)構(gòu)變形響應呈現(xiàn)出明顯的周期性變化特征。其中，葉片揮舞變形的變化周期與波浪周期相近，約為9.7 s，這主要是受到周期性波動的氣動推力的影響。而葉片擺振變形的波動周期接近5 s，與風輪旋轉(zhuǎn)周期相近，這表明引起葉片擺振變形周期性波動的主要因素是葉片自身重力。從葉片扭轉(zhuǎn)變形的時空分布可以看出，其變化周期同時包含波浪周期和風輪旋轉(zhuǎn)周期的特性。進一步地，取浮式風機計算穩(wěn)定后200～300 s內(nèi)的數(shù)據(jù)進行統(tǒng)計分析，計算葉尖位置結(jié)構(gòu)變形響應的均方根(RMS)和標準差(STD)，如表4所示。RMS主要反映葉片結(jié)構(gòu)變形響應的有效值，而STD主要反映葉片結(jié)構(gòu)變形響應的波動程度。葉片揮舞變形的RMS可達3.09 m，顯著大于葉片沿擺振方向的彎曲變形。但葉片擺振變形的STD略大于葉片揮舞變形的STD，說明擺振變形的波動程度更加顯著。此外，葉片扭轉(zhuǎn)變形的RMS為3.07°，這會顯著改變?nèi)~片的局部攻角，進而影響風力機的氣動載荷。周期性變化的結(jié)構(gòu)變形響應會對葉片的疲勞壽命產(chǎn)生明顯損傷。

表4 葉尖位置結(jié)構(gòu)變形響應的均方根和標準差

3.3 風力機氣動載荷響應

由于葉片結(jié)構(gòu)變形響應與風力機氣動載荷之間的耦合作用，葉片變形會對浮式風機氣動載荷產(chǎn)生影響。為此，分別對浮式風機在考慮和不考慮葉片變形兩種情況下的耦合響應進行了數(shù)值模擬，并對比了兩種不同情況下風力機的氣動載荷系數(shù)，如圖9所示。

圖9 浮式風機氣動載荷系數(shù)的時歷曲線

圖例中“剛性葉片”表示不考慮葉片變形時的情況，“彈性葉片”表示考慮葉片變形時的情況。無論是否考慮葉片結(jié)構(gòu)變形響應，浮式風機的氣動載荷系數(shù)均呈現(xiàn)出明顯的周期性變化的特征，變化周期接近波浪周期，這主要是由于平臺運動引起的入流風速周期性變化導致的。在考慮葉片變形時，浮式風機的氣動載荷系數(shù)明顯減小。其中，氣動功率系數(shù)CP的RMS降低了約7.8%，而氣動推力系數(shù)CT的降幅更大，達到了11.8%。此外，CP的STD降低了約2.8%，CT的STD增大了約4.2%。這表明，葉片結(jié)構(gòu)變形會導致浮式風機氣動載荷平均值和波動幅值的降低。相較于氣動功率，氣動推力受到葉片結(jié)構(gòu)變形的響應更加顯著，同時，氣動載荷平均值比其波動幅值受葉片變形的影響更大。

葉片結(jié)構(gòu)變形對浮式風機氣動載荷的影響主要包括以下兩個方面：一方面，葉片彎曲變形會使得風輪迎風面積減小，從而降低了風力機受到的氣動推力，導致氣動推力的RMS降低；另一方面，葉片結(jié)構(gòu)變形顯著改變了風力機的入流條件，從而對浮式風機的氣動載荷產(chǎn)生影響。由風力機氣動載荷的計算方式可知，軸向入流風速和葉片局部攻角是影響風力機氣動載荷的關(guān)鍵因素。圖10和圖11分別對比了浮式風機的葉片(Blade #1)在是否考慮葉片變形情況下軸向入流風速Ua和局部攻角α?？梢园l(fā)現(xiàn)，在考慮葉片變形時，相對位置(x/R)在0.5～0.9范圍內(nèi)的局部葉片的入流風速有一定程度的增大，這有利于提高風力機的氣動載荷，由于扭轉(zhuǎn)變形的影響，該范圍內(nèi)葉片的局部攻角反而明顯減小，尤其在葉尖區(qū)域附近(0.7

圖10 浮式風機葉片(Blade #1)軸向入流風速(Ua)時空分布云圖

圖11 浮式風機葉片(Blade #1)攻角時空分布云圖

進一步地，葉片變形引起的氣動載荷變化會對浮式風機的葉根彎矩和偏航力矩產(chǎn)生影響。圖12對比了不同葉片變形情況下葉根切向彎矩Moop、葉根軸向彎矩Mtrq以及偏航力矩Myaw的時歷曲線。與氣動載荷的變化類似，在考慮葉片變形時，Moop、Mtrq和Myaw的RMS都有所減小，其中，Myaw的降幅最大，達到了49.8%，Mtrq的降幅最小，為7.8%，Moop的降幅為14.3%。考慮到Moop和Mtrq分別由氣動推力和氣動升力積分得到，進一步表明葉片變形對氣動推力的影響更加顯著。此外，葉片結(jié)構(gòu)變形對葉根彎矩及偏航力矩的波動幅度的影響不大。

圖12 風力機葉片(Blade #1)葉根彎矩及偏航力矩的時歷曲線

3.4 浮式平臺水動力響應

由于風力機和浮式平臺之間的耦合作用，風力機的氣動載荷會對浮式平臺及系泊系統(tǒng)的水動力響應產(chǎn)生影響。表5對作用在浮式平臺上的氣動力載荷和水動力載荷進行了總結(jié)?？梢园l(fā)現(xiàn)，氣動力載荷沿x方向的水平分量遠大于水動力載荷，這會顯著增大浮式平臺的縱蕩位移。此外，氣動力在載荷繞y軸的縱搖力矩也明顯大于水動力載荷在相應方向上的分量，這會增加浮式平臺的縱搖角。但氣動力載荷作用在平臺上的縱搖力矩的波高幅值遠小于水動力載荷，因此其對平臺縱搖運動的幅值影響較小。此外，氣動力載荷沿z軸方向的水平力和繞z軸的力矩會增大浮式平臺的垂蕩和艏搖運動幅值。進一步地，對不同葉片變形情況下浮式平臺的六自由度運動響應進行分析，如圖13所示。在考慮葉片結(jié)構(gòu)變形時，浮式平臺縱蕩、縱搖運動響應的平均值有所減小，但波動程度變化不大。其中，浮式平臺縱蕩位移的平均值降低了約12.1%，平臺縱搖角的平均值下降了約11.9%，這主要是浮式風機在葉片變形時氣動推力的降低導致的。此外，葉片變形使得浮式風機的偏航力矩有所減小，從而導致浮式平臺的艏搖運動響應的平均值和波動幅值均有所下降，但變化幅值不大。此外，由于氣動載荷在沿y軸和z軸方向的分量較小，因此葉片結(jié)構(gòu)變形引起的氣動載荷變化對浮式平臺的橫蕩、橫搖和垂蕩運動的影響較小。

表5 浮式平臺受到的氣動力載荷與水動力載荷統(tǒng)計表

平臺運動響應的變化會進一步引起系泊纜繩張力的改變。圖14顯示了不同的葉片變形下系泊纜繩的張力隨時間的變化情況。由于#2和#3系泊纜繩沿x軸方向成對稱布置，其系泊張力的大小及變化趨勢幾乎一致，因此僅對#1和#2系泊纜繩的張力進行對比分析。對比平臺縱蕩運動的時歷曲線，可以發(fā)現(xiàn)，#2系泊纜繩張力的變化趨勢與平臺縱蕩位移的變化趨勢一致，而#1系泊纜繩張力的變化趨勢與之正好相反。系泊纜繩的張力變化主要受到平臺縱蕩位移的影響，在考慮葉片變形時，平臺縱蕩位移的平均值有所減小，相應的#1和#2系泊纜繩的波動幅值分別下降了約4.1%和2.3%。

圖14 系泊纜繩張力的時歷曲線

3.5 浮式風機尾流場特性

圖15 輪轂高度平面內(nèi)時均尾流速度云圖

圖16 輪轂高度平面內(nèi)渦量場云圖

進一步地，圖17對比了不同葉片變形狀態(tài)下輪轂高度平面內(nèi)的湍動能分布云圖。在葉片所在位置，由于風速的劇烈變化，此處的湍動能較大。在風輪后方，湍動能最初主要分布在尾流區(qū)與周圍環(huán)境之間的剪切層內(nèi)。隨著尾流距離的增加，尾渦逐漸翻卷失穩(wěn)，流動不穩(wěn)定性增加，尾流中的湍動能隨之增大，湍動能的分布范圍也逐漸從剪切層向尾流中心區(qū)域發(fā)展。對比不同葉片變形情況下的湍動能分布，在考慮葉片變形時，尾流中湍動能明顯增大的位置從x/D=3附近后移到x/D=4附近，這與尾渦翻卷失穩(wěn)位置的變化是一致的，表明尾渦翻卷失穩(wěn)是引起尾流中湍動能增強的主要因素之一。此外，在葉片變形情況下，尾流場中的湍動能明顯小于不考慮葉片變形時尾流中的湍動能，表明葉片變形會使得風輪后方尾流中的湍動能有所減小。

圖17 輪轂高度平面內(nèi)湍動能分布云圖

4 結(jié) 語

基于兩相流CFD求解器naoe-FOAM-SJTU，結(jié)合彈性致動線模型和等效梁理論，建立了浮式風機氣動—水動—氣彈性耦合響應計算模型，并對均勻風和規(guī)則波作用下浮式風機的耦合響應特性開展了數(shù)值模擬，探究了葉片結(jié)構(gòu)變形響應、風力機氣動載荷、浮式平臺運動響應以及系泊張力之間的耦合作用，對并浮式風機的復雜尾流場特性進行了分析。結(jié)果表明，在氣動載荷和自身重力的作用下，風力機葉片會產(chǎn)生顯著的結(jié)構(gòu)變形，葉片揮舞變形的平均值大于葉片擺振變形，而葉片揮舞變形波動幅值則小于葉片擺振變形。受氣動載荷的影響，葉片揮舞變形呈現(xiàn)出周期性變化特性，變化周期與波浪周期相近，而葉片擺振變形主要受到自身重力的影響，其變化周期與風輪旋轉(zhuǎn)周期相近。葉片扭轉(zhuǎn)變化會顯著降低浮式風機的氣動載荷，包括氣動功率、氣動推力、葉根彎矩以及偏航力矩。相較于氣動功率，氣動推力受葉片變形的影響更加顯著。由于平臺運動響應的影響，風力機氣動載荷呈現(xiàn)出大幅度周期性波動的特征，波動周期與波浪周期相近，而在氣動載荷的影響下，浮式平臺的縱蕩位移和縱搖角會顯著增加，垂蕩和艏搖運動幅值也隨之增大。此外，葉片變形引起的風力機氣動載荷的下降，會導致浮式支撐平臺縱蕩、縱搖以及艏搖運動響應平均值的減小，而系泊張力的波動幅值也隨著平臺縱蕩位移的減小而有所下降。進一步地，葉片變形使得浮式風機在輪轂高度平面內(nèi)的尾流速度損失有所減小，并導致尾流場中尾渦失穩(wěn)的位置向后移動，同時會使得尾流中的湍動能水平有所降低。