基于本征正交分解的波浪力隨機(jī)場(chǎng)降維模擬

劉章軍， 劉 磊， 汪 峰

(1. 武漢工程大學(xué) 土木工程與建筑學(xué)院，武漢 430074；2. 三峽大學(xué) 土木與建筑學(xué)院，湖北 宜昌 443002)

隨著我國(guó)經(jīng)濟(jì)快速發(fā)展，海洋資源的開發(fā)日益重要，對(duì)海洋工程結(jié)構(gòu)(如海上平臺(tái)、離岸碼頭和海上風(fēng)電機(jī)裝置等)安全性的要求越來(lái)越高。此類結(jié)構(gòu)下部支撐一般為墩柱結(jié)構(gòu)，在復(fù)雜的海洋環(huán)境中，將受到波浪荷載、風(fēng)荷載、輪渡撞擊和地震等作用。其中，波浪荷載是主要作用力之一，其不僅隨時(shí)間變化(具有動(dòng)態(tài)特性)，而且具有明顯的隨機(jī)性和空間相干性，因而采用波浪力隨機(jī)場(chǎng)加以描述[1]。

在波浪隨機(jī)過(guò)程或隨機(jī)場(chǎng)研究中，Pierson[2]建立了波面位移隨機(jī)過(guò)程與波高譜之間的關(guān)系，開辟了應(yīng)用波高譜研究波浪隨機(jī)過(guò)程的先河。Borgman[3]對(duì)Morison方程進(jìn)行了線性化處理，并利用譜分析法研究了作用在柱樁上的隨機(jī)波浪力。Di Paola等[4]結(jié)合波高與波速的關(guān)系，利用本征正交分解(POD)方法模擬了波速隨機(jī)場(chǎng)。Ren等[5]進(jìn)行了結(jié)構(gòu)在隨機(jī)波浪沖擊下的數(shù)值模擬分析，并與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比。Dong等[6]結(jié)合波浪譜模擬了重力式網(wǎng)箱在隨機(jī)海浪力作用下的力學(xué)特性。上述方法均是基于波高譜來(lái)實(shí)現(xiàn)波浪隨機(jī)過(guò)程或隨機(jī)場(chǎng)模擬。在波浪方向譜研究方面，Longuet-higgins等[7]應(yīng)用傅里葉級(jí)數(shù)估計(jì)海浪方向譜，并提出了方向譜的分布函數(shù)。Ren等[8]對(duì)方向譜的觀測(cè)方法進(jìn)行了改進(jìn)，并進(jìn)行了波面模擬；文獻(xiàn)[9-11]利用方向譜對(duì)波浪進(jìn)行了三維數(shù)值模擬；趙振民等[12]則將波面的三維仿真模擬用于船舶振動(dòng)研究。

上述波浪隨機(jī)過(guò)程或隨機(jī)場(chǎng)模擬均是以Monte Carlo方法為基礎(chǔ)，為保證模擬精度，往往需要對(duì)一系列隨機(jī)變量進(jìn)行大量的隨機(jī)抽樣，不僅極大地增加了波浪隨機(jī)場(chǎng)模擬的計(jì)算量，而且對(duì)復(fù)雜海洋工程結(jié)構(gòu)的非線性隨機(jī)動(dòng)力反應(yīng)分析帶來(lái)巨大困難。為克服這一困難，文獻(xiàn)[13-14]通過(guò)引入隨機(jī)函數(shù)的約束，實(shí)現(xiàn)了僅用一個(gè)基本隨機(jī)變量即可精確地模擬一維單變量隨機(jī)過(guò)程。為此，本文將隨機(jī)函數(shù)思想引入到平穩(wěn)多變量隨機(jī)過(guò)程模擬的POD方法中，結(jié)合線性化的Morison方程和P-M波浪譜，實(shí)現(xiàn)用兩個(gè)基本隨機(jī)變量即可模擬作用于小尺度樁上的波浪力隨機(jī)場(chǎng)，且僅需數(shù)百條代表性樣本即可在概率密度層次上描述波浪力隨機(jī)場(chǎng)的概率特性，為結(jié)合概率密度演化理論[15-16]研究海洋工程結(jié)構(gòu)在隨機(jī)波浪力場(chǎng)作用下的動(dòng)力響應(yīng)與可靠度分析奠定基礎(chǔ)。

1 基于正交隨機(jī)變量的本征正交分解

本征正交分解是時(shí)-空隨機(jī)場(chǎng)模擬的最主要方法之一。在一般情況下，由于無(wú)法直接獲取特征問(wèn)題的解析解或封閉解，需要將連續(xù)的1V-2D時(shí)-空隨機(jī)場(chǎng)f0(x,t)離散化為nV-1D隨機(jī)過(guò)程X(t)。

假設(shè)X(t)=[X1(t),X2(t),…,Xn(t)]T是一個(gè)零均值的nV-1D平穩(wěn)隨機(jī)過(guò)程，其雙邊的功率譜密度函數(shù)矩陣SX(ω)可表示為

(1)

式中：Sii(ω)是第i個(gè)分量過(guò)程Xi(t)的自功率譜密度函數(shù)；Sij(ω)(i≠j)是第i個(gè)分量過(guò)程Xi(t)與第j個(gè)分量過(guò)程Xj(t)的互功率譜密度函數(shù)。

一般地，功率譜密度函數(shù)矩陣SX(ω)是一個(gè)非負(fù)定的Hermitian矩陣， 其特征值λi(ω)是頻率ω的非負(fù)實(shí)函數(shù)， 特征向量ψi(ω)=[φ1i(ω),φ2i(ω),…,φni(ω)]T一般是頻率ω的復(fù)函數(shù)。對(duì)于任意給定的頻率ω，存在如下的特征問(wèn)題

Ψ*T(ω)Ψ(ω)=In×n,SX(ω)Ψ(ω)=Ψ(ω)Λ(ω)

(2)

式中： 符號(hào)“*”和“T”分別表示共軛和轉(zhuǎn)置；Λ(ω)=diag{λ1(ω),λ2(ω),…,λn(ω)}為特征值矩陣；Ψ(ω)=[ψ1(ω),ψ2(ω),…,ψn(ω)]為特征函數(shù)矩陣。

對(duì)于零均值的nV-1D平穩(wěn)隨機(jī)過(guò)程X(t)， 可分解為n個(gè)互不相干的子隨機(jī)向量過(guò)程之和[17]

(3)

式中： 子隨機(jī)向量過(guò)程xi(t)的表達(dá)式為

(4)

在式(4)中，Zi(ω)(i=1,2,…,n)是零均值的復(fù)隨機(jī)過(guò)程，其增量滿足如下的基本條件

E[dZi(ω)]=0

(5a)

(5b)

(5c)

式中：δij和δωω′均為Kronecker符號(hào)。

進(jìn)一步地，在式(4)中，若令

(6)

式中：Pik是一組零均值的標(biāo)準(zhǔn)正交復(fù)隨機(jī)變量，滿足如下的基本條件

(7)

于是，子隨機(jī)向量過(guò)程xi(t)可近似地表達(dá)為有限級(jí)數(shù)的復(fù)形式

(8)

式中：ωk為離散的頻率點(diǎn)，N為頻率截?cái)囗?xiàng)數(shù)， Δω為頻率步長(zhǎng)。為了盡量避免零頻率點(diǎn)對(duì)模擬結(jié)果的影響，本文定義ωk如下

(9)

(10)

式中：Rik和Iik是一組零均值的正交隨機(jī)變量，滿足如下的基本條件

E[Rik]=E[Iik]=0

(11a)

E[RikIjm]=0

(11b)

(11c)

最后，將式(10)代入式(3)中得到平穩(wěn)多變量隨機(jī)過(guò)程X(t)的本征正交分解

(12)

式(12)即為基于正交隨機(jī)變量的本征正交分解。 在實(shí)際應(yīng)用中，由于正交隨機(jī)變量Rik和Iik僅需滿足基本條件式(11)， 其概率分布并未給定，因而不能直接采用式(12)來(lái)模擬平穩(wěn)多變量隨機(jī)過(guò)程。

2 基于本征正交分解的降維模擬

在基于正交隨機(jī)變量的本征正交分解式(12)中，若將正交隨機(jī)變量Rik和Iik定義為

Rik=cosφik,Iik=sinφik

(13)

式中：φik(i=1,2,…,n,k=1,2,…,N)為區(qū)間[0,2π)上相互獨(dú)立的均勻分布隨機(jī)相位角。顯然，式(13)滿足正交隨機(jī)變量的基本條件式(11)。

同時(shí)，將特征向量ψi(ωk)的實(shí)部ηi(ωk)和虛部ρi(ωk)分別寫為

ηi(ωk)=|ψi(ωk)|cos ?i(ωk)

(14a)

ρi(ωk)=|ψi(ωk)|sin ?i(ωk)

(14b)

于是，將式(13)和式(14)代入式(12)中，得到基于隨機(jī)相位角的本征正交分解模擬式[18]

(15)

式中：X(C)(t)為傳統(tǒng)的本征正交分解模擬過(guò)程。在式(15)中，由于隨機(jī)相位角的概率分布已知，通過(guò)隨機(jī)相位角的一組隨機(jī)抽樣即可生成樣本函數(shù)，因而在工程模擬中得到廣泛應(yīng)用。

從上述式(15)的推導(dǎo)過(guò)程可知，基于隨機(jī)相位角的本征正交分解式(15)是基于正交隨機(jī)變量的本征正交分解式(12)的一個(gè)特例。事實(shí)上，通過(guò)定義正交隨機(jī)變量Rik和Iik為式(13)的形式，使得式(15)中隨機(jī)變量的數(shù)量比式(12)減少一半。這表明，通過(guò)施加某種約束(即定義式(13))可以有效地減少平穩(wěn)隨機(jī)向量過(guò)程模擬的隨機(jī)度(隨機(jī)變量的數(shù)量)。

進(jìn)一步，在基于正交隨機(jī)變量的本征正交分解式(12)中，可利用隨機(jī)函數(shù)對(duì)正交隨機(jī)變量Rik和Iik施加更強(qiáng)的約束，從而實(shí)現(xiàn)平穩(wěn)隨機(jī)向量過(guò)程的降維模擬。為此，將滿足基本條件式(11)的正交隨機(jī)變量Rik和Iik均定義為r維隨機(jī)變量Θ=(Θ1,Θ2,…,Θr)的正交函數(shù)形式， 即Rik=hik(Θ)，Iik=gik(Θ)， 其中hik(·)和gik(·)均為確定性的正交函數(shù)。于是，式(11)可改寫為

(16a)

(16b)

(16c)

(16d)

(16e)

式中：i,j=1,2,…,n和k,m=1,2,…,N；Ω為r維隨機(jī)變量Θ的定義區(qū)間；pΘ(θ)為r維隨機(jī)向量Θ的聯(lián)合概率密度。

(17)

3 基于線性化Morison方程的波浪力隨機(jī)場(chǎng)

在隨機(jī)波浪下，水質(zhì)點(diǎn)的水平速度u和水平加速度a可分別表示為

(18)

(19)

ω2=gktanh (kh)

(20)

式中：g為重力加速度，可取g=9.81 m/s2。

對(duì)于直立的小直徑圓柱樁，可采用Morison公式來(lái)計(jì)算隨機(jī)波浪力[20]

F(z,t)=FD(z,t)+FI(z,t)

(21)

式中：FD(z,t)和FI(z,t)分別為拖曳力和慣性力，其表達(dá)式為

(22a)

(22b)

式中：ρ為海水的密度；D為樁柱的直徑；CD和CI分別為拖曳力系數(shù)和慣性力系數(shù)，文獻(xiàn)[21]建議CD=1.2及CI=2.0。

(23)

式中：σu(z)為水質(zhì)點(diǎn)水平速度u(z,t)的均方差。 對(duì)于海浪隨機(jī)過(guò)程η(t)采用P-M譜時(shí)，σu(z)可近似寫為[22]

(24)

式中：H1/3為三分之一有效波高。

于是，隨機(jī)波浪力的等價(jià)線性表達(dá)式為

(25)

在式(23)中， 對(duì)于任意點(diǎn)z=zi處的拖曳力隨機(jī)過(guò)程FD(zi,t)和z=zj處的隨機(jī)過(guò)程FD(zj,t)， 其功率譜密度函數(shù)可寫為

SDij(ω)=Sη(ω)Gi(ω)Gj(ω),i,j=1,2,…,n

(26)

式中：SDij(ω)為隨機(jī)過(guò)程FD(zi,t)和FD(zj,t)的功率譜密度，Sη(ω)為海浪隨機(jī)過(guò)程η(t)的功率譜密度，Gi(ω)的表達(dá)式為

(27)

因此，n維拖曳力平穩(wěn)隨機(jī)向量過(guò)程FD(t)=[FD(z1,t),FD(z2,t),…,FD(zn,t)]T的功率譜密度矩陣為

(28)

(29)

式中：Li(ω)的表達(dá)式為

(30)

最后，n維波浪力平穩(wěn)隨機(jī)向量過(guò)程F(t)=[F(z1,t),F(z2,t),…,F(zn,t)]T的功率譜密度矩陣可表示為

SF(ω)=Sη(ω)×

(31)

在上述Morison公式中，得到的拖曳力和慣性力的功率譜密度矩陣是一個(gè)非負(fù)定的Hermitian矩陣，滿足本征正交分解的要求。

4 計(jì)算實(shí)例

為了驗(yàn)證本文方法的有效性，采用文獻(xiàn)[23]中直立小圓柱樁進(jìn)行分析，取樁底為原點(diǎn)，海底平面水平向?yàn)閤軸，豎直方向?yàn)閦軸，模擬10個(gè)等間距點(diǎn)處的波浪力隨機(jī)過(guò)程，并與Monte Carlo方法的模擬結(jié)果進(jìn)行比較，波浪力沿z軸的分布如圖1所示。

圖1 波浪力沿z軸的分布

在此算例中， 波高譜采用P-M譜

(32)

式中：g=9.81 m/s2為重力加速度；v為海面以上19.5 m處的風(fēng)速，取v=16.24 m/s。

本文算例的具體參數(shù)及取值如表1所示。

表1 計(jì)算參數(shù)及取值

圖2為本文方法與Monte Carlo法生成位置z3、z7和z9三點(diǎn)處(見圖1所示)的波浪力代表性樣本曲線(第100條)的比較。由圖可知，代表性樣本曲線具有平穩(wěn)隨機(jī)波浪過(guò)程的基本特性，且本文方法與Monte Carlo法生成的曲線在幅值、波長(zhǎng)變化等方面具有相似性。

表2為采用本文方法和Monte Carlo方法分別生成10個(gè)點(diǎn)處波浪力樣本集合的均值及標(biāo)準(zhǔn)差相對(duì)誤差的平均值，其中，各點(diǎn)處波浪力樣本集合的均值相對(duì)誤差和標(biāo)準(zhǔn)差相對(duì)誤差的定義見文獻(xiàn)[16]。表3為采用本文方法和Monte Carlo方法分別生成10個(gè)點(diǎn)處波浪力樣本集合所需的總時(shí)間。從表2中可知，Monte Carlo方法生成的233條和1 000條樣本的均值相對(duì)誤差均遠(yuǎn)大于本文方法生成的233條代表性樣本均值相對(duì)誤差；對(duì)于標(biāo)準(zhǔn)差相對(duì)誤差，本文方法生成233條代表性樣本的標(biāo)準(zhǔn)差相對(duì)誤差雖略大于Monte Carlo方法生成1 000條的樣本標(biāo)準(zhǔn)差相對(duì)誤差，但其較小于Monte Carlo方法生成233條樣本的標(biāo)準(zhǔn)差相對(duì)誤差。同時(shí)，從表3中可知，Monte Carlo方法生成233條樣本的總時(shí)間略小于本文方法生成233條代表性樣本總時(shí)間，但其生成1 000條樣本的總時(shí)間約為本文方法生成233條代表性樣本總時(shí)間的3倍。因此，從模擬的精度和效率兩方面綜合考慮，本文方法生成的233條代表性樣本遠(yuǎn)優(yōu)于Monte Carlo方法生成的233條樣本，也優(yōu)于Monte Carlo方法生成的1 000條樣本。此外，本文方法僅需兩個(gè)基本隨機(jī)變量即可精細(xì)地表達(dá)波浪力隨機(jī)場(chǎng)，生成的代表性樣本數(shù)量少，且構(gòu)成一個(gè)完備的概率集，因而與概率密度演化理論具有一致性，這為結(jié)合概率密度演化理論進(jìn)行復(fù)雜工程結(jié)構(gòu)的非線性隨機(jī)動(dòng)力響應(yīng)分析和可靠度評(píng)價(jià)提供了基礎(chǔ)。

圖2 隨機(jī)波浪力代表性樣本

圖3和圖4分別給出了本文方法在z3和z7兩點(diǎn)處生成的233條波浪力代表性樣本的自相關(guān)函數(shù)、互相關(guān)函數(shù)與目標(biāo)值比較。為了清晰地反映模擬值與目標(biāo)值的擬合程度，圖3和圖4均僅給出了時(shí)間區(qū)間(-50 s～50 s)的比較結(jié)果，同時(shí)對(duì)自相關(guān)函數(shù)和互相關(guān)函數(shù)進(jìn)行了歸一化處理。從圖3和圖4中可知，點(diǎn)z3和z7處模擬的自相關(guān)函數(shù)、互相關(guān)函數(shù)與目標(biāo)值均擬合一致，進(jìn)一步表明了本文方法的有效性。

表2 均值和標(biāo)準(zhǔn)差的相對(duì)誤差

表3 模擬效率比較

圖3 自相關(guān)函數(shù)比較

圖4 互相關(guān)函數(shù)的比較

5 結(jié) 論

利用波高譜，結(jié)合線性化的Morison方程，建立了波浪力隨機(jī)場(chǎng)的功率譜密度函數(shù)矩陣。在基于正交隨機(jī)變量的本征正交分解上，引入正交隨機(jī)變量的隨機(jī)函數(shù)形式，實(shí)現(xiàn)了波浪力隨機(jī)場(chǎng)模擬的降維表達(dá)。數(shù)值算例表明，本文方法具有如下特點(diǎn)：

(1) 結(jié)合波高譜與Morison方程，構(gòu)造了波浪力多變量隨機(jī)過(guò)程。與直接積分法計(jì)算整個(gè)樁柱上的波浪力相比，本文方法可精細(xì)化模擬作用在小尺度樁上每一點(diǎn)處的波浪力隨機(jī)過(guò)程，更符合實(shí)際情況。

(2) 基于本征正交分解的降維模擬方法，能夠極大地降低波浪力多變量隨機(jī)過(guò)程模擬的隨機(jī)度(隨機(jī)變量的數(shù)量)。在本文方法中，僅需2個(gè)基本隨機(jī)變量即可模擬波浪力多變量隨機(jī)過(guò)程，從而可利用數(shù)論方法選取基本隨機(jī)變量的代表性點(diǎn)集，實(shí)現(xiàn)僅用數(shù)百條代表性樣本即可在概率密度層次上描述波浪力多變量隨機(jī)過(guò)程的概率特性。

(3) 與Monte Carlo方法相比，本文方法所需的基本隨機(jī)變量最少，生成的代表性樣本數(shù)量少，且每條代表性樣本具有給定的賦得概率，所有的代表性樣本構(gòu)成一個(gè)完備概率集，從而可結(jié)合概率密度演化理論進(jìn)行復(fù)雜海洋工程結(jié)構(gòu)隨機(jī)波浪力作用的非線性動(dòng)力反應(yīng)分析。