湯群益，孫震洲，陳杰峰，金 磊

(1. 浙江省深遠(yuǎn)海風(fēng)電技術(shù)研究重點(diǎn)實(shí)驗(yàn)室，浙江 杭州 311122； 2. 中國(guó)電建集團(tuán)華東勘測(cè)設(shè)計(jì)研究院有限公司，浙江 杭州 311122； 3. 三峽新能源山東昌邑發(fā)電有限公司，山東 濰坊 261300)

隨著海上風(fēng)電建設(shè)的推進(jìn)，海上風(fēng)電逐漸由近海走向深海，為了適應(yīng)海上風(fēng)電的發(fā)展趨勢(shì)，海上風(fēng)電場(chǎng)中的海工結(jié)構(gòu)物結(jié)構(gòu)型式也由固定式轉(zhuǎn)變?yōu)楦∈?。作為電力輸送、匯集的樞紐，浮式海上升壓站的安全性對(duì)海上風(fēng)電場(chǎng)的健康運(yùn)營(yíng)至關(guān)重要。相比于固定式升壓站，浮式海上升壓站在6自由度上的響應(yīng)更為明顯，這也給升壓站中的電氣設(shè)備帶來(lái)較大的安全隱患[1]。因此，需要在設(shè)計(jì)階段對(duì)浮式海上升壓站結(jié)構(gòu)進(jìn)行動(dòng)力響應(yīng)分析，進(jìn)而優(yōu)化其水動(dòng)力性能，保障浮式海上升壓站的平穩(wěn)運(yùn)行。

基于勢(shì)流理論，同船體等浮式結(jié)構(gòu)相同，浮式海上升壓站的運(yùn)動(dòng)常采用Cummins方程進(jìn)行描述[2]。為了求解Cummins方程得到浮式結(jié)構(gòu)的動(dòng)力響應(yīng)，常用的方法包括時(shí)域法和頻域法。頻域法通過(guò)傅里葉變換，在頻域內(nèi)求解Cummins方程，避免卷積項(xiàng)積分。但在運(yùn)算時(shí)應(yīng)用傅里葉變換，會(huì)導(dǎo)致傅里葉變換漏頻、能量泄露等局限性突顯[3]，并且無(wú)法計(jì)算浮式結(jié)構(gòu)瞬態(tài)響應(yīng)。時(shí)域法又可分為直接時(shí)域法和間接時(shí)域法。直接時(shí)域法直接在時(shí)域內(nèi)構(gòu)建求解速度勢(shì)的初邊值偏微分方程，通過(guò)求解時(shí)域格林函數(shù)來(lái)分析浮式結(jié)構(gòu)物的動(dòng)力響應(yīng)。直接時(shí)域提出時(shí)間較早，但在求解浮式結(jié)構(gòu)的水動(dòng)力過(guò)程中，面臨計(jì)算量大、分析效率低等問(wèn)題。相比直接時(shí)域法，間接時(shí)域法采用頻域勢(shì)流理論計(jì)算浮式平臺(tái)的水動(dòng)力參數(shù)及波浪力，通過(guò)傅里葉變換將上述頻域水動(dòng)力參數(shù)和波浪力轉(zhuǎn)化到時(shí)域，從而在時(shí)域內(nèi)求解浮式結(jié)構(gòu)的動(dòng)力響應(yīng)。間接時(shí)域法能夠分析浮式結(jié)構(gòu)的瞬時(shí)動(dòng)力響應(yīng)，并且比直接時(shí)域法的計(jì)算效率高。然而，間接時(shí)域法求解浮式結(jié)構(gòu)的動(dòng)力響應(yīng)時(shí)，依賴于卷積運(yùn)算，從而消耗大量計(jì)算資源[4]。同時(shí)，間接時(shí)域法由于高頻位置的水動(dòng)力系數(shù)難以精確求得，對(duì)最終的動(dòng)力響應(yīng)求解造成難以控制的誤差[5]。針對(duì)這些問(wèn)題，一種思路是通過(guò)狀態(tài)空間模型來(lái)代替線性卷積項(xiàng)，從而提高浮式結(jié)構(gòu)動(dòng)力響應(yīng)分析的計(jì)算速度和精度。對(duì)于此，學(xué)者們展開(kāi)了大量研究。Schmiechen[6]將狀態(tài)空間模型和船體的瞬態(tài)響應(yīng)建立起聯(lián)系。Xia等[7]也在相關(guān)海洋結(jié)構(gòu)物的水動(dòng)力分析中使用了狀態(tài)空間模型并進(jìn)行了深入研究。Sutulo和Guedes[8]進(jìn)一步提出可以用狀態(tài)空間形式來(lái)表達(dá)輻射力，進(jìn)而代替?zhèn)鹘y(tǒng)時(shí)域方程中的卷積項(xiàng)。Taghipour等[9]應(yīng)用時(shí)域方程直接計(jì)算卷積方法和狀態(tài)空間方法做了相應(yīng)的算例分析，并給出了方法的詳細(xì)介紹。然而，這些方法大多通過(guò)對(duì)興波阻尼進(jìn)行余弦變換得到延遲函數(shù)，進(jìn)而求解對(duì)應(yīng)的狀態(tài)空間模型。此過(guò)程采用梯形積分法計(jì)算延遲函數(shù)，額外引入計(jì)算誤差。

為了解決卷積項(xiàng)導(dǎo)致浮式海上升壓站動(dòng)力響應(yīng)分析效率低及采用梯形積分法引起計(jì)算誤差的問(wèn)題，提出了一種基于狀態(tài)空間模型的浮式海上升壓站動(dòng)力響應(yīng)分析方法，同時(shí)運(yùn)用附加質(zhì)量和興波阻尼在頻域內(nèi)計(jì)算延遲函數(shù)對(duì)應(yīng)的狀態(tài)空間，從而提高浮式海上升壓站動(dòng)力響應(yīng)的分析效率和精度。通過(guò)將文中方法應(yīng)用于日本福島浮式海上風(fēng)電示范項(xiàng)目中的浮式升壓站模型[10]，并與商業(yè)軟件SESAM[11]計(jì)算結(jié)果進(jìn)行對(duì)比，檢驗(yàn)文中方法的正確性和有效性。

1 基礎(chǔ)理論

1.1 浮式平臺(tái)運(yùn)動(dòng)方程

在勢(shì)流理論的假設(shè)條件下，只考慮一階水動(dòng)力作用，無(wú)航速海上浮式結(jié)構(gòu)的運(yùn)動(dòng)方程可以用Cummins方程表示：

(1)

1.2 Oglivie關(guān)系式

Cuminns方程中附加質(zhì)量A和延遲函數(shù)K(t)是與海上浮式結(jié)構(gòu)水動(dòng)力參數(shù)相關(guān)的量，為探討它們之間的數(shù)學(xué)關(guān)系，對(duì)式(1)進(jìn)行傅里葉變換[12]：

(2)

Ogilvie通過(guò)對(duì)船體運(yùn)動(dòng)預(yù)測(cè)的研究，建立了隨頻率變化的附加質(zhì)量A(ω)及興波阻尼B(ω)與延遲函數(shù)K(t)的關(guān)系式，即[10]：

(3)

(4)

延遲函數(shù)K(t)的傅里葉變換通過(guò)A(ω)和B(ω)可表示為：

(5)

由式(3)可知，式(1)中附加質(zhì)量為：

(6)

從而，式(1)可表示為：

(7)

2 狀態(tài)空間模型頻域算法

2.1 頻域線性回歸擬合

(8)

對(duì)式(8)進(jìn)行Laplace變換，即可得到該系統(tǒng)的傳遞函數(shù)。傳遞函數(shù)與脈沖響應(yīng)函數(shù)為一對(duì)Laplace變換對(duì)，從而得到：

(9)

式(5)中的水動(dòng)力參數(shù)A(ω)和B(ω)能夠通過(guò)SESAM軟件提取，由于海上浮式結(jié)構(gòu)的水動(dòng)力參數(shù)無(wú)法通過(guò)解析的方式得到解析解，只能通過(guò)數(shù)值方法計(jì)算其離散值。式(5)和式(9)可以通過(guò)關(guān)系式s=jω建立聯(lián)系，根據(jù)離散值在頻域內(nèi)估算浮式結(jié)構(gòu)的傳遞函數(shù)表達(dá)式，即：

(10)

式中：

(11)

由式(10)可知傳遞函數(shù)的有理分式形式是關(guān)于θ的函數(shù)，求解式(10)中θ的常用方法包括非線性最小二乘擬合法[14]、擬線性頻域回歸法[15]及權(quán)重迭代頻域擬合法[16]。這幾種方法在實(shí)際應(yīng)用中有各自的優(yōu)缺點(diǎn)，文中采用權(quán)重迭代法對(duì)延遲函數(shù)的Laplace變換進(jìn)行頻域線性回歸擬合，其計(jì)算公式為：

(12)

其中，

(13)

通過(guò)幾次迭代后，式(12)會(huì)很快收斂，即θk≈θk-1，從而估算出式(10)中有理分式的分子和分母系數(shù)，即得到系統(tǒng)傳遞函數(shù)表達(dá)式。

2.2 狀態(tài)空間模型構(gòu)建

傳遞函數(shù)的有理分式可以等價(jià)地轉(zhuǎn)化為極值—留數(shù)的和式形式，即：

(14)

式中：λik,l為該輸入輸出系統(tǒng)的極值，γik,l為對(duì)應(yīng)的留數(shù)。

式(14)中的極值為分母Q(s)=0的根，可以通過(guò)求解下式獲得。

sn+qn-1sn-1+……+q1s+q0=0

(15)

將式(15)的根s=λik,l代入下列極限公式中，對(duì)應(yīng)的留數(shù)為：

(16)

由式(15)～(16)得到的卷積項(xiàng)對(duì)應(yīng)系統(tǒng)的極值和留數(shù)，從而可以構(gòu)建該系統(tǒng)的狀態(tài)空間模型，即:

(17)

其中，

(18)

2.3 動(dòng)力響應(yīng)計(jì)算

通過(guò)上述方法得到各延遲函數(shù)對(duì)應(yīng)的狀態(tài)空間模型，將各自由度的狀態(tài)空間模型進(jìn)行組裝，得到浮式結(jié)構(gòu)的運(yùn)動(dòng)方程卷積項(xiàng)整體狀態(tài)空間模型，并代入式(7)中，得到：

(19)

其中，Z(t)為組裝后的狀態(tài)向量，并且有：

(20)

其中，A′、B′、C′為Aik、Bik、Cik組裝后的狀態(tài)空間矩陣。

為求解式(19)和(20)，計(jì)算浮式結(jié)構(gòu)的動(dòng)力響應(yīng)，對(duì)式(19)和(20)進(jìn)行變形，并聯(lián)立得到新的狀態(tài)空間方程，即：

(21)

式中：v(t)為浮式結(jié)構(gòu)的速度向量。

采用四階龍格—庫(kù)塔法對(duì)式(21)進(jìn)行積分，進(jìn)而得到浮式結(jié)構(gòu)在外荷載作用的動(dòng)力響應(yīng)。

圖1 浮式海上升壓站Fig. 1 Floating offshore substation

3 數(shù)值算例

采用的數(shù)值模型為日本福島浮式風(fēng)電示范項(xiàng)目的改進(jìn)Spar浮式海上升壓站模型，如圖1所示。該浮式升壓站總高110 m，上部為甲板和塔臺(tái)結(jié)構(gòu)，下部浮式基礎(chǔ)為立柱式平臺(tái)，由中央立柱和三層艙體組成。三層艙體為八邊形柱體，由中部和上部艙體提供浮力，下部艙體壓載。浮式海上升壓站工作水域水深120 m，吃水為50 m。根據(jù)浮式海上升壓站的物理尺寸，采用SESAM軟件建立水動(dòng)力模型，計(jì)算結(jié)構(gòu)的水動(dòng)力參數(shù)，進(jìn)而通過(guò)提出的基于狀態(tài)空間模型動(dòng)力響應(yīng)分析方法計(jì)算浮式結(jié)構(gòu)的動(dòng)力響應(yīng)。

3.1 水動(dòng)力參數(shù)

浮式海上升壓站的附加質(zhì)量A(ω)和興波阻尼B(ω)由SESAM軟件提取。根據(jù)浮式海上升壓站的尺寸可知，該浮式結(jié)構(gòu)的幾何模型同時(shí)關(guān)于xoz平面和yoz平面對(duì)稱，從而使其水動(dòng)力參數(shù)具有特殊的對(duì)稱性質(zhì)。根據(jù)勢(shì)流理論，A(ω)和B(ω)分別對(duì)應(yīng)輻射力的實(shí)部和虛部，從而具有相同的對(duì)稱性質(zhì)，以A(ω)為例，有如下關(guān)系式：

A11(ω)=A22(ω)，A15(ω)=A51(ω)，A24(ω)=A42(ω)，A44(ω)=A55(ω)

(22)

因此，在計(jì)算浮式升壓站的動(dòng)力響應(yīng)時(shí)，只有A11(ω)、A15(ω)、A24(ω)、A33(ω)、A44(ω)和A66(ω)等6個(gè)位置的水動(dòng)力參數(shù)參與計(jì)算。這6個(gè)位置的水動(dòng)力參數(shù)隨入射波頻率變化如圖2所示。

圖2 附加質(zhì)量A(ω)和B(ω)隨頻率變化Fig. 2 Curve of added mass A(ω) and B(ω) with frequency

圖3 估算延遲函數(shù)傅里葉變換與原始值的實(shí)部對(duì)比Fig. 3 Real part comparison between estimated delay function Fourier transform and original value

圖4 估算延遲函數(shù)傅里葉變換與原始值的虛部對(duì)比Fig. 4 Imaginary part comparison between estimated delay function Fourier transform and original value

3.2 浮式升壓站動(dòng)力響應(yīng)分析

圖5 浮式海上升壓站RAOFig. 5 RAO of the floating offshore substation

為了驗(yàn)證方法的正確性，通過(guò)商業(yè)軟件SESAM的Wasim模塊計(jì)算浮式海上升壓站在波浪荷載作用下的動(dòng)力響應(yīng)，將其作為參考值評(píng)估文中方法計(jì)算結(jié)果的正確性與有效性。分析浮式海上升壓站的動(dòng)力響應(yīng)之前，先通過(guò)SESAM計(jì)算浮式結(jié)構(gòu)的RAO。當(dāng)入射波的方向?yàn)?°時(shí)，浮式海上升壓站平臺(tái)的RAO如圖5所示，即縱蕩、垂蕩和橫搖自由度上有響應(yīng)，橫蕩、縱搖和艏搖自由度上的響應(yīng)為零。因此，只對(duì)比浮式海上升壓站的縱蕩、垂蕩和橫搖3個(gè)自由度的動(dòng)力響應(yīng)。

實(shí)際海域中的波浪多為不規(guī)則波，采用Jonswap譜模擬浮式海上升壓站所在海域的不規(guī)則波，其參數(shù)如表1所示。入射波的方向?yàn)?°，從工況1到工況2，海況條件由溫和變惡劣，兩種工況下作用于浮式海上升壓站的波浪力如圖6和7所示，波浪力在橫蕩、縱搖和艏搖自由度上為零。

表1 各工況入射波浪參數(shù)Tab. 1 Incident wave parameters under various working conditions

圖6 工況1波浪力Fig. 6 Wave force of working condition 1

圖7 工況2波浪力Fig. 7 Wave force of working condition 2

通過(guò)SESAM中的Wasim模塊提取浮式海上升壓站的質(zhì)量矩陣M、附加質(zhì)量矩陣A(∞)和靜水恢復(fù)力系數(shù)矩陣C，其結(jié)果如下：

(23)

(24)

(25)

將上式及上節(jié)估算的狀態(tài)空間模型代入式(21)中，并運(yùn)用四階龍格—庫(kù)塔法求解該狀態(tài)空間方程，從而得到浮式海上升壓站在波浪荷載作用下的動(dòng)力響應(yīng)。在工況1和工況2條件下，采用文中方法計(jì)算得到浮式海上升壓站縱蕩、垂蕩和橫搖自由度的響應(yīng)與Wasim計(jì)算結(jié)果對(duì)比如圖8和10所示。從對(duì)比結(jié)果可知，在上述兩種工況下，文中基于狀態(tài)空間模型的浮式海上升壓站動(dòng)力響應(yīng)計(jì)算方法得到的結(jié)果與商業(yè)軟件SESAM計(jì)算結(jié)果吻合較好，二者計(jì)算結(jié)果的差異如圖9和11所示，證明文中方法的正確性。

圖8 工況1動(dòng)力響應(yīng)對(duì)比Fig. 8 Dynamic response comparison of working condition 1

圖9 工況1文中方法與WASIM計(jì)算結(jié)果差異Fig. 9 Difference between the method in this paper and Wasim calculation results of working condition 1

圖10 工況2動(dòng)力響應(yīng)對(duì)比Fig. 10 Dynamic response comparison of working condition 2

圖11 工況2文中方法與WASIM計(jì)算結(jié)果差異Fig. 11 Difference between the method in this paper and Wasim calculation results of working condition 2

同時(shí)，采用狀態(tài)空間模型代替Cummins方程中的卷積項(xiàng)，避免時(shí)域積分時(shí)卷積項(xiàng)運(yùn)算消耗大量計(jì)算資源的情況，從而提高了浮式海上升壓站動(dòng)力響應(yīng)分析的效率。

4 結(jié) 語(yǔ)

提出了一種新的浮式海上升壓站動(dòng)力響應(yīng)分析方法，該方法根據(jù)頻域內(nèi)線性回歸擬合，通過(guò)附加質(zhì)量和興波阻尼，估算延遲函數(shù)對(duì)應(yīng)傳遞函數(shù)的有理分式形式，進(jìn)而計(jì)算延遲函數(shù)的極值和留數(shù)，再由極值和留數(shù)構(gòu)建延遲函數(shù)對(duì)應(yīng)的狀態(tài)空間模型，從而使用狀態(tài)空間模型代替Cummins方程中的卷積項(xiàng)，最后通過(guò)四階龍格—庫(kù)塔法計(jì)算浮式海上升壓站的動(dòng)力響應(yīng)。與傳統(tǒng)的時(shí)域積分方法相比，文中方法通過(guò)狀態(tài)空間模型代替Cummins方程卷積項(xiàng)，不同于SESAM中的頻域方程和直接時(shí)域方法，為間接時(shí)域方法的延伸，避免了卷積項(xiàng)計(jì)算導(dǎo)致的誤差積累，同時(shí)提高了動(dòng)力響應(yīng)分析效率。

文中以日本福島浮式海上風(fēng)電示范項(xiàng)目的浮式升壓站為研究對(duì)象，提取其附加質(zhì)量和興波阻尼，估算延遲函數(shù)的狀態(tài)空間模型，從而計(jì)算浮式升壓站在兩種不規(guī)則波作用下的動(dòng)力響應(yīng)，并與SESAM軟件計(jì)算結(jié)果對(duì)比，二者結(jié)果吻合較好，說(shuō)明文中方法的正確性。