劉章軍， 劉 磊， 汪 峰

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

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

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

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

1 基于正交隨機變量的本征正交分解

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

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

(1)

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

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

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

(2)

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

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

(3)

式中： 子隨機向量過程xi(t)的表達式為

(4)

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

E[dZi(ω)]=0

(5a)

(5b)

(5c)

式中：δij和δωω′均為Kronecker符號。

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

(6)

式中：Pik是一組零均值的標準正交復隨機變量，滿足如下的基本條件

(7)

于是，子隨機向量過程xi(t)可近似地表達為有限級數(shù)的復形式

(8)

式中：ωk為離散的頻率點，N為頻率截斷項數(shù)， Δω為頻率步長。為了盡量避免零頻率點對模擬結果的影響，本文定義ωk如下

(9)

(10)

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

E[Rik]=E[Iik]=0

(11a)

E[RikIjm]=0

(11b)

(11c)

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

(12)

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

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

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

Rik=cosφik,Iik=sinφik

(13)

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

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

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

(14a)

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

(14b)

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

(15)

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

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

進一步，在基于正交隨機變量的本征正交分解式(12)中，可利用隨機函數(shù)對正交隨機變量Rik和Iik施加更強的約束，從而實現(xiàn)平穩(wěn)隨機向量過程的降維模擬。為此，將滿足基本條件式(11)的正交隨機變量Rik和Iik均定義為r維隨機變量Θ=(Θ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維隨機變量Θ的定義區(qū)間；pΘ(θ)為r維隨機向量Θ的聯(lián)合概率密度。

(17)

3 基于線性化Morison方程的波浪力隨機場

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

(18)

(19)

ω2=gktanh (kh)

(20)

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

對于直立的小直徑圓柱樁，可采用Morison公式來計算隨機波浪力[20]

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

(21)

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

(22a)

(22b)

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

(23)

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

(24)

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

于是，隨機波浪力的等價線性表達式為

(25)

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

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

(26)

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

(27)

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

(28)

(29)

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

(30)

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

SF(ω)=Sη(ω)×

(31)

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

4 計算實例

為了驗證本文方法的有效性，采用文獻[23]中直立小圓柱樁進行分析，取樁底為原點，海底平面水平向為x軸，豎直方向為z軸，模擬10個等間距點處的波浪力隨機過程，并與Monte Carlo方法的模擬結果進行比較，波浪力沿z軸的分布如圖1所示。

圖1 波浪力沿z軸的分布

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

(32)

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

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

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

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

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

圖2 隨機波浪力代表性樣本

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

表2 均值和標準差的相對誤差

表3 模擬效率比較

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

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

5 結 論

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

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

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

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