海上浮式風(fēng)機(jī)的一體化建模及其整體可靠性

陳建兵，宋玉鵬

(1. 同濟(jì)大學(xué) 土木工程防災(zāi)國家重點(diǎn)實(shí)驗(yàn)室，上海 200092；2. 同濟(jì)大學(xué) 土木工程學(xué)院，上海 200092；3. 南京工業(yè)大學(xué) 土木工程學(xué)院，南京 211816)

0 引言

深遠(yuǎn)海風(fēng)能具有資源豐富、風(fēng)力穩(wěn)定、環(huán)境友好等諸多優(yōu)點(diǎn)，近年來在全球范圍內(nèi)獲得了廣泛關(guān)注和積極發(fā)展[1]。為了降低深遠(yuǎn)海風(fēng)能開發(fā)成本，發(fā)電機(jī)組采用漂浮式基礎(chǔ)結(jié)構(gòu)已成為共識[2]。當(dāng)前，深遠(yuǎn)海風(fēng)能開發(fā)主要面臨機(jī)組服役條件十分惡劣和由于離岸較遠(yuǎn)而造成高昂維修保養(yǎng)費(fèi)用等挑戰(zhàn)[3-4]。此外，與陸上風(fēng)能和近海風(fēng)能開發(fā)相比，深遠(yuǎn)海風(fēng)能開發(fā)中支撐結(jié)構(gòu)的費(fèi)用急劇增長[1]。上述因素使得當(dāng)前浮式風(fēng)機(jī)尚未完全實(shí)現(xiàn)商業(yè)化，限制了深遠(yuǎn)海風(fēng)能的大規(guī)模開發(fā)。為此，應(yīng)當(dāng)采用基于整體可靠度的結(jié)構(gòu)設(shè)計方法，實(shí)現(xiàn)兼顧浮式風(fēng)機(jī)結(jié)構(gòu)安全性與經(jīng)濟(jì)性的定量設(shè)計，為深遠(yuǎn)海風(fēng)能的大規(guī)模開發(fā)提供關(guān)鍵技術(shù)保障。因此，浮式風(fēng)機(jī)結(jié)構(gòu)的整體可靠度分析具有重要意義。

當(dāng)前，復(fù)雜結(jié)構(gòu)的整體可靠度評估通常采用數(shù)值分析方法，主要涉及三個方面的問題：結(jié)構(gòu)響應(yīng)分析的數(shù)值模型、隨機(jī)荷載的建模、結(jié)構(gòu)可靠度分析方法[5]。浮式風(fēng)機(jī)結(jié)構(gòu)的精細(xì)化數(shù)值分析模型是其整體可靠度分析的關(guān)鍵基礎(chǔ)[6]。由于浮式風(fēng)機(jī)結(jié)構(gòu)中同時存在多種復(fù)雜耦合效應(yīng)[7-8]，例如槳葉-機(jī)艙-塔筒-浮體-錨索等部件之間的相互作用、剛體運(yùn)動與彈性變形耦合、氣動-結(jié)構(gòu)-機(jī)電-水動耦合、槳葉轉(zhuǎn)動效應(yīng)和幾何非線性效應(yīng)等，使得浮式風(fēng)機(jī)結(jié)構(gòu)的一體化動力學(xué)建模成為一項(xiàng)難題[9]。雖然針對該建模問題當(dāng)前已有較多研究，但是在實(shí)際工程中仍主要采用分離式設(shè)計方法，難以合理考慮上述復(fù)雜的耦合效應(yīng)，也難以進(jìn)行結(jié)構(gòu)的整體可靠度分析。一般來說，浮式風(fēng)機(jī)結(jié)構(gòu)的一體化建模主要基于多體動力學(xué)理論[9]，采用Lagrange方法或者Kane方法建立結(jié)構(gòu)的整體動力學(xué)方程，并可以進(jìn)一步分為多剛體動力學(xué)模型和多柔性體動力學(xué)模型[10]。多剛體動力學(xué)模型將浮式風(fēng)機(jī)中的所有部件均考慮為剛體，因而無法獲得柔性部件（槳葉、塔筒等）的彈性變形和內(nèi)力響應(yīng)，具有較大的局限性[10]。但是，由于多剛體動力學(xué)模型建模簡便、分析效率較高，因此多用于浮式風(fēng)機(jī)結(jié)構(gòu)的概念設(shè)計和初始設(shè)計階段[11]。多柔性體模型能夠同時考慮系統(tǒng)中的剛體運(yùn)動和彈性變形響應(yīng)，其中一般采用有限元方法或者模態(tài)疊加方法建立柔性部件的彈性模型[7]。雖然多柔性體模型能夠全面反映系統(tǒng)的力學(xué)特性和復(fù)雜的耦合效應(yīng)，但是該類模型建模較為復(fù)雜，同時存在分析效率不高的問題，難以應(yīng)用于需要進(jìn)行大量數(shù)值模擬分析的問題中，例如結(jié)構(gòu)的可靠度評估與優(yōu)化設(shè)計等[12]。

浮式風(fēng)機(jī)在服役期間的主要環(huán)境荷載為風(fēng)荷載和波浪荷載，且兩種荷載具有相關(guān)性[13]。因此，在浮式風(fēng)機(jī)結(jié)構(gòu)的動力分析與可靠度評估中應(yīng)當(dāng)合理考慮風(fēng)和波浪荷載的聯(lián)合作用[14]。實(shí)測數(shù)據(jù)表明，10 min內(nèi)的風(fēng)速時程和20 min至3～6 h內(nèi)的波浪時程均可以認(rèn)為是平穩(wěn)過程，并可以采用平均風(fēng)速、有效波高、譜峰周期、風(fēng)-浪夾角等物理量描述上述風(fēng)況和海況[13]。為了刻畫結(jié)構(gòu)在服役期內(nèi)所經(jīng)歷的所有工況的概率大小，可基于結(jié)構(gòu)所在場址的長期觀測數(shù)據(jù)建立上述物理量的聯(lián)合概率分布，這也是當(dāng)前海洋結(jié)構(gòu)設(shè)計規(guī)范中所推薦的方法[13]。在建立風(fēng)-浪參數(shù)的聯(lián)合概率分布模型時，可采用條件概率分布建模方法[15]、Nataf變換方法[16]和Copula方法[17-18]。條件概率分布建模方法是當(dāng)前海洋結(jié)構(gòu)設(shè)計規(guī)范中推薦使用的方法[13]。該方法首先確定主導(dǎo)隨機(jī)變量，并建立其邊緣概率分布模型，然后根據(jù)主導(dǎo)隨機(jī)變量的取值將數(shù)據(jù)劃分為一系列區(qū)間，進(jìn)而估計主導(dǎo)隨機(jī)變量取不同值條件下的其他隨機(jī)變量的概率分布。但是，當(dāng)隨機(jī)變量的個數(shù)較多時，該方法需要進(jìn)行多次條件概率分布建模。此外，在某些條件下（例如高風(fēng)速條件），觀測數(shù)據(jù)較少，將導(dǎo)致該條件下的概率分布建模結(jié)果誤差較大[18]。Nataf變換方法首先基于等概率變換原則，將原始隨機(jī)變量變換為服從標(biāo)準(zhǔn)高斯分布的隨機(jī)變量，進(jìn)而利用多維高斯分布模型建立變換后的隨機(jī)變量的聯(lián)合概率分布，最后再利用反變換獲得原始隨機(jī)變量的聯(lián)合概率分布[16]。該方法雖然簡便易行，但是當(dāng)隨機(jī)變量為非線性、非高斯相關(guān)時，例如風(fēng)-浪參數(shù)間的相關(guān)結(jié)構(gòu)通常呈現(xiàn)這種特點(diǎn)，該方法將引入較大誤差。Copula方法是近年來發(fā)展的多維隨機(jī)變量聯(lián)合概率分布建模的新方法，對于非線性相關(guān)問題十分有效[17]。該方法在聯(lián)合概率分布建模中，將隨機(jī)變量的邊緣概率分布建模與變量間的概率相關(guān)結(jié)構(gòu)建模分開處理，具有簡便、靈活、精確等優(yōu)點(diǎn)。該方法目前已在多個領(lǐng)域獲得廣泛關(guān)注[19]，但在海洋工程中尚應(yīng)用較少。

結(jié)構(gòu)可靠度分析方法主要包括一次和二次可靠度方法、基于跨越過程理論的首次超越破壞可靠度方法、隨機(jī)模擬方法、概率密度演化理論等[5,20]。一次和二次可靠度方法常用于靜力問題的分析，具有較高的效率，然而該方法需要將系統(tǒng)的極限狀態(tài)函數(shù)展開為一階或二階形式，對于強(qiáng)非線性系統(tǒng)（例如浮式風(fēng)機(jī)結(jié)構(gòu)）會產(chǎn)生較大誤差[9]?；诳缭竭^程理論的首次超越破壞可靠度方法適用于動力系統(tǒng)，但通常引入Poisson假定，具有一定的近似性，尤其對于大失效概率問題其誤差較大[5]。隨機(jī)模擬方法主要包括Monte Carlo模擬及其改進(jìn)方法以及子集模擬等方法，這類方法不受問題的維度、非線性程度等因素的限制，應(yīng)用范圍較廣，但通常具有隨機(jī)收斂和效率不高等問題[20]。為了將這類方法應(yīng)用于浮式風(fēng)機(jī)的可靠度分析中，通常發(fā)展浮式風(fēng)機(jī)結(jié)構(gòu)分析的代理模型以提高效率，但同時也引入了新的誤差[21]。概率密度演化理論給出了隨機(jī)動力系統(tǒng)響應(yīng)的概率密度函數(shù)的演化方程，揭示了系統(tǒng)響應(yīng)隨機(jī)性演化的內(nèi)在物理規(guī)律[5]，該方法在分析系統(tǒng)的可靠度時同樣不受問題類型、問題維度和非線性程度等因素的影響，能給出系統(tǒng)響應(yīng)量的概率密度函數(shù)演化過程，同時具有較高的分析效率，因此適用于浮式風(fēng)機(jī)的可靠度分析，但目前在此領(lǐng)域的應(yīng)用研究較少。

為了分析浮式風(fēng)機(jī)結(jié)構(gòu)的整體可靠度，首先引入結(jié)合多剛體動力學(xué)理論和有限元方法建立起來的浮式風(fēng)機(jī)結(jié)構(gòu)一體化分析模型（StoDRAFOWT模型）。該模型不僅能同時獲得浮式風(fēng)機(jī)的剛體運(yùn)動和彈性變形響應(yīng)，而且具有建模簡便、易于求解等優(yōu)點(diǎn)[9]。然后，基于我國南海某場址長期風(fēng)-浪數(shù)據(jù)并引入Copula方法，建立風(fēng)-浪聯(lián)合作用的多參數(shù)非線性相關(guān)的聯(lián)合概率分布模型，進(jìn)而采用概率密度演化理論分析浮式風(fēng)機(jī)結(jié)構(gòu)的整體可靠度。由此，建立了風(fēng)-浪聯(lián)合作用下浮式風(fēng)機(jī)結(jié)構(gòu)整體可靠度分析的基本框架。針對風(fēng)機(jī)結(jié)構(gòu)設(shè)計中需要考慮的額定風(fēng)速和極端風(fēng)剪不利工況，采用上述方法進(jìn)行了單柱式浮式風(fēng)機(jī)結(jié)構(gòu)的整體可靠度分析。

1 浮式風(fēng)機(jī)一體化動力學(xué)建模

1.1 基本建模思想與坐標(biāo)系統(tǒng)

如前所述，浮式風(fēng)機(jī)的多剛體動力學(xué)模型雖然分析效率較高，但是無法分析結(jié)構(gòu)的彈性變形。而多柔性體動力學(xué)模型雖然能夠全面分析結(jié)構(gòu)的動力響應(yīng)，但是建模復(fù)雜、分析效率不高。本文首先介紹將多剛體動力學(xué)模型與有限元方法相結(jié)合的建模策略，不僅能夠全面分析結(jié)構(gòu)響應(yīng)，而且建模較為簡便，分析效率較高[9]。已有研究表明[10-11]，對于當(dāng)前兆瓦級的浮式風(fēng)機(jī)系統(tǒng)，結(jié)構(gòu)的彈性變形對系統(tǒng)剛體運(yùn)動響應(yīng)的影響可以忽略，因此可以首先采用多剛體動力學(xué)模型分析系統(tǒng)的剛體運(yùn)動響應(yīng)。而系統(tǒng)的剛體運(yùn)動對結(jié)構(gòu)的彈性變形具有重要影響。為此，采用Kane方法[22-23]推導(dǎo)了考慮系統(tǒng)剛體運(yùn)動效應(yīng)的結(jié)構(gòu)有限元模型，從而能夠全面獲得結(jié)構(gòu)的動力響應(yīng)。

圖1為單柱式浮式風(fēng)機(jī)的示意圖和坐標(biāo)系的定義。為描述浮式風(fēng)機(jī)系統(tǒng)的狀態(tài)，定義三組笛卡爾坐標(biāo)系，分別為慣性坐標(biāo)系Ox0y0z0、浮體-塔筒的局部坐標(biāo)系Oxyz、 槳葉的局部坐標(biāo)系Ox′y′z′。 其中Ox0y0z0坐標(biāo)系的原點(diǎn)位于海平面位置，x0軸與風(fēng)速方向一致，z0軸豎直向上。在初始時刻浮式風(fēng)機(jī)未發(fā)生位移時，坐標(biāo)系Oxyz與Ox0y0z0重 合。坐標(biāo)系Ox′y′z′的原點(diǎn)位于輪轂中心，x′軸與x軸 平行，z′軸由葉根指向葉尖。

圖1 單柱式浮式風(fēng)機(jī)示意圖與坐標(biāo)系統(tǒng)Fig. 1 Schematic diagram of a spar-type floating offshore wind turbine and the coordinate system

1.2 一體化結(jié)構(gòu)分析模型

為了建立浮式風(fēng)機(jī)系統(tǒng)的多剛體動力學(xué)模型，通?；贚agrange方程得到考慮轉(zhuǎn)子轉(zhuǎn)動和浮體六自由度剛體運(yùn)動的動力方程[10]，即包含七個自由度。由于現(xiàn)代風(fēng)力機(jī)通常配備槳距控制系統(tǒng)，為此，將變槳距控制方程與上述七自由度動力方程耦合，得到考慮變速、變槳距控制的浮式風(fēng)機(jī)多剛體動力學(xué)模型[9,24]：

式中，M0和C0分別為浮式風(fēng)機(jī)的質(zhì)量矩陣和阻尼矩陣，其表達(dá)式可參見文獻(xiàn)[24]；Irx為轉(zhuǎn)子相對于x′軸的轉(zhuǎn)動慣量；vp和 ωp=(ωx,ωy,ωz)T分 別為浮體在Oxyz坐標(biāo)系中的平動和轉(zhuǎn)動速度； ωb和 φ分別為槳葉的轉(zhuǎn)動速度和轉(zhuǎn)動角度；β為槳葉的槳距角；c0、G、 τd和 τp是與風(fēng)力機(jī)槳距控制和轉(zhuǎn)矩控制有關(guān)的參數(shù)； ω0為低轉(zhuǎn)速軸的額定轉(zhuǎn)速；FG為重力在Oxyz坐標(biāo)系下的效應(yīng)；MG為 在Oxyz坐標(biāo)系下重力關(guān)于坐標(biāo)系原點(diǎn)的力矩；FD為 外部荷載在Oxyz坐 標(biāo)系下的效應(yīng)；MD為在Oxyz坐標(biāo)系下外部荷載關(guān)于坐標(biāo)系原點(diǎn)的力矩；FM為低轉(zhuǎn)速軸上的凈扭矩。

通過求解方程式（1）即可獲得浮式風(fēng)機(jī)中浮體的六自由度剛體運(yùn)動響應(yīng)、槳葉轉(zhuǎn)速和槳葉的槳距角。為了分析浮式風(fēng)機(jī)中柔性部件的彈性變形，基于Kane方程分別推導(dǎo)得到塔筒單元和槳葉單元的有限元表達(dá)式。

在坐標(biāo)系Oxyz中，可以得到考慮剛體運(yùn)動效應(yīng)的塔筒單元的動力方程[9]：

式中：Mt和 μ分別為單元的質(zhì)量矩陣和結(jié)點(diǎn)位移向量；Ct、Kt和Ft分別為單元的剛度矩陣、阻尼矩陣和荷載向量，表達(dá)式為：

式中：ρ為單元的材料密度；V為單元體積；N?為形函數(shù)矩陣；上標(biāo)“T”表示轉(zhuǎn)置矩陣； ω?為浮體轉(zhuǎn)動速度ωp對 應(yīng)的反對稱矩陣[24]；R為坐標(biāo)轉(zhuǎn)換矩陣，將坐標(biāo)系Oxyz中的矢量轉(zhuǎn)換至坐標(biāo)系Ox0y0z0中 ；ro為Oxyz的坐標(biāo)系原點(diǎn)在坐標(biāo)系Ox0y0z0中 的位矢；rk為單元中任意一點(diǎn)在坐標(biāo)系Oxyz中的位矢。

類似地，在坐標(biāo)系Ox′y′z′中，可以得到考慮剛體運(yùn)動和旋轉(zhuǎn)效應(yīng)的槳葉單元的動力方程[9]：

式中：Mb為 單元的質(zhì)量矩陣；Cb、Kb和Fb分別為單元的剛度矩陣、阻尼矩陣和荷載向量，其表達(dá)式與式（3）、式（4）、式（5）類似。其中，由于剛體運(yùn)動和旋轉(zhuǎn)效應(yīng)產(chǎn)生的槳葉單元的附加阻尼矩陣、附加剛度矩陣和附加荷載向量的表達(dá)式為：

其中：T為坐標(biāo)轉(zhuǎn)換矩陣，將坐標(biāo)系Ox′y′z′中的矢量轉(zhuǎn)換至坐標(biāo)系Oxyz中；為 槳葉轉(zhuǎn)動速度 ( ωb,0,0)T對應(yīng)的 反 對 稱 矩 陣； ¨ro為Oxyz的 坐 標(biāo) 系 原 點(diǎn) 在 坐 標(biāo) 系Ox0y0z0中 的加速度；ro′為Ox′y′z′的坐標(biāo)系原點(diǎn)在坐標(biāo)系Oxyz中的位矢；rq為 單元中任意一點(diǎn)在坐標(biāo)系Ox′y′z′中的位矢。

值得注意的是，槳葉截面為非對稱形狀，并且槳葉截面存在沿翼展方向變化的預(yù)扭轉(zhuǎn)角，同時還存在槳距控制。因此，槳葉在風(fēng)輪平面內(nèi)、外兩個方向的振動相互耦合，在建模過程中應(yīng)當(dāng)合理考慮。此外，槳葉在風(fēng)機(jī)正常運(yùn)轉(zhuǎn)工況下的變形較大，應(yīng)當(dāng)在槳葉的三維梁單元模型中考慮幾何非線性效應(yīng)。這是與塔筒有限元模型的不同之處。

基于塔筒和槳葉的有限元模型，并結(jié)合前述浮式風(fēng)機(jī)的剛體運(yùn)動響應(yīng)，可建立浮式風(fēng)機(jī)柔性構(gòu)件的整體有限元模型，進(jìn)而分析浮式風(fēng)機(jī)結(jié)構(gòu)彈性變形響應(yīng)。

1.3 荷載計算

在浮式風(fēng)機(jī)結(jié)構(gòu)的動力響應(yīng)分析中，主要應(yīng)考慮槳葉和塔筒的氣動荷載、浮體的水動力荷載和錨泊荷載。本文采用葉素動量理論方法[25]計算槳葉的氣動荷載，并考慮了Prandtl葉尖損失因子和高推力系數(shù)條件下的Glauert修正[25]。塔筒的氣動荷載可基于準(zhǔn)定常假定計算[25]。

在計算大型浮體的水動力荷載時，應(yīng)當(dāng)合理考慮輻射效應(yīng)和繞射效應(yīng)[26]。本文關(guān)注的是單柱式浮式風(fēng)機(jī)，浮體的截面尺寸相對較小，因此可以簡化考慮這兩種效應(yīng)。在實(shí)際中可采用相對形式的Morison公式計算單柱式浮體的水動力荷載[10,26]，此時，單位長度浮體的水動力荷載為[27]：

式中：ρw為 海水密度；dp為浮體外徑；Cd和CA分別為拖曳力系數(shù)和附加質(zhì)量系數(shù)；vf和af分別為該浮體微元的速度和加速度向量；vw和aw分別為水質(zhì)點(diǎn)的速度和加速度，可采用Airy波浪理論、Stokes波浪理論或其他波浪理論 計 算[13,25]； ∥·∥表示 向 量 的幅值。應(yīng) 指 出，求解結(jié)構(gòu)的動力方程（式1）時，通常將式（13）中的移至動力方程的左邊，由此產(chǎn)生浮體的附加質(zhì)量矩陣。

在計算浮式風(fēng)機(jī)的錨泊荷載時，需要建立錨索的力學(xué)分析模型，主要包括錨索的靜力模型、擬靜力模型和動力模型[28]。由于靜力和擬靜力模型在某些工況下會低估錨索的疲勞荷載，因此推薦采用錨索的動力模型[28]。本文采用錨索的集中質(zhì)量動力模型[28]計算錨泊荷載。該模型不僅考慮錨索的受拉彈性變形，而且將錨索的動力效應(yīng)和水動力效應(yīng)以及與海底的接觸效應(yīng)等也考慮在內(nèi)。錨索底部固定在海床上，錨索頂部的邊界條件由浮體的位置確定。通過求解錨索的動力方程，可以獲得錨索對浮體的荷載效應(yīng)。

1.4 模型驗(yàn)證

為了檢驗(yàn)上述模型的正確性，分別利用該模型和當(dāng)前風(fēng)力機(jī)結(jié)構(gòu)分析的常用軟件FAST，分析風(fēng)-浪聯(lián)合作用下的海上單柱式浮式風(fēng)機(jī)結(jié)構(gòu)的動力響應(yīng)，并進(jìn)行對比。浮式風(fēng)機(jī)結(jié)構(gòu)參數(shù)采用美國國家可再生能源實(shí)驗(yàn)室（NREL）給出的5 MW海上風(fēng)力發(fā)電機(jī)組參數(shù)[29]和OC3項(xiàng)目中海上單柱式浮式風(fēng)機(jī)基礎(chǔ)參數(shù)[30]，關(guān)鍵參數(shù)如表1所示。

表1 浮式風(fēng)力發(fā)電機(jī)組關(guān)鍵參數(shù)Table 1 Key parameters of the floating offshore wind turbine

本文僅給出一種工況條件下的分析對比結(jié)果，其他工況條件下的對比結(jié)果見文獻(xiàn)[24]。此處的工況條件為：定常風(fēng)速場，且輪轂高度處的風(fēng)速為15 m/s；規(guī)則波浪，波高為6 m，波浪周期為10 s。分析結(jié)果如圖2所示。結(jié)果對比表明，本文模型的分析結(jié)果與FAST的分析結(jié)果吻合良好，尤其是在響應(yīng)進(jìn)入平穩(wěn)階段以后。因此，可以利用本文模型進(jìn)行浮式風(fēng)機(jī)結(jié)構(gòu)的可靠度分析。

圖2 風(fēng)-浪作用下浮式風(fēng)機(jī)結(jié)構(gòu)響應(yīng)的對比結(jié)果Fig. 2 Comparisons of the responses of floating offshore wind turbine under wind and wave excitations

2 風(fēng)-浪聯(lián)合作用的概率模型

由于浮式風(fēng)機(jī)在服役期間同時受到風(fēng)荷載和波浪荷載作用，因此在進(jìn)行浮式風(fēng)機(jī)的動力響應(yīng)與可靠度分析時，應(yīng)當(dāng)合理考慮風(fēng)-浪聯(lián)合作用。如前所述，當(dāng)前海洋結(jié)構(gòu)的設(shè)計規(guī)范中建議采用風(fēng)-浪參數(shù)的聯(lián)合概率分布模型考慮風(fēng)-浪聯(lián)合作用[13]。本文采用Copula方法建立風(fēng)-浪參數(shù)的聯(lián)合概率分布模型。

在實(shí)際工程中，通常需要考慮平均風(fēng)速Vw、有效波高Hs、 波浪譜峰周期Tp以及風(fēng)-浪夾角等物理量。本文主要關(guān)注額定風(fēng)速條件下發(fā)生極端風(fēng)剪時浮式風(fēng)機(jī)結(jié)構(gòu)的可靠度。根據(jù)風(fēng)力機(jī)結(jié)構(gòu)的設(shè)計規(guī)范，在分析該工況下風(fēng)力機(jī)結(jié)構(gòu)的動力響應(yīng)時，取風(fēng)-浪夾角為0°[31]。因此，需要建立額定風(fēng)速條件下Hs和Tp的條件聯(lián)合概率分布模型。

在建立風(fēng)-浪參數(shù)的聯(lián)合概率分布模型時，應(yīng)當(dāng)基于場址的長期風(fēng)-浪觀測數(shù)據(jù)。本文采用歐洲中尺度天氣預(yù)報中心ERA-Interim數(shù)據(jù)庫[32]提供的風(fēng)-浪參數(shù)再分析數(shù)據(jù)，獲取了我國南海場址（21°N，113°E）處1979年至2018年（共計40年）的風(fēng)-浪再分析數(shù)據(jù)。這些數(shù)據(jù)均為30 min內(nèi)觀測數(shù)據(jù)的平均值，時間分辨率為6 h。通過對數(shù)據(jù)預(yù)處理，可得到輪轂處平均風(fēng)速、有效波高和譜峰周期的長期數(shù)據(jù)。

為了建立額定風(fēng)速條件下Hs和Tp的條件聯(lián)合概率分布模型，在上述風(fēng)-浪數(shù)據(jù)的基礎(chǔ)上選取額定風(fēng)速附近的數(shù)據(jù)進(jìn)行概率建模。本文中，風(fēng)力機(jī)的額定風(fēng)速為11.4 m/s，因此取風(fēng)速在11.3 m/s ～11.5 m/s區(qū)間內(nèi)的風(fēng)-浪數(shù)據(jù)進(jìn)行Hs和Tp的聯(lián)合概率分布建模，共有681個數(shù)據(jù)點(diǎn)。

根據(jù)Copula聯(lián)合概率分布建模的基本理論[17]，對于任意兩個隨機(jī)變量X1和X2，其聯(lián)合概率密度函數(shù)的表達(dá)式為：

式中：F1(x1) 和F2(x2) 分 別 為X1和X2的邊緣概率分布函數(shù)；f1(x1)和f2(x2)為相應(yīng)的邊緣概率密度函數(shù)；c12(·,·)為二元Copula密度函數(shù)，可從常用的Copula函數(shù)中選用，例如高斯Copula、t Copula、Clayton Copula、Gumbel Copula、Frank Copula等[17]。

從式（14）中可以看到，在利用Copula方法建立多維隨機(jī)變量的聯(lián)合概率分布模型時，需要分別建立各隨機(jī)變量的邊緣概率分布模型及變量間的相關(guān)結(jié)構(gòu)模型。由于邊緣概率分布模型及變量間的相關(guān)結(jié)構(gòu)模型可以分開處理，因此該方法具有靈活、便捷、精確等優(yōu)點(diǎn)。

基于上述681個額定風(fēng)速條件下的數(shù)據(jù)點(diǎn)，首先采用常用的參數(shù)概率分布模型，建立Hs和Tp的邊緣概率分布模型，并進(jìn)行假設(shè)檢驗(yàn)，進(jìn)而基于AIC準(zhǔn)則[19]確定最優(yōu)的邊緣概率分布模型。結(jié)果表明，Hs和Tp的最優(yōu)邊緣分布均為Gamma分布，結(jié)果見表2和圖3。表中， Γ (·)為Gamma函數(shù)。

圖3 額定風(fēng)速下有效波高和譜峰周期的邊緣概率分布Fig. 3 Marginal probability distributions of Hs and Tp under the rated wind speed condition

進(jìn)一步，采用常用的Copula函數(shù)建立Hs和Tp的相關(guān)結(jié)構(gòu)模型。備選的Copula函數(shù)除了上述常用的Copula函數(shù)外，還包括按照如下方式構(gòu)建的非對稱Copula函數(shù)[33]：

在確定備選的Copula函數(shù)后，采用極大似然方法估算Copula函數(shù)中的待定參數(shù)，然后根據(jù)AIC準(zhǔn)則選取最優(yōu)Copula函數(shù)并進(jìn)行擬合優(yōu)度檢驗(yàn)[17，19]。對于上述額定風(fēng)速條件下的風(fēng)浪數(shù)據(jù)，識別的最優(yōu)結(jié)果為按照式（15）構(gòu)造的非對稱高斯Copula模型，見表2。表中，C(u,v) 表示高斯Copula， Φλ(·)為相關(guān)系數(shù)為 λ(? 1≤λ≤1) 的二維標(biāo)準(zhǔn)正態(tài)分布函數(shù)， Φ(·)為一維標(biāo)準(zhǔn)正態(tài)分布函數(shù)， Φ?1(·) 為 Φ(·)的反函數(shù)。進(jìn)而，基于式（14），可以得到額定風(fēng)速條件下Hs和Tp的聯(lián)合概率密度函數(shù)，如圖4所示。

圖4 額定風(fēng)速下有效波高和譜峰周期的聯(lián)合概率密度函數(shù)Fig. 4 Joint probability density function of Hs and Tp under the rated wind speed condition

表2 波浪參數(shù)概率分布建模結(jié)果Table 2 Probabilistic modeling results of wave parameters

應(yīng)指出，本文采用Copula方法建立了二維隨機(jī)環(huán)境變量的條件聯(lián)合概率分布模型。對于高維隨機(jī)環(huán)境變量，可采用C-vine Copula方法便捷地建立其聯(lián)合概率分布模型[18-19]。

3 浮式風(fēng)機(jī)結(jié)構(gòu)整體可靠性分析

3.1 概率密度演化理論簡介

本文采用概率密度演化理論進(jìn)行額定風(fēng)速下發(fā)生極端風(fēng)剪時浮式風(fēng)機(jī)的可靠度分析。為此，先簡要介紹概率密度演化理論。

不失一般性，對于多自由度動力系統(tǒng)，考慮系統(tǒng)的隨機(jī)因素后，動力方程可寫為：

式中：Y(t)為 系統(tǒng)的狀態(tài)響應(yīng)量；Y0為初始時刻的系統(tǒng)狀態(tài)； Θ=(Θ1,...,Θs)為s維隨機(jī)向量，刻畫了系統(tǒng)中結(jié)構(gòu)參數(shù)和外部激勵的隨機(jī)性。 Θ的聯(lián)合概率密度函數(shù)記為pΘ(θ)。

式(16)的解可寫為Y(t)=H(Θ,t)，該系統(tǒng)中其他物理量的響應(yīng)Z(t)=(Z1(t),...,Zm(t))均 可通過Y(t)和(t)獲得?；诟怕适睾阍?，李杰和陳建兵[5]建立了廣義概率密度演化方程。當(dāng)只關(guān)注一個物理響應(yīng)量Z(t)時，廣義概率密度演化方程具有如下形式[5]：

且具有如下的初值條件：

式 中：pZΘ(z,θ,t)為 (Θ,Z(t))的 聯(lián) 合 概 率 密 度 函 數(shù)；δ(·)為 Dirac函數(shù)；z0為Z(t)的初值。

因此，Z(t)的概率密度函數(shù)為：

式中， ?Θ為 隨機(jī)向量 Θ的分布域。

為了求解廣義概率密度演化方程，已發(fā)展了點(diǎn)演化和群演化兩類數(shù)值求解方法。本文采用點(diǎn)演化方法，具體的求解步驟可參考文獻(xiàn)[5]。

基于概率密度演化理論發(fā)展了兩種結(jié)構(gòu)首次超越破壞可靠度分析方法，即，基于吸收邊界條件的可靠度分析方法和基于極值分布的可靠度分析方法[5]。本文采用基于極值分布的可靠度分析方法。首先構(gòu)造隨機(jī)過程Z(Θ,t)的等價極值事件：

式中，T為結(jié)構(gòu)響應(yīng)的分析時間。

進(jìn)而構(gòu)造如下的虛擬隨機(jī)過程[5]：

由于Z(Θ,t) 的 極值分布等價于t=1時 刻Zvir(Θ,t)的概率分布，因此可先通過求解廣義概率密度演化方程獲得Zvir(Θ,t)的 概率密度演化過程，然后獲得Z(Θ,t)的極值分布，即：

由此，可得到結(jié)構(gòu)的可靠度為：

式中，b為Z(Θ,t)的失效閾值。

3.2 算例

針對風(fēng)力機(jī)結(jié)構(gòu)設(shè)計規(guī)范[31]中定義的額定風(fēng)速條件下發(fā)生極端風(fēng)剪的工況，進(jìn)行了單柱式浮式風(fēng)機(jī)結(jié)構(gòu)的整體可靠度分析。此處僅考慮豎直方向的正極端風(fēng)剪，此時風(fēng)剖面為[31]：

式中：Vhub為輪轂高度處的平均風(fēng)速；z和zhub分別為高度和輪轂高度；α=0.2為風(fēng)剖面指數(shù)律的冪指數(shù)；T0=12s為 極端風(fēng)剪的持續(xù)時間；t為 時間；D為風(fēng)輪直徑； β=6.4 ； σ1為 脈 動 風(fēng) 速 標(biāo) 準(zhǔn) 差 的 代 表 值，σ1=Iref(0.75Vhub+5.6)， 其中Iref為 湍 流 強(qiáng) 度； Λ1為湍流積分尺度； σ1和 Λ1的取值可根據(jù)規(guī)范[31]確定。當(dāng)Vhub=15m/s、Iref=0.12 時，5 MW風(fēng)力機(jī)在t=T0/2時的極端風(fēng)剖面和正常風(fēng)剖面對比圖見圖5。

圖5 正常風(fēng)剖面和極端風(fēng)剖面對比Fig. 5 Comparison of normal and extreme wind profiles

在分析結(jié)構(gòu)的可靠度時，應(yīng)當(dāng)充分、合理考慮系統(tǒng)中的隨機(jī)源。本文采用的概率密度演化理論在分析結(jié)構(gòu)可靠度時能夠同時全面考慮結(jié)構(gòu)自身和外部激勵的隨機(jī)性[5]。作為示例，本文不僅將額定風(fēng)速條件下的有效波高Hs和 譜峰周期Tp作為隨機(jī)變量，而且考慮了脈動風(fēng)速、波面高程以及塔筒阻尼比 ξT的隨機(jī)性。本文基于譜表達(dá)方法并分別采用波數(shù)-頻率聯(lián)合功率譜[34-35]和JONSWAP譜[13]模擬隨機(jī)脈動風(fēng)場和隨機(jī)波浪場。為了降低隨機(jī)場模擬中的隨機(jī)變量（相位）個數(shù)，采用演化相位譜方法[36]，由此分別引入4個和1個相位演化時間，分別記為Te1～Te4和Te5。相位演化時間為隨機(jī)變量，基于實(shí)測風(fēng)速數(shù)據(jù)和波浪數(shù)據(jù)可以識別其概率分布[36]。相位演化時間和塔筒阻尼比的概率分布信息如表3所示[9]。Hs和Tp的概率分布信息由第2節(jié)方法確定。

表3 隨機(jī)變量的概率分布信息Table 3 Probability distributions of random variables

在求解廣義概率密度演化方程時，首先需要將概率空間剖分為Nq個子域，并在每個子域內(nèi)選取一個代表點(diǎn)[5]。然后對每個代表點(diǎn)進(jìn)行確定性動力響應(yīng)分析以獲得概率密度演化方程（偏微分方程）中的時變系數(shù)，進(jìn)而求解概率密度演化方程，并獲得響應(yīng)量的概率密度演化結(jié)果[5]。本文共考慮8個隨機(jī)變量，采用基于點(diǎn)集GF-偏差最小的點(diǎn)集優(yōu)選策略[37]在八維概率空間中選取了329個代表點(diǎn)。圖6中的紅點(diǎn)展示的是Hs和Tp的代表點(diǎn)，圖中的灰線為Hs和Tp聯(lián)合概率密度函數(shù)的等值線。

圖6 額定風(fēng)速下隨機(jī)動力分析中的波浪參數(shù)代表點(diǎn)Fig. 6 Representative points of wave parameters for the stochastic analysis under rated wind speed condition

在選定代表點(diǎn)之后，需要模擬每種情況下的隨機(jī)風(fēng)場和波浪場，并分析浮式風(fēng)機(jī)結(jié)構(gòu)的動力響應(yīng)。作為示例，圖7給出了Hs=2m、Tp=8.4s工況下浮式風(fēng)機(jī)結(jié)構(gòu)動力響應(yīng)的典型樣本。同時分析了不考慮極端風(fēng)剪時的結(jié)構(gòu)響應(yīng)并進(jìn)行對比。動力響應(yīng)的時間分析長度為1500 s，并從第500 s開始施加極端風(fēng)剪條件。為了避免初始瞬態(tài)響應(yīng)的影響，僅展示后1200 s的動力激勵和結(jié)構(gòu)動力響應(yīng)結(jié)果。圖中，“EWS”表示極端風(fēng)剪。

圖7 浮式風(fēng)機(jī)結(jié)構(gòu)在特定工況下的響應(yīng)Fig. 7 Responses of the floating offshore wind turbine under specified load conditions

從圖7（a）可以看出，在施加極端風(fēng)剪條件后，z=30m 處和z=150m處的風(fēng)速差異明顯增大，由此導(dǎo)致塔筒和槳葉的位移響應(yīng)發(fā)生變化（如圖7d、圖7f所示）。從圖7（d）還可以看出，考慮極端風(fēng)剪后，該樣本中塔頂位移響應(yīng)的極值有所增大，因此會影響結(jié)構(gòu)的整體可靠度結(jié)果。圖7（c）和圖7（f）顯示，大約在700～750 s時間段內(nèi)，結(jié)構(gòu)響應(yīng)有突然的“跌落”，這是由于啟動了變槳控制，降低了氣動荷載（如圖7（j）所示）。圖7（g）～圖7（j）表明，極端風(fēng)剪對浮式風(fēng)機(jī)的縱蕩響應(yīng)、縱搖響應(yīng)、槳葉轉(zhuǎn)速、槳葉的槳距角變化均影響較小。

采用概率密度演化理論可以獲得結(jié)構(gòu)響應(yīng)極值的概率分布。以塔頂順風(fēng)向位移為例，通過求解廣義概率密度演化方程，得到其響應(yīng)極值的概率密度函數(shù)如圖8所示。在給定的失效閾值條件下，當(dāng)塔頂位移響應(yīng)超過該閾值時，即認(rèn)為發(fā)生失效。由此可以得到不同失效閾值條件下的失效概率，如圖9所示。

圖8 塔頂順風(fēng)向位移響應(yīng)極值的概率密度函數(shù)Fig. 8 Probability density function of the extreme value of the fore-aft displacement response at the tower top

圖9 失效概率隨失效閾值變化關(guān)系Fig. 9 Relationship between failure probability and failure threshold

作為示例，本文在分析浮式風(fēng)機(jī)結(jié)構(gòu)的可靠度時，僅考慮了塔頂位移。在實(shí)際中，應(yīng)考慮多個指標(biāo)，例如塔頂加速度、塔底應(yīng)力以及浮式風(fēng)機(jī)的縱搖響應(yīng)等。概率密度演化理論能夠同時考慮多個失效指標(biāo)[5]。當(dāng)同時采用多個失效指標(biāo)分析時，可靠度結(jié)果可能發(fā)生變化。

4 結(jié)論

浮式風(fēng)機(jī)結(jié)構(gòu)的整體可靠性分析是進(jìn)行基于可靠度優(yōu)化設(shè)計的關(guān)鍵基礎(chǔ)，對于大規(guī)模深遠(yuǎn)海風(fēng)能的開發(fā)具有重要意義。本文針對海上單柱式浮式風(fēng)機(jī)結(jié)構(gòu)，首先介紹了一體化耦合動力學(xué)分析模型，然后建立了環(huán)境參數(shù)的聯(lián)合概率分布模型以合理考慮風(fēng)-浪聯(lián)合作用，進(jìn)而分析浮式風(fēng)機(jī)結(jié)構(gòu)在極端工況條件下的動力可靠度。主要研究結(jié)論如下：

1）對于當(dāng)前兆瓦級浮式風(fēng)機(jī)，柔性部件的彈性變形對浮式風(fēng)機(jī)的剛體運(yùn)動響應(yīng)幾乎沒有影響，因此在建模過程中可以忽略這一因素，從而使建模方法與分析過程更為簡便。本文基于這一思想，將多剛體動力學(xué)理論與有限元方法結(jié)合起來考察了浮式風(fēng)機(jī)結(jié)構(gòu)的一體化模型，不僅能夠彌補(bǔ)一般多剛體動力學(xué)模型的不足，而且利用Kane方法便捷地考慮了系統(tǒng)中剛體運(yùn)動對柔性部件變形的影響。雖然從理論上看，全耦合多柔性體模型（例如FAST）的建模方法考慮得更為全面，但是從實(shí)際工程應(yīng)用的角度來看，本文所提的方法更為便捷。兩種模型的對比分析結(jié)果表明了本文模型具有良好的精度。

2）在浮式風(fēng)機(jī)結(jié)構(gòu)的動力響應(yīng)分析中，應(yīng)當(dāng)合理考慮風(fēng)-浪聯(lián)合作用?；趫鲋返拈L期風(fēng)-浪觀測數(shù)據(jù)并采用Copula方法，能夠精確、便捷地建立風(fēng)-浪多參數(shù)非線性相關(guān)的聯(lián)合概率分布模型，從而合理地考慮風(fēng)和波浪對結(jié)構(gòu)的荷載效應(yīng)。

3）基于概率密度演化理論的可靠度分析方法適用于海上大型浮式風(fēng)機(jī)結(jié)構(gòu)的整體可靠性分析。該方法不僅能夠同時考慮外部激勵和結(jié)構(gòu)自身的隨機(jī)性，而且能夠獲得結(jié)構(gòu)響應(yīng)的概率密度演化過程，同時具有很高的分析效率。

本文研究建立了海上大型浮式風(fēng)機(jī)結(jié)構(gòu)整體可靠性分析的基本框架，為進(jìn)行基于整體可靠度的結(jié)構(gòu)優(yōu)化設(shè)計奠定了基礎(chǔ)。由于我國東南海域的風(fēng)電場常年遭受臺風(fēng)侵襲[38]，因此在未來工作中需要進(jìn)一步分析在臺風(fēng)等極端環(huán)境條件下的浮式風(fēng)機(jī)結(jié)構(gòu)整體可靠度，并開展基于整體可靠度的結(jié)構(gòu)優(yōu)化設(shè)計等研究。