陳嘉豪, 劉格梁, 胡志強(qiáng)

(1. 上海交通大學(xué) 海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室; 高新船舶與深海開(kāi)發(fā)裝備協(xié)同創(chuàng)新中心，上海 200240; 2. 紐卡斯?fàn)柎髮W(xué) 工程學(xué)院, 紐卡斯?fàn)?NE1 7RU，英國(guó);3. 中國(guó)能源建設(shè)集團(tuán)廣東省電力設(shè)計(jì)研究院有限公司，廣州 510663)

近十年來(lái)，由于各國(guó)日漸重視可再生能源的轉(zhuǎn)型發(fā)展，海上風(fēng)電技術(shù)及其相關(guān)產(chǎn)業(yè)在全球呈現(xiàn)蓬勃發(fā)展的態(tài)勢(shì)[1].而隨著近海優(yōu)質(zhì)風(fēng)電資源的日益開(kāi)發(fā)，海上風(fēng)電產(chǎn)業(yè)也如同過(guò)去油氣產(chǎn)業(yè)所經(jīng)歷的變革一樣，逐漸從陸地走向海洋，從淺海過(guò)渡到中等或更深的水域，其相應(yīng)的支撐平臺(tái)形式也從固定式逐漸過(guò)渡為浮式[2].但由于海上浮式風(fēng)機(jī)具有復(fù)雜的多物理場(chǎng)耦合特性，如何準(zhǔn)確、快速地對(duì)海上浮式風(fēng)機(jī)的耦合動(dòng)力特性進(jìn)行數(shù)值預(yù)報(bào)，成為一直以來(lái)的研究難點(diǎn)與熱點(diǎn)問(wèn)題之一[3].

目前，海上浮式風(fēng)機(jī)的數(shù)值分析主要分為頻域計(jì)算和時(shí)域計(jì)算.早期研究人員主要借鑒離岸油氣工業(yè)的經(jīng)驗(yàn)，利用頻域分析法對(duì)海上浮式風(fēng)機(jī)進(jìn)行設(shè)計(jì)與分析[4-5].雖然，頻域分析法能夠快速簡(jiǎn)便地計(jì)算海上結(jié)構(gòu)物的頻域動(dòng)力響應(yīng)特性，但由于海上浮式風(fēng)機(jī)涉及空氣動(dòng)力學(xué)、水動(dòng)力學(xué)、結(jié)構(gòu)動(dòng)力學(xué)、錨泊動(dòng)力學(xué)及控制系統(tǒng)等技術(shù)，其動(dòng)力響應(yīng)具有高度的非線性特性.頻域分析法無(wú)法全面、準(zhǔn)確地分析海上浮式風(fēng)機(jī)的耦合動(dòng)力響應(yīng)特性.在這之后，國(guó)際上一些研究人員改寫(xiě)或調(diào)用固定式風(fēng)機(jī)或機(jī)械行業(yè)的計(jì)算軟件，用以時(shí)域計(jì)算并分析海上浮式風(fēng)機(jī)的耦合動(dòng)力特性，如FAST[6]、 HAWC2[7]、 Bladed[8]、 Adams、 Simpack等軟件.在國(guó)內(nèi)，唐友剛等[9]利用水動(dòng)力計(jì)算軟件(HydroD)對(duì)自主設(shè)計(jì)的新型海上浮式風(fēng)機(jī)在頻域范圍內(nèi)進(jìn)行動(dòng)力學(xué)仿真.Ma等[10]利用FAST軟件對(duì)立柱型浮式平臺(tái)耦合動(dòng)力響應(yīng)進(jìn)行時(shí)域計(jì)算分析.趙文超等[11]利用計(jì)算流體動(dòng)力學(xué)(CFD)滑移網(wǎng)格技術(shù)對(duì)海上浮式風(fēng)機(jī)的氣動(dòng)性能做了仿真計(jì)算.劉利琴等[12]針對(duì)垂直軸的海上浮式風(fēng)機(jī)開(kāi)發(fā)了氣動(dòng)-水動(dòng)力耦合數(shù)值仿真程序.當(dāng)前，國(guó)內(nèi)外關(guān)于海上浮式風(fēng)機(jī)的全耦合時(shí)域仿真程序尚在研究中，且國(guó)內(nèi)對(duì)于海上浮式風(fēng)機(jī)相關(guān)理論的研究及數(shù)值程序的開(kāi)發(fā)與國(guó)外同行尚有一定差距.

本文介紹自主研發(fā)的用于海上浮式風(fēng)機(jī)的氣動(dòng)-水動(dòng)-結(jié)構(gòu)-錨鏈-控制全耦合時(shí)域數(shù)值仿真程序(DARwind)的相關(guān)理論方法.在氣動(dòng)計(jì)算方面，利用非定常葉素-動(dòng)量理論計(jì)算風(fēng)機(jī)的氣動(dòng)載荷；在水動(dòng)力計(jì)算方面，利用Airy線性波浪理論、勢(shì)流理論以及Morison方程計(jì)算相關(guān)的水動(dòng)力載荷；在錨泊模型方面，采用準(zhǔn)靜態(tài)的懸鏈線錨泊模型計(jì)算錨鏈張力；在多體動(dòng)力學(xué)方面，利用Kane動(dòng)力學(xué)方程構(gòu)建系統(tǒng)的動(dòng)力學(xué)方程，并求解各自由度的加速度信息；在柔性動(dòng)力學(xué)方面，對(duì)于如塔架和槳葉等柔性體，采用高階變形模型以考慮非線性的剛?cè)狁詈闲?yīng)；對(duì)比測(cè)試DARwind程序與FAST軟件，以驗(yàn)證DARwind程序在海上浮式風(fēng)機(jī)數(shù)值仿真方面的有效性.本文可為海上浮式風(fēng)機(jī)耦合動(dòng)力理論及其相關(guān)的數(shù)值仿真程序的研發(fā)提供一定的參考依據(jù).

1 理論模型

1.1 水動(dòng)力模型

DARwind程序中，利用Airy線性波浪理論、勢(shì)流理論以及Morison方程計(jì)算相關(guān)水動(dòng)力載荷[13].部分水動(dòng)力參數(shù)，如：復(fù)原力系數(shù)、波浪激勵(lì)力傳遞函數(shù)、波浪輻射力傳遞函數(shù)、二階波浪力的二次傳遞函數(shù)(QTF)矩陣均由三維頻域勢(shì)流軟件，如WAMIT、Sesam等，進(jìn)行頻域計(jì)算并生成.這些頻域水動(dòng)力參數(shù)將作為DARwind程序的讀入數(shù)據(jù)，再依據(jù)實(shí)際計(jì)算的波幅信息、平臺(tái)的位移、速度和加速度等參數(shù)，實(shí)時(shí)生成時(shí)域水動(dòng)力載荷.其中，黏性阻力采用Morison經(jīng)驗(yàn)公式進(jìn)行修正.

波浪激勵(lì)力FW包含波浪入射力和繞射力，考慮線性化的一階波浪力，并采用諧波疊加法進(jìn)行計(jì)算，即

(1)

式中：上標(biāo)W為波浪激勵(lì)；fn為第n個(gè)單元規(guī)則波的頻率；εn為第n個(gè)單元規(guī)則波的相位；An為第n個(gè)單元規(guī)則波的波幅；HW(ifn)表示頻率為fn的波浪力的頻率響應(yīng)函數(shù)；t為時(shí)間.

考慮自由表面記憶效應(yīng)，波浪輻射力采用間接時(shí)域計(jì)算方法，利用頻域勢(shì)流理論計(jì)算可以獲得輻射水動(dòng)力參數(shù)μjk、λjk，再經(jīng)過(guò)下面的公式，可以得到時(shí)域輻射力

(2)

j=1,2,…,6

(3)

沿濕表面對(duì)平臺(tái)在水中受到的靜壓力積分可以獲得浮體的靜水力載荷，分為平臺(tái)正浮時(shí)的浮力和平臺(tái)發(fā)生位移時(shí)的靜水回復(fù)力，其中靜水回復(fù)力的表達(dá)式為

(4)

對(duì)于平臺(tái)的細(xì)長(zhǎng)構(gòu)件，由流體黏性引起的載荷較為明顯，勢(shì)流理論無(wú)法準(zhǔn)確地計(jì)算該載荷.此時(shí)，可根據(jù)Morison經(jīng)驗(yàn)公式對(duì)勢(shì)流載荷進(jìn)行黏性修正.圓柱體在lx處離散的切片長(zhǎng)度為dl的流體黏性力FV表達(dá)式為

(5)

式中：上標(biāo)V為黏性；CD為黏性系數(shù)，一般由實(shí)驗(yàn)測(cè)得；D為圓柱直徑；dl為圓柱離散的切片厚度；νw為水質(zhì)點(diǎn)在切片處垂直于切片軸向的速度分量，其值等于波浪運(yùn)動(dòng)導(dǎo)致的水質(zhì)點(diǎn)的速度與海流引起的速度在該方向的速度分量之和；νs為切片的運(yùn)動(dòng)速度在垂直于切片軸向的速度分量.黏性力矩陣包含集中力和力矩，沿著圓柱總長(zhǎng)度LC進(jìn)行積分，表達(dá)式為

(6)

因此，總的水動(dòng)力FH根據(jù)線性理論，為上述FW、FR、FS和FV之和，則有

FH=FW+FR+FS+FV

(7)

1.2 錨泊動(dòng)力模型

錨泊系統(tǒng)對(duì)海上浮式風(fēng)機(jī)平臺(tái)起到基本的定位和約束作用，同時(shí)影響著浮式結(jié)構(gòu)物的運(yùn)動(dòng)響應(yīng)特性.DARwind程序采用準(zhǔn)靜態(tài)的懸鏈線錨泊模型[14]計(jì)算錨泊對(duì)平臺(tái)的約束載荷.其中，除錨泊的伸長(zhǎng)變形外的一些動(dòng)態(tài)特性，如慣性力、阻尼力、錨鏈變形彎矩均被忽略不計(jì).

圖1 懸鏈線錨泊示意圖Fig.1 The schematic diagram of catenary mooring lines

懸鏈線錨泊模型如圖1所示.其中，錨點(diǎn)受到水平方向張力HA和垂直方向張力VA的作用.假定躺底段(錨鏈與海床接觸的部分)長(zhǎng)為lb；錨鏈任意一點(diǎn)的張力為T(mén)s′，水平方向的張力為Hs′，垂直方向的張力為Vs′，其距離錨點(diǎn)的水平距離為xs′，垂直距離為zs′，從躺底點(diǎn)算起的錨鏈自然長(zhǎng)度為ls′；導(dǎo)纜孔的張力為T(mén)F，該張力水平方向的分力為HF，垂直方向的分力為VF；導(dǎo)纜孔到海底的垂直距離為zF，錨鏈整體的水平投影長(zhǎng)度(水平距離)為xF；錨鏈的總長(zhǎng)度為L(zhǎng)，濕密度為w, 截面面積為A，軸向剛度為EA.

當(dāng)錨鏈與海底接觸時(shí)，即lb≠0，錨點(diǎn)處此時(shí)沒(méi)有垂向張力作用于錨鏈.根據(jù)導(dǎo)纜孔和錨點(diǎn)之間的水平距離和垂向距離，平臺(tái)導(dǎo)纜孔處對(duì)錨鏈的水平張力和垂向張力有如下關(guān)系式：

(8)

(9)

當(dāng)lb=0時(shí)，錨點(diǎn)對(duì)錨鏈的垂向張力等于導(dǎo)纜孔對(duì)錨鏈的垂向張力減去錨泊質(zhì)量，關(guān)系式為

VA=VF-wL

(10)

此時(shí)，導(dǎo)纜孔處的水平張力和垂向張力的關(guān)系式為

(11)

(12)

只要獲得錨鏈當(dāng)前時(shí)刻的錨點(diǎn)和導(dǎo)纜孔的位置(一般情況下，錨點(diǎn)的位置保持不變，導(dǎo)纜孔的位置需要根據(jù)當(dāng)前的平臺(tái)位置信息獲得)，通過(guò)上述非線性方程組即可求得導(dǎo)纜孔對(duì)平臺(tái)的約束張力.其中，非線性方程的數(shù)值求解可以借助迭代技術(shù)，如牛頓迭代法等進(jìn)行數(shù)值求解.

1.3 空氣動(dòng)力模型

葉素-動(dòng)量法[15]結(jié)合了一維的動(dòng)量理論和二維的葉素理論.動(dòng)量理論是基于理想的無(wú)窮槳葉正對(duì)來(lái)流時(shí)的致動(dòng)盤(pán)假定推導(dǎo)而獲得的.葉素理論則考慮了槳葉的實(shí)際形狀，將槳葉劃分為許多二維的翼型剖面，稱(chēng)之為“葉素”.根據(jù)每一個(gè)葉素上實(shí)際入流風(fēng)速計(jì)算獲得局部入流攻角，然后根據(jù)入流攻角從翼型氣動(dòng)系數(shù)的數(shù)據(jù)庫(kù)中插值獲得當(dāng)前翼型剖面的氣動(dòng)系數(shù)，如升力系數(shù)、阻力系數(shù)和轉(zhuǎn)矩系數(shù)等.根據(jù)這些系數(shù)，結(jié)合相對(duì)入流風(fēng)速即可計(jì)算獲得該剖面的氣動(dòng)載荷.最后，再通過(guò)積分(疊加)各個(gè)槳葉剖面葉素的氣動(dòng)載荷，計(jì)算總的風(fēng)機(jī)氣動(dòng)載荷.令由動(dòng)量理論得到的氣動(dòng)載荷和由葉素理論得到的氣動(dòng)載荷在環(huán)形控制體(見(jiàn)圖2(a))內(nèi)保持相等，即可得到各不同徑長(zhǎng)的環(huán)形控制體內(nèi)的槳葉葉素的誘導(dǎo)系數(shù)，其表達(dá)式為

(13)

(14)

式中：a為軸向誘導(dǎo)系數(shù)；a′為切向誘導(dǎo)系數(shù)；σ=crbBblade/(2πrb)為風(fēng)輪實(shí)度，其中Bblade為槳葉數(shù)，rb為葉素的徑向長(zhǎng)度，crb為葉素局部弦長(zhǎng)；φ為入流角；CL和CD分別是葉素升力和阻力系數(shù).

圖2 氣動(dòng)模型示意圖Fig.2 The schematic diagram of aerodynamic model

在經(jīng)典的葉素動(dòng)量理論中，假定葉輪只有轉(zhuǎn)動(dòng)，且來(lái)流是空間均勻和時(shí)間定常的.但對(duì)于實(shí)際的海上浮式風(fēng)機(jī)而言，空氣流動(dòng)是時(shí)空變化的，浮式平臺(tái)的運(yùn)動(dòng)也會(huì)引起槳葉的運(yùn)動(dòng)，加劇其相對(duì)入流風(fēng)速的變化，氣動(dòng)誘導(dǎo)系數(shù)也隨時(shí)空所變化，因此必須考慮非定常的局部氣動(dòng)模型，如圖2所示.其中，α為葉素的局部入流攻角.為了簡(jiǎn)化模型，各葉素的氣動(dòng)性能只考慮小范圍控制體內(nèi)的影響.其研究的單元?dú)鈩?dòng)影響區(qū)域?yàn)橐粋€(gè)扇形區(qū)域，如圖2(c)所示.該扇形區(qū)域面積是圖2(a)所示的環(huán)形控制體區(qū)域面積的1/3，即在當(dāng)前時(shí)刻每個(gè)槳葉的葉素氣動(dòng)計(jì)算只考慮該槳葉附近區(qū)域的影響.

當(dāng)不考慮槳葉柔性變形和平臺(tái)運(yùn)動(dòng)所帶來(lái)的葉素運(yùn)動(dòng)，葉素的相對(duì)入流風(fēng)速如圖2(b) 所示，其表達(dá)式為

(15)

式中：v0為上游未受干擾的入流風(fēng)速大??；Ω為風(fēng)輪的轉(zhuǎn)度；vrel為最終合成的相對(duì)入流風(fēng)速矢量；在非定常的局部葉素-動(dòng)量模型中，DARwind程序考慮了局部槳葉的葉素因平臺(tái)運(yùn)動(dòng)、塔架變形等引起的運(yùn)動(dòng)速度，以及由槳葉本身變形導(dǎo)致的葉素運(yùn)動(dòng)速度.葉素-動(dòng)量法雖然計(jì)算速度快，且過(guò)去一直作為風(fēng)能領(lǐng)域常用的可靠計(jì)算方法，但由于其并不能完全描述風(fēng)場(chǎng)及風(fēng)機(jī)的細(xì)節(jié)信息，且計(jì)算精度在特定場(chǎng)合可能有一定的誤差，故需要引入一些修正方法[16]，例如：通常情況下，風(fēng)機(jī)槳盤(pán)面并不正對(duì)來(lái)流風(fēng)，此時(shí)稱(chēng)之為“偏航”狀態(tài).對(duì)于海上浮式風(fēng)機(jī)而言，由于浮式平臺(tái)的大幅度運(yùn)動(dòng)，這種偏航狀態(tài)更為常見(jiàn)，故需要考慮偏航狀態(tài)下的氣動(dòng)修正；另外，由于動(dòng)量理論假定葉片數(shù)量為無(wú)窮多，但實(shí)際的風(fēng)輪葉片數(shù)量是有限的，風(fēng)機(jī)尾流渦系與無(wú)窮槳葉的理想風(fēng)輪的尾流場(chǎng)有差別，此時(shí)需要采用Prandtl葉尖損失模型進(jìn)行修正；當(dāng)軸向誘導(dǎo)因子約大于0.4時(shí)，簡(jiǎn)單的動(dòng)量理論不再適用，此時(shí)需要使用 Glarert 經(jīng)驗(yàn)公式[16-17]進(jìn)行適當(dāng)?shù)男拚?，以使?shù)值計(jì)算結(jié)果與實(shí)際結(jié)果更為符合.上述所提及的修正理論均被采用并編入DARwind程序的氣動(dòng)模型中，程序中的葉素-動(dòng)量法是一個(gè)反復(fù)迭代求解誘導(dǎo)系數(shù)直至收斂的過(guò)程，其計(jì)算速度關(guān)系到程序整體的仿真速度，所以此模塊還需要采用一些算法的優(yōu)化技術(shù)以提高迭代計(jì)算的收斂速度.

1.4 動(dòng)力學(xué)模型

海上浮式風(fēng)機(jī)的多個(gè)不同部件擁有不同的特性，可各自根據(jù)相對(duì)變形量簡(jiǎn)化為剛性體或柔性體，因此海上浮式風(fēng)機(jī)模型為一個(gè)剛?cè)狁詈系亩囿w系統(tǒng).采用浮動(dòng)坐標(biāo)系[18]描述剛?cè)狁詈隙囿w系統(tǒng)的運(yùn)動(dòng)關(guān)系.該方法采用兩套坐標(biāo)系描述每個(gè)部件的位置、方向及變形：① 全局坐標(biāo)系，用以描述各部件的位置和方向；② 各部件的隨體坐標(biāo)系，用以描述部件的柔性變形場(chǎng).其中對(duì)于柔性體，如塔架和槳葉，其變形場(chǎng)采用模態(tài)疊加法進(jìn)行有限模態(tài)的近似離散，變形量U如下所示：

U=ΘQ

(16)

式中：Θ為模態(tài)形狀函數(shù)矩陣；Q為模態(tài)坐標(biāo)(廣義坐標(biāo))矩陣.研究表明[19-20]：當(dāng)柔性體進(jìn)行大范圍運(yùn)動(dòng)時(shí)，尤其是高速旋轉(zhuǎn)運(yùn)動(dòng)，用傳統(tǒng)的線性化方法將會(huì)產(chǎn)生較大的數(shù)值誤差.因此在DARwind程序中，柔性體(塔架和槳葉)均被近似為三維的Euler-Bernoulli 模型.當(dāng)發(fā)生大范圍的運(yùn)動(dòng)時(shí)，梁的側(cè)向變形會(huì)導(dǎo)致軸向的耦合收縮效應(yīng)，并且隨著運(yùn)動(dòng)頻率的增加，柔性體的剛度將會(huì)隨之增加，該效應(yīng)稱(chēng)為剛?cè)峄旌隙囿w的“動(dòng)力剛化”.三維梁模型的高階非線性耦合變形量可簡(jiǎn)化為

(17)

式中：ush是由梁結(jié)構(gòu)側(cè)向變形引起的軸向收縮量；uy0和uz0分別是梁側(cè)向變形量關(guān)于局部隨體坐標(biāo)系y，z方向的分量.當(dāng)考慮高階非線性耦合變形效應(yīng)時(shí)，上述的柔性變形可以表示為

(18)

其中非線性高階耦合效應(yīng)矩陣H為

(19)

(20)

式中：下標(biāo)x為對(duì)應(yīng)物理量的坐標(biāo)分量.通過(guò)上述的浮動(dòng)坐標(biāo)系法及模態(tài)離散，可以描述系統(tǒng)各部件的各個(gè)質(zhì)點(diǎn)之間的運(yùn)動(dòng)學(xué)模型.研究中需要構(gòu)建系統(tǒng)的動(dòng)力學(xué)控制方程，以求解系統(tǒng)在外力和系統(tǒng)約束作用下各自由度的信息.采用Kane動(dòng)力學(xué)方程[21]構(gòu)建系統(tǒng)動(dòng)力學(xué)控制方程.這種方法兼具矢量力學(xué)和分析力學(xué)特點(diǎn)，既避免了出現(xiàn)理想約束反力，也避免了對(duì)動(dòng)力學(xué)函數(shù)進(jìn)行求導(dǎo)，減少了計(jì)算工作量，且易于拓展和推廣.

一個(gè)質(zhì)點(diǎn)系統(tǒng)，關(guān)于第j個(gè)廣義速度的Kane方程可表示為

Fj+F*j=0

(21)

廣義主動(dòng)力Fj以及廣義慣性力F*j分別計(jì)算如下：

(24)

(25)

式中：Bki和Dki分別為部件i的平動(dòng)和轉(zhuǎn)動(dòng)的加速度余項(xiàng)，其解析表達(dá)式較為復(fù)雜，可通過(guò)線速度和角速度關(guān)于時(shí)間的一階導(dǎo)數(shù)推導(dǎo)而得，需要注意各個(gè)部件的坐標(biāo)系之間的轉(zhuǎn)換關(guān)系.

對(duì)于剛?cè)崴欧囿w耦合的海上浮式風(fēng)機(jī)系統(tǒng)，其Kane方程可表示為

(28)

Nr+Nf+Ns=N

將式(28)中的各部件的廣義主動(dòng)力和慣性力分別進(jìn)行累計(jì)計(jì)算，則有

(29)

(30)

式中：Fi和Mi分別為部件i的動(dòng)力學(xué)參考點(diǎn)的力和力矩；mi和Ii分別為部件i的質(zhì)量和轉(zhuǎn)動(dòng)慣量；式(29)右端第一項(xiàng)為剛體和柔性體的廣義主動(dòng)力，比起剛體，柔性體中存在額外的結(jié)構(gòu)變形能Ei；式(29)右端的第二項(xiàng)為控制器的廣義主動(dòng)力或內(nèi)力；式(30)右端為剛體和柔性體的廣義慣性力和力矩.在海上浮式風(fēng)機(jī)系統(tǒng)中，浮式支撐平臺(tái)、機(jī)艙、轉(zhuǎn)軸和輪轂結(jié)構(gòu)均視為剛體，塔架和槳葉均視為柔性體.需將柔性體劃分為許多微小單元，對(duì)每個(gè)微單元分別建立Kane方程，故柔性體的存在增加了動(dòng)力學(xué)方程的復(fù)雜度和計(jì)算量.

[IiDi+ωi×(Iiωi)]}=0

(31)

令

(32)

[IiDi+ωi×(Iiωi)]}

(33)

(34)

1.5 控制系統(tǒng)模型

DARwind程序采用發(fā)電機(jī)轉(zhuǎn)矩控制器和槳距角控制器.當(dāng)風(fēng)速小于額定風(fēng)速時(shí)，以?xún)?yōu)化風(fēng)能捕獲效率為目標(biāo)，通過(guò)調(diào)節(jié)發(fā)電機(jī)轉(zhuǎn)矩進(jìn)而調(diào)節(jié)風(fēng)機(jī)的最佳葉尖速比；當(dāng)風(fēng)速達(dá)到或超過(guò)額定風(fēng)速后，則主要通過(guò)改變槳距角控制系統(tǒng)調(diào)節(jié)氣動(dòng)性能，進(jìn)而達(dá)到穩(wěn)定功率的目的.

圖3 發(fā)電機(jī)轉(zhuǎn)矩控制Fig.3 Generator-torque control

槳距角控制器采用比例-積分統(tǒng)一變槳機(jī)制，其表達(dá)式為

(35)

式中：θ為槳距角；Δθ為當(dāng)前需要改變的槳距角度；t為槳距角控制器的執(zhí)行時(shí)間；ΔΩ為當(dāng)前風(fēng)輪轉(zhuǎn)速與目標(biāo)轉(zhuǎn)速之差；KP=KPc/[1+(θ/θc)]為槳距角控制器的比例項(xiàng)，其中θc為當(dāng)槳距角為此值時(shí)，功率對(duì)槳距角的導(dǎo)數(shù)為槳距角等于0時(shí)所對(duì)應(yīng)的導(dǎo)數(shù)的2倍，KPc為θc槳距角下的比例項(xiàng)；KI=KIc/ [1+(θ/θc)]為槳距角控制器的積分項(xiàng)，KIc為θc槳距角下的積分項(xiàng).二者的具體取值需要根據(jù)實(shí)時(shí)的槳距角動(dòng)態(tài)進(jìn)行調(diào)整，其推導(dǎo)過(guò)程可參考文獻(xiàn)[22].

1.6 程序計(jì)算流程

運(yùn)用上述理論模型構(gòu)建用于海上浮式風(fēng)機(jī)動(dòng)力性能預(yù)報(bào)的DARwind程序，計(jì)算流程如下：

(1) 程序開(kāi)始時(shí)，讀入標(biāo)記參數(shù)、結(jié)構(gòu)信息及環(huán)境信息等，用以初始化計(jì)算模型.

(2) 在每一個(gè)計(jì)算時(shí)間步下，程序需要根據(jù)平臺(tái)的位置、速度信息及平臺(tái)參考點(diǎn)處的波浪場(chǎng)信息生成相應(yīng)的水動(dòng)力載荷；根據(jù)槳葉的變形、運(yùn)動(dòng)及局部槳葉在當(dāng)前時(shí)空的風(fēng)場(chǎng)風(fēng)速，利用非穩(wěn)態(tài)葉素-動(dòng)量法計(jì)算氣動(dòng)載荷；根據(jù)平臺(tái)的位置信息計(jì)算導(dǎo)纜孔的位置；利用懸鏈線錨泊模型計(jì)算當(dāng)前錨鏈對(duì)平臺(tái)的約束力.另外，程序還考慮了重力、浮力及控制器的作用力等參數(shù).

(3) 經(jīng)過(guò)步驟(2)求得所有環(huán)境載荷后，通過(guò)Kane方程構(gòu)建系統(tǒng)動(dòng)力學(xué)控制方程，求得當(dāng)前各自由度的加速度信息，并利用4階Runge-Kutta時(shí)間積分法遞推下一時(shí)刻各自由度的位置及速度信息.

(4) 程序依次循環(huán)直到預(yù)設(shè)的結(jié)束時(shí)間.

相比于其他類(lèi)似程序，DARwind程序有兩方面的創(chuàng)新[23]：一方面，該程序可以非常簡(jiǎn)便地根據(jù)不同計(jì)算目的，對(duì)海上浮式風(fēng)機(jī)系統(tǒng)建立不同的多體系統(tǒng)方案.例如：當(dāng)需要快速獲得海上浮式風(fēng)機(jī)系統(tǒng)的動(dòng)力特性時(shí)，可將整個(gè)系統(tǒng)作為單體模型；當(dāng)需要適中的計(jì)算速度和精度時(shí)，可以將整個(gè)系統(tǒng)作為多剛體結(jié)構(gòu)，忽略柔性結(jié)構(gòu)的影響；當(dāng)需要較為準(zhǔn)確的系統(tǒng)模型時(shí)，可將海上浮式風(fēng)機(jī)系統(tǒng)作為剛?cè)狁詈隙囿w模型.另一方面，DARwind程序中創(chuàng)新地采用了高階的剛?cè)狁詈夏Ｐ?，考慮了剛體運(yùn)動(dòng)對(duì)柔性體動(dòng)力性能的影響.

2 對(duì)比驗(yàn)證

DARwind程序?qū)⑴c美國(guó)國(guó)家可再生能源實(shí)驗(yàn)室(NREL)的風(fēng)機(jī)仿真FAST 程序[6]進(jìn)行對(duì)比分析，以驗(yàn)證DARwind程序在海上浮式風(fēng)機(jī)數(shù)值仿真方面的有效性.采用OC3-Hywind Spar海上浮式風(fēng)機(jī)作為研究對(duì)象，如圖4所示，關(guān)于該海上浮式風(fēng)機(jī)更詳細(xì)的信息可以參考文獻(xiàn)[24].

圖4 OC3-Hywind Spar海上浮式風(fēng)機(jī)[24]Fig.4 OC3-Hywind Spar floating offshore wind turbine[24]

2.1 固有頻率測(cè)試

在相同的模型參數(shù)設(shè)置下，DARwind和FAST程序在自由衰減運(yùn)動(dòng)中測(cè)得的固有頻率以及無(wú)因次阻尼系數(shù)如表1所示.由表1可知, 除了縱搖運(yùn)動(dòng)具有微小的差異之外，DARwind以及FAST程序的自由衰減運(yùn)動(dòng)的特性都是非常接近的.該結(jié)果表明DARwind程序能準(zhǔn)確地反映海上浮式風(fēng)機(jī)系統(tǒng)的質(zhì)量、慣性、靜水回復(fù)力、錨鏈剛度及水動(dòng)力阻尼等特性.

表1 DARwind和FAST程序中的運(yùn)動(dòng)固有頻率與無(wú)因次阻尼系數(shù)

Tab.1 Natural frequency and damping coefficient of platform motions in DARwind and FAST softwares

運(yùn)動(dòng)模態(tài)固有頻率/HzDARwindFAST無(wú)因次阻尼系數(shù)DARwindFAST縱蕩0.00810.00800.48770.5286垂蕩0.03230.03230.24230.2415縱搖0.03140.03340.27590.3300首搖0.12180.12090.27230.2820

2.2 風(fēng)工況

海上浮式風(fēng)機(jī)是用以捕獲風(fēng)能的離岸設(shè)備，對(duì)其氣動(dòng)性能及結(jié)構(gòu)變形的預(yù)測(cè)對(duì)海上浮式風(fēng)機(jī)而言同樣也是需要重點(diǎn)考慮的因素.對(duì)DARwind程序進(jìn)行穩(wěn)態(tài)風(fēng)工況測(cè)試，測(cè)試結(jié)果如表2所示.數(shù)值計(jì)算中，風(fēng)速空間均勻且時(shí)間定常；轉(zhuǎn)速為風(fēng)輪轉(zhuǎn)速；槳距角設(shè)置為0°.當(dāng)不考慮控制系統(tǒng)的作用且平臺(tái)固定時(shí)，該工況的測(cè)試結(jié)果如表3所示.由表3可知，DARwind 程序的氣動(dòng)載荷值(氣動(dòng)推力和氣動(dòng)轉(zhuǎn)矩)比FAST程序的氣動(dòng)載荷值略大，這可能是由于兩個(gè)數(shù)值軟件在氣動(dòng)理論方法上的差異所導(dǎo)致的.DARwind程序采用的是基于葉素扇形影響區(qū)域的局部非穩(wěn)態(tài)氣動(dòng)模型；而FAST程序則采用的是基于風(fēng)輪整體的穩(wěn)態(tài)氣動(dòng)模型.在塔架和槳葉的變形方面，相比于FAST程序，DARwind程序的計(jì)算結(jié)果相對(duì)較小，但兩者差別不大.這主要是由于DARwind程序采用了高階的剛?cè)狁詈夏Ｐ?見(jiàn)式(17)～(20))，考慮了剛體運(yùn)動(dòng)對(duì)柔性體的影響.當(dāng)柔性體發(fā)生大范圍運(yùn)動(dòng)(如高速轉(zhuǎn)動(dòng))時(shí)，柔性體的剛度會(huì)隨之提高，變形也會(huì)隨之減小.總體而言，無(wú)論是氣動(dòng)性能還是結(jié)構(gòu)彈性動(dòng)力特性，兩個(gè)程序的計(jì)算結(jié)果都比較接近，這驗(yàn)證了上述柔性體變形理論及葉素-動(dòng)量理論在海上漂浮式風(fēng)機(jī)時(shí)域計(jì)算的有效性.

表2 穩(wěn)態(tài)風(fēng)工況測(cè)試Tab.2 Steady wind cases

表3 氣動(dòng)載荷與葉尖變形量Tab.3 Aerodynamic loads and blade-tip deflection

2.3 風(fēng)波聯(lián)合工況

海上浮式風(fēng)機(jī)工作于有風(fēng)及波浪載荷共同作用的海洋環(huán)境下，因此，對(duì)浮式平臺(tái)運(yùn)動(dòng)特性的預(yù)報(bào)顯得尤為重要.下文將驗(yàn)證DARwind程序在風(fēng)波聯(lián)合工況下的運(yùn)動(dòng)特性.測(cè)試風(fēng)速為額定風(fēng)速(11.4 m/s)；風(fēng)輪保持額定轉(zhuǎn)速(12.1 r/min)；槳距角統(tǒng)一設(shè)置為0°；不考慮控制器的作用；波浪模型參考“北海聯(lián)合海浪計(jì)劃”(Jonswap)波浪譜，有義波高Hs為2.0 m；譜峰周期Tp為8.0 s；譜形參數(shù)r設(shè)置為3.3.DARwind 和 FAST 程序的縱蕩及縱搖的功率譜密度對(duì)比如圖5所示.由圖5可知，兩種數(shù)值程序在風(fēng)波聯(lián)合工況下的運(yùn)動(dòng)特性吻合得較好，且均能反映出縱蕩及縱搖運(yùn)動(dòng)中明顯的共振響應(yīng)、波浪能量區(qū)的強(qiáng)迫運(yùn)動(dòng)響應(yīng)、縱蕩及縱搖的耦合運(yùn)動(dòng)效應(yīng).

圖5 平臺(tái)運(yùn)動(dòng)功率譜密度Fig.5 Power spectral density of platform motion

2.4 控制系統(tǒng)模型

在海上浮式風(fēng)機(jī)正常運(yùn)轉(zhuǎn)的過(guò)程中，控制系統(tǒng)起著重要的作用，不但能夠調(diào)節(jié)海上浮式風(fēng)機(jī)的氣動(dòng)載荷與發(fā)電功率，而且對(duì)海上浮式風(fēng)機(jī)的運(yùn)動(dòng)響應(yīng)有著重要的影響.在變槳距控制過(guò)程中，通過(guò)槳葉的槳距角、風(fēng)機(jī)功率等參數(shù)的變化情況，驗(yàn)證 DARwind 程序控制系統(tǒng)的有效性.測(cè)試中，海上浮式風(fēng)機(jī)的平臺(tái)受載荷運(yùn)動(dòng)的影響，對(duì)控制系統(tǒng)的要求比較高，有別于傳統(tǒng)的固定式風(fēng)機(jī).在該控制系統(tǒng)下，槳距角隨時(shí)間的變化情況如圖6(a)所示，發(fā)電機(jī)功率隨時(shí)間的變化情況如圖6(b)所示.由圖6(a)可知，上述測(cè)試過(guò)程中的槳距角隨時(shí)間的變化呈現(xiàn)出波動(dòng)，這是由于平臺(tái)運(yùn)動(dòng)導(dǎo)致槳葉的相對(duì)入流風(fēng)速急劇變化導(dǎo)致的.此時(shí)，為了保持發(fā)電機(jī)功率的穩(wěn)定，槳距角也需要進(jìn)行相應(yīng)的調(diào)整.對(duì)比這兩個(gè)程序的計(jì)算結(jié)果，總體上DARwind程序的槳距角變化波動(dòng)稍微大一些，但是兩者的變化情況基本相同.圖6(b)中的發(fā)電機(jī)功率也有類(lèi)似的情況.其主要原因是由于在額定工況時(shí)，F(xiàn)AST程序的轉(zhuǎn)矩控制依然在進(jìn)行調(diào)整，而DARwind程序在額定工況下選擇穩(wěn)定發(fā)電機(jī)轉(zhuǎn)矩，以減少軸系疲勞損傷.此時(shí)，變槳距控制器成為了DARwind程序唯一的發(fā)電機(jī)功率調(diào)節(jié)手段，所以其槳距角的波動(dòng)會(huì)相對(duì)大一些.

圖6 控制系統(tǒng)下的槳距角及發(fā)電機(jī)功率隨時(shí)間的變化情況Fig.6 Blade pitch angle and generator power vary with time

3 結(jié)語(yǔ)

本文介紹了海上浮式風(fēng)機(jī)時(shí)域耦合仿真程序——DARwind程序的相關(guān)理論，并做了相應(yīng)的數(shù)值驗(yàn)證分析.研究結(jié)果表明：所采用的理論方法在海上浮式風(fēng)機(jī)時(shí)域耦合計(jì)算中是有效的，DARwind程序能夠有效地仿真并準(zhǔn)確地預(yù)報(bào)海上浮式風(fēng)機(jī)復(fù)雜的時(shí)域耦合動(dòng)力響應(yīng)特性.然而，當(dāng)前的DARwind程序仍有進(jìn)一步發(fā)展和優(yōu)化的空間，如：可以進(jìn)一步提升氣動(dòng)模型精度，尤其對(duì)于海上浮式風(fēng)機(jī)大范圍搖擺運(yùn)動(dòng)導(dǎo)致的復(fù)雜尾流場(chǎng)而言；當(dāng)前的控制系統(tǒng)也可以采用更為高級(jí)的控制方式，如神經(jīng)網(wǎng)絡(luò)控制、模糊控制等現(xiàn)代智能化控制方式.