黃 揚(yáng)，趙偉文，萬(wàn)德成, 2

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

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

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

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

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

1 浮式風(fēng)機(jī)氣動(dòng)—水動(dòng)—?dú)鈴椥择詈嫌?jì)算方法

1.1 風(fēng)力機(jī)氣動(dòng)性能計(jì)算模型

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

圖1 彈性致動(dòng)線模型原理示意

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

Urel=Uin+UB+UM+US

(1)

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

φ=arctan(Ua/Ut)

(2)

α=φ-θtotal

(3)

θtotal=θpitch+θtwist+θtorsion

(4)

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

(5)

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

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

(6)

(7)

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

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

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

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

圖2 風(fēng)力機(jī)葉片結(jié)構(gòu)計(jì)算模型示意

(8)

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

(9)

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

C=a0M+a1K

(10)

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

(11)

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

1.3 浮式風(fēng)機(jī)氣動(dòng)—水動(dòng)—?dú)鈴椥择詈享憫?yīng)求解

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

圖3 浮式風(fēng)機(jī)氣動(dòng)—水動(dòng)—?dú)鈴椥择詈嫌?jì)算模型

浮式風(fēng)機(jī)氣動(dòng)—水動(dòng)—?dú)鈴椥择詈舷耥憫?yīng)的計(jì)算流程如圖4所示。為了體現(xiàn)各模塊之間的耦合作用，在平臺(tái)運(yùn)動(dòng)響應(yīng)的計(jì)算中，風(fēng)力機(jī)的氣動(dòng)載荷和系泊張力會(huì)作為外載荷作用在平臺(tái)上，而平臺(tái)運(yùn)動(dòng)響應(yīng)會(huì)反過來(lái)影響系泊點(diǎn)的位置以及風(fēng)力機(jī)葉片的位置和相對(duì)風(fēng)速，實(shí)現(xiàn)風(fēng)力機(jī)—平臺(tái)—系泊系統(tǒng)的耦合作用。此外，風(fēng)力機(jī)的氣動(dòng)載荷作為外載荷會(huì)加劇葉片的結(jié)構(gòu)變形響應(yīng)，而葉片變形帶來(lái)的位置和相對(duì)風(fēng)速的變化也會(huì)影響到風(fēng)力機(jī)的氣動(dòng)載荷，從而實(shí)現(xiàn)風(fēng)力機(jī)氣動(dòng)響應(yīng)和葉片結(jié)構(gòu)變形之間的耦合，并最終實(shí)現(xiàn)浮式風(fēng)機(jī)氣動(dòng)—水動(dòng)—?dú)鈴椥择詈享憫?yīng)模擬。

圖4 浮式風(fēng)機(jī)氣動(dòng)—水動(dòng)—?dú)鈴椥择詈享憫?yīng)計(jì)算流程

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

(12)

(13)

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

2 計(jì)算設(shè)置

2.1 浮式風(fēng)機(jī)模型

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

表1 OC3 Hywind Spar型浮式風(fēng)機(jī)主要參數(shù)

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

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

圖5 計(jì)算域布置情況

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

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

2.3 風(fēng)浪耦合入流條件

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

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

3 算例結(jié)果分析

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

在開展浮式風(fēng)機(jī)氣動(dòng)—水動(dòng)—?dú)鈴椥择詈享憫?yīng)數(shù)值模擬之前，首先對(duì)計(jì)算網(wǎng)格進(jìn)行收斂性驗(yàn)證。采用三套不同密度的網(wǎng)格開展收斂性驗(yàn)證計(jì)算，具體網(wǎng)格信息如表2所示。計(jì)算風(fēng)速同樣設(shè)置為11.4 m/s，采用均勻入流條件，不考慮風(fēng)力機(jī)的控制系統(tǒng)。同時(shí)為了減少計(jì)算穩(wěn)定所需的時(shí)間，不放開浮式平臺(tái)的自由度。計(jì)算采用的時(shí)間步長(zhǎng)為0.01 s，小于文獻(xiàn)[21]中對(duì)計(jì)算收斂所需的最小時(shí)間步，同時(shí)也滿足CFL條件。

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

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

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

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

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

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

圖8 風(fēng)力機(jī)葉片(Blade #1)結(jié)構(gòu)變形響應(yīng)

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

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

表4 葉尖位置結(jié)構(gòu)變形響應(yīng)的均方根和標(biāo)準(zhǔn)差

3.3 風(fēng)力機(jī)氣動(dòng)載荷響應(yīng)

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

圖9 浮式風(fēng)機(jī)氣動(dòng)載荷系數(shù)的時(shí)歷曲線

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

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

圖10 浮式風(fēng)機(jī)葉片(Blade #1)軸向入流風(fēng)速(Ua)時(shí)空分布云圖

圖11 浮式風(fēng)機(jī)葉片(Blade #1)攻角時(shí)空分布云圖

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

圖12 風(fēng)力機(jī)葉片(Blade #1)葉根彎矩及偏航力矩的時(shí)歷曲線

3.4 浮式平臺(tái)水動(dòng)力響應(yīng)

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

表5 浮式平臺(tái)受到的氣動(dòng)力載荷與水動(dòng)力載荷統(tǒng)計(jì)表

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

圖14 系泊纜繩張力的時(shí)歷曲線

3.5 浮式風(fēng)機(jī)尾流場(chǎng)特性

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

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

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

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

4 結(jié) 語(yǔ)

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