一種基于狀態空間模型的浮式海上升壓站平臺動力響應計算方法研究

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

(1. 浙江省深遠海風電技術研究重點實驗室，浙江 杭州 311122； 2. 中國電建集團華東勘測設計研究院有限公司，浙江 杭州 311122； 3. 三峽新能源山東昌邑發電有限公司，山東 濰坊 261300)

隨著海上風電建設的推進，海上風電逐漸由近海走向深海，為了適應海上風電的發展趨勢，海上風電場中的海工結構物結構型式也由固定式轉變為浮式。作為電力輸送、匯集的樞紐，浮式海上升壓站的安全性對海上風電場的健康運營至關重要。相比于固定式升壓站，浮式海上升壓站在6自由度上的響應更為明顯，這也給升壓站中的電氣設備帶來較大的安全隱患[1]。因此，需要在設計階段對浮式海上升壓站結構進行動力響應分析，進而優化其水動力性能，保障浮式海上升壓站的平穩運行。

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

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

1 基礎理論

1.1 浮式平臺運動方程

在勢流理論的假設條件下，只考慮一階水動力作用，無航速海上浮式結構的運動方程可以用Cummins方程表示：

(1)

1.2 Oglivie關系式

Cuminns方程中附加質量A和延遲函數K(t)是與海上浮式結構水動力參數相關的量，為探討它們之間的數學關系，對式(1)進行傅里葉變換[12]：

(2)

Ogilvie通過對船體運動預測的研究，建立了隨頻率變化的附加質量A(ω)及興波阻尼B(ω)與延遲函數K(t)的關系式，即[10]：

(3)

(4)

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

(5)

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

(6)

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

(7)

2 狀態空間模型頻域算法

2.1 頻域線性回歸擬合

(8)

對式(8)進行Laplace變換，即可得到該系統的傳遞函數。傳遞函數與脈沖響應函數為一對Laplace變換對，從而得到：

(9)

式(5)中的水動力參數A(ω)和B(ω)能夠通過SESAM軟件提取，由于海上浮式結構的水動力參數無法通過解析的方式得到解析解，只能通過數值方法計算其離散值。式(5)和式(9)可以通過關系式s=jω建立聯系，根據離散值在頻域內估算浮式結構的傳遞函數表達式，即：

(10)

式中：

(11)

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

(12)

其中，

(13)

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

2.2 狀態空間模型構建

傳遞函數的有理分式可以等價地轉化為極值—留數的和式形式，即：

(14)

式中：λik,l為該輸入輸出系統的極值，γik,l為對應的留數。

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

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

(15)

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

(16)

由式(15)～(16)得到的卷積項對應系統的極值和留數，從而可以構建該系統的狀態空間模型，即:

(17)

其中，

(18)

2.3 動力響應計算

通過上述方法得到各延遲函數對應的狀態空間模型，將各自由度的狀態空間模型進行組裝，得到浮式結構的運動方程卷積項整體狀態空間模型，并代入式(7)中，得到：

(19)

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

(20)

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

為求解式(19)和(20)，計算浮式結構的動力響應，對式(19)和(20)進行變形，并聯立得到新的狀態空間方程，即：

(21)

式中：v(t)為浮式結構的速度向量。

采用四階龍格—庫塔法對式(21)進行積分，進而得到浮式結構在外荷載作用的動力響應。

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

3 數值算例

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

3.1 水動力參數

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

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

(22)

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

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

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

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

3.2 浮式升壓站動力響應分析

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

為了驗證方法的正確性，通過商業軟件SESAM的Wasim模塊計算浮式海上升壓站在波浪荷載作用下的動力響應，將其作為參考值評估文中方法計算結果的正確性與有效性。分析浮式海上升壓站的動力響應之前，先通過SESAM計算浮式結構的RAO。當入射波的方向為0°時，浮式海上升壓站平臺的RAO如圖5所示，即縱蕩、垂蕩和橫搖自由度上有響應，橫蕩、縱搖和艏搖自由度上的響應為零。因此，只對比浮式海上升壓站的縱蕩、垂蕩和橫搖3個自由度的動力響應。

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

表1 各工況入射波浪參數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

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

(23)

(24)

(25)

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

圖8 工況1動力響應對比Fig. 8 Dynamic response comparison of working condition 1

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

圖10 工況2動力響應對比Fig. 10 Dynamic response comparison of working condition 2

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

同時，采用狀態空間模型代替Cummins方程中的卷積項，避免時域積分時卷積項運算消耗大量計算資源的情況，從而提高了浮式海上升壓站動力響應分析的效率。

4 結 語

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

文中以日本福島浮式海上風電示范項目的浮式升壓站為研究對象，提取其附加質量和興波阻尼，估算延遲函數的狀態空間模型，從而計算浮式升壓站在兩種不規則波作用下的動力響應，并與SESAM軟件計算結果對比，二者結果吻合較好，說明文中方法的正確性。