隨機(jī)風(fēng)場(chǎng)模擬的連續(xù)本征正交分解-隨機(jī)函數(shù)方法

劉章軍， 劉增輝

(1. 防災(zāi)減災(zāi)湖北省重點(diǎn)實(shí)驗(yàn)室(三峽大學(xué))， 湖北 宜昌 443002； 2. 三峽大學(xué) 土木與建筑學(xué)院， 湖北 宜昌 443002)

在風(fēng)工程中，一般認(rèn)為脈動(dòng)風(fēng)速隨機(jī)場(chǎng)的概率特性在時(shí)間上是不變的，在空間上是變化的，因而可將其看成是一個(gè)典型的非完全均勻時(shí)-空隨機(jī)場(chǎng)。在應(yīng)用Monte Carlo方法模擬時(shí)-空隨機(jī)場(chǎng)時(shí)，主要有譜表示法[1-3]和本征正交分解(Proper Othogonal Decomposition, POD)法[4-6]，以及譜表示與本征正交分解相結(jié)合的混合方法[7]。在具體實(shí)施中，一般將多維單變量(1V-mD)的連續(xù)時(shí)-空隨機(jī)場(chǎng)轉(zhuǎn)化為一維多變量(nV-1D)的離散隨機(jī)場(chǎng)(隨機(jī)向量)過程。值得說明的是，應(yīng)用譜表示法模擬一維多變量隨機(jī)向量過程是基于功率譜密度矩陣的Cholesky分解；而POD法則是基于功率譜密度矩陣或協(xié)方差矩陣的特征分解。僅在特定情況下，即互功率譜密度函數(shù)的特征問題存在封閉的解析解時(shí)，POD法模擬多維單變量隨機(jī)場(chǎng)則是基于互功率譜密度函數(shù)的特征分解[8]。顯然，基于功率譜密度矩陣或協(xié)方差矩陣的本征正交分解是離散形式；而基于互功率譜密度函數(shù)的POD法則是連續(xù)形式，因而可有效地提高互功率譜密度函數(shù)特征分解的效率。

上述模擬方法中，無論是譜表示方法還是POD法，都需要成千上萬個(gè)隨機(jī)變量來實(shí)現(xiàn)對(duì)脈動(dòng)風(fēng)速隨機(jī)場(chǎng)的模擬，從而極大地增加了計(jì)算工作量。為了克服這一局限性，文獻(xiàn)[9]基于物理的建模思想，建立了脈動(dòng)風(fēng)速場(chǎng)的物理隨機(jī)函數(shù)模型，實(shí)現(xiàn)了用若干個(gè)基本隨機(jī)變量表達(dá)脈動(dòng)風(fēng)速場(chǎng)。文獻(xiàn)[10]則基于隨機(jī)函數(shù)的思想，建立了一維單變量隨機(jī)過程的正交展開-隨機(jī)函數(shù)模型，從而實(shí)現(xiàn)了僅用一個(gè)基本隨機(jī)變量對(duì)原隨機(jī)過程在二階統(tǒng)計(jì)意義上的精確模擬。本文進(jìn)一步將隨機(jī)函數(shù)的思想引入到脈動(dòng)風(fēng)速隨機(jī)場(chǎng)中，結(jié)合基于互功率譜密度函數(shù)的POD法[8]，建立脈動(dòng)風(fēng)速隨機(jī)場(chǎng)(1V-2D)的連續(xù)POD-隨機(jī)函數(shù)模型，實(shí)現(xiàn)僅用兩個(gè)基本隨機(jī)變量即可表達(dá)脈動(dòng)風(fēng)速隨機(jī)場(chǎng)的目的。同時(shí)，生成的脈動(dòng)風(fēng)速代表性時(shí)程可構(gòu)成一個(gè)完備的概率集，在本質(zhì)上與概率密度演化理論[11-12]具有統(tǒng)一性，這為應(yīng)用概率密度演化理論進(jìn)行結(jié)構(gòu)風(fēng)振響應(yīng)和抗風(fēng)可靠度分析奠定基礎(chǔ)。

1 基于互功率譜密度函數(shù)的本征正交分解方法

設(shè)f0(x,t)是一個(gè)零均值的1V-2D連續(xù)時(shí)-空隨機(jī)場(chǎng)，其定義在時(shí)間變量t和坐標(biāo)為x的一維空間域D上。假定f0(x,t)是一個(gè)非完全均勻的時(shí)-空隨機(jī)場(chǎng)，即f0(x,t)關(guān)于時(shí)間變量t是平穩(wěn)的，關(guān)于空間變量x是有限能量的。對(duì)于非完全均勻的連續(xù)時(shí)-空隨機(jī)場(chǎng)f0(x,t)，可以表示為Fourier-Stieltjes積分形式

(1)

式中：Z(x,ω)是一個(gè)正交增量的復(fù)隨機(jī)場(chǎng)，其頻率增量dZ(x,ω)=Z(x,ω+dω)-Z(x,ω)滿足如下的條件

E[dZ(x,ω)]=0, dZ(x,-ω)=dZ*(x,ω)

(2a)

E[dZ(x,ω)dZ*(x′,ω′)]=

(2b)

式中:E[·]為數(shù)學(xué)期望；“*”為取共軛復(fù)數(shù)；Sf0(x,x′,ω)為隨機(jī)過程f0(x,t)與f0(x′,t)的雙邊互功率譜密度函數(shù)。設(shè)λk(ω)與ψk(x,ω)(k=1,2,…)分別為互功率譜密度函數(shù)Sf0(x,x′,ω)的特征值和特征函數(shù)，它們是第二類Fredholm積分方程的非平凡解

(3)

注意到，互功率譜密度函數(shù)Sf0(x,x′,ω)是一個(gè)有界的、埃爾米特的、非負(fù)定的函數(shù)，其特征值λk(ω)是非負(fù)的實(shí)函數(shù)，特征函數(shù)ψk(x,ω)一般是頻率ω的復(fù)函數(shù)，且具有如下的正交性

(4)

λi(ω)δij

(5)

式中：δij為Kronecker符號(hào)。事實(shí)上，特征函數(shù)集{ψk(x,ω),k=1,2,…}的完備性保證了互功率譜密度函數(shù)的譜分解

(6)

一般地，Sf0(x,x′,ω)存在有限或無限個(gè)特征值，可將特征值按從大到小的順序排列，并取前n階展開項(xiàng)來近似代替式(6)，其展開精度計(jì)算如下

(7)

式中：Π(n)為展開精度；n為展開項(xiàng)數(shù)。因此，式(6)可近似寫為

(8)

結(jié)合式(8)與式(2)，可知

(9)

(10)

于是，將式(9)代入式(1)中，連續(xù)的時(shí)-空隨機(jī)場(chǎng)f0(x,t)可以表示為

f0(x,t)=

(11)

式(11)即為連續(xù)時(shí)-空隨機(jī)場(chǎng)f0(x,t)的POD形式。POD將連續(xù)時(shí)-空隨機(jī)場(chǎng)f0(x,t)表達(dá)為前n階分量之和的形式，這些分量稱為連續(xù)時(shí)-空隨機(jī)場(chǎng)f0(x,t)的本征模態(tài)，其中特征函數(shù)ψk(x,ω)確定了關(guān)于空間變量x的模態(tài)形狀，特征值λk(ω)則表征了各階模態(tài)的能量。

對(duì)于非完全均勻的連續(xù)時(shí)-空隨機(jī)場(chǎng)f0(x,t)，其數(shù)值模擬可以通過式(11)的頻率離散形式來實(shí)現(xiàn)

f(x,t)=

(12)

式中：f(x,t)為模擬的時(shí)-空隨機(jī)場(chǎng)；ωm=mΔω，ωu=NΔω為截?cái)鄨A頻率；N為頻率截?cái)囗?xiàng)數(shù)；Pkm=Wk(ωm)為一組零均值的正交復(fù)隨機(jī)變量，滿足如下的正交性

(13)

(Rkmcosωmt+Ikmsinωmt)+Zk(x,ωm)×

(Rkmsinωmt-Ikmcosωmt)]

(14)

ψk(x,ωm)=χk(x,ωm)-iZk(x,ωm)

(15a)

Pkm=Rkm-iIkm

(15b)

式中：Rkm和Ikm為零均值的實(shí)正交隨機(jī)變量，滿足如下的基本條件

E[Rkm]=E[Ikm]=0，E[RkmIpq]=0，

(16)

當(dāng)特征函數(shù)ψk(x,ω)(k=1,2,…,n)為實(shí)函數(shù)時(shí)，則式(14)可進(jìn)一步簡(jiǎn)化為

Rkmcosωmt+Ikmsinωmt

(17)

式(14)或式(17)即為基于互功率譜密度函數(shù)的POD模擬公式。

2 正交隨機(jī)變量集的隨機(jī)函數(shù)表達(dá)

(18a)

(18b)

pΘ1(θ1)pΘ2(θ2)dθ1dθ2=0

(18c)

(18d)

(18e)

式中：i,r=1,2,…,n，j,s=1,2,…,N；Ω1,Ω2分別為基本隨機(jī)變量Θ1和Θ2的定義區(qū)間，pΘ1(θ1)和pΘ2(θ2)分別為基本隨機(jī)變量Θ1和Θ2的概率密度函數(shù)。

由式(18)，可定義如下隨機(jī)函數(shù)形式

i=1,2,…,n;j=1,2,…,N

(19)

3 脈動(dòng)風(fēng)場(chǎng)的互功率譜密度函數(shù)及其特征問題的封閉解

一般地，假定脈動(dòng)風(fēng)場(chǎng)是一個(gè)零均值的非完全均勻的連續(xù)隨機(jī)場(chǎng)v(x,t)，定義在0≤x≤L和0

(20)

式中：Sv(x,x′,ω)為隨機(jī)過程v(x,t)和v(x′,t)的雙邊互功率譜密度函數(shù)，m2/s；S0(ω)為隨機(jī)過程v(x,t)的雙邊自功率譜密度函數(shù)，m2/s；c為水平方向上的衰減因子，可取c=10；U為給定地面高度的平均風(fēng)速，m/s；ω為圓頻率，rad/s。

x,x′∈[0,L]

(21)

結(jié)合式(3)和式(21)，互功率譜密度函數(shù)的特征問題可表示為

ψk(x′,ω)dx′

(22)

式(22)特征問題的封閉解為

特征值

(23)

特征函數(shù)

(24)

式中：參數(shù)μk必須滿足如下的條件

(25)

(26)

可見，參數(shù)μk僅依賴于參數(shù)α(ω)；因此，可根據(jù)參數(shù)α(ω)的離散值來求解參數(shù)μk，從而獲得特征值和特征函數(shù)的數(shù)值解。

從特征函數(shù)的解析式(24)可知，特征函數(shù)ψk(x,ω)關(guān)于空間域x∈[0,L]的中點(diǎn)x=L/2具有對(duì)稱性，即：當(dāng)tan(μk/2)=α/μk時(shí)，特征函數(shù)關(guān)于中點(diǎn)x=L/2是正對(duì)稱的；當(dāng)tan(μk/2)=-μk/α?xí)r，特征函數(shù)關(guān)于中點(diǎn)x=L/2是反對(duì)稱的，這與結(jié)構(gòu)動(dòng)力學(xué)中的結(jié)構(gòu)振型具有類似的物理意義。事實(shí)上，特征函數(shù)表征了隨機(jī)風(fēng)場(chǎng)關(guān)于空間域的各階模態(tài)形狀[14]。

4 數(shù)值算例

以Kaimal脈動(dòng)風(fēng)速譜作為算例，對(duì)于給定的地面高度z=zd，其雙邊自功率譜密度函數(shù)的表達(dá)式為

(27)

式中：U為給定地面高度z=zd的平均風(fēng)速；u*為氣流的剪切速度；其計(jì)算公式為

(28)

式中：z0為地面粗糙長(zhǎng)度，本文取z0=0.03 m。

計(jì)算分析中，結(jié)構(gòu)的水平跨徑L=150 m，地面高度zd=35 m，平均風(fēng)速U=45 m/s。此時(shí)，雙邊自功率譜密度函數(shù)可寫為

(29)

為簡(jiǎn)便之，僅以水平坐標(biāo)x=50 m，80 m和120 m三點(diǎn)處的脈動(dòng)風(fēng)速進(jìn)行分析。

在隨機(jī)風(fēng)場(chǎng)的數(shù)值模擬中，取截?cái)囝l率ωu=2π rad/s，頻率步長(zhǎng)Δω=0.01 rad/s，則頻率截?cái)囗?xiàng)數(shù)N=628；持時(shí)T=600 s，時(shí)間步長(zhǎng)Δt=0.1 s。同時(shí)，根據(jù)式(7)可確定POD模擬公式(17)中的展開項(xiàng)數(shù)n，為保證計(jì)算精度，本文要求展開精度Π(n)≥90%。

圖1給出特征值展開項(xiàng)數(shù)與展開精度的關(guān)系，其中圖1(a)為前5階特征值的大小分布；圖1(b)為前20階特征值所對(duì)應(yīng)的展開精度。從圖1可知，第1階特征值占有約60%的能量，且前5階特征值即可使展開精度達(dá)到90.2%，為此取展開項(xiàng)數(shù)n=5。

(a) 前五階特征值

(b) 前20階特征值對(duì)應(yīng)的展開精度

(a) 50 m處的第100條代表性時(shí)程

(b) 80 m處的第200條代表性時(shí)程

(c) 120 m處的第300條代表性時(shí)程

圖3為610條脈動(dòng)風(fēng)速代表性時(shí)程集合自功率譜與目標(biāo)自功率譜的比較。從圖3可知，模擬的自功率譜與目標(biāo)譜在低頻部分?jǐn)M合較好，但隨著頻率增大，其擬合誤差將逐漸增大。事實(shí)上，從特征值的分布圖1(a)可知，前幾階特征值在低頻部分的能量占絕對(duì)優(yōu)勢(shì)；然而，在高頻部分，各階特征值幾乎相等，因而各階模態(tài)所占的能量比較接近，當(dāng)僅采用前幾階本征模態(tài)模擬原脈動(dòng)風(fēng)速隨機(jī)場(chǎng)時(shí)，勢(shì)必將導(dǎo)致高頻部分的誤差較大。

圖3 模擬自功率譜密度函數(shù)與目標(biāo)自功率譜的比較

圖4分別給出了x=50 m與x=80 m處、x=50 m與x=120 m處以及x=80 m與x=120 m處610條脈動(dòng)風(fēng)速代表性時(shí)程的集合互功率譜與目標(biāo)互功率譜的比較。從圖4可知，模擬的互功率譜與目標(biāo)互功率譜擬合較好，進(jìn)一步證明本方法的有效性。

(a) x=50 m與x=80 m

(b) x=50 m與x=120 m

(c) x=80 m與x=120 m

5 結(jié) 論

在基于互功率譜密度函數(shù)的本征正交分解方法上，通過引入正交隨機(jī)變量集的隨機(jī)函數(shù)表達(dá)形式，提出了脈動(dòng)風(fēng)速隨機(jī)場(chǎng)模擬的連續(xù)本征正交分解-隨機(jī)函數(shù)方法。以Kaimal脈動(dòng)風(fēng)速譜為例進(jìn)行了脈動(dòng)風(fēng)速隨機(jī)場(chǎng)的模擬分析。研究表明，本方法具有如下特點(diǎn)：

(1) 基于互功率譜密度函數(shù)的本征正交分解方法具有明確的物理意義，僅用少數(shù)幾階本征模態(tài)將脈動(dòng)風(fēng)速隨機(jī)場(chǎng)表達(dá)為連續(xù)形式，同時(shí)給出了互功率譜密度函數(shù)特征問題的封閉解，實(shí)現(xiàn)了對(duì)脈動(dòng)風(fēng)速隨機(jī)場(chǎng)的高效降階處理。

(2) 基于互功率譜密度函數(shù)的本征正交分解方法，對(duì)于脈動(dòng)風(fēng)速隨機(jī)場(chǎng)，若空間域D定義在迎風(fēng)面的水平方向上，由于水平方向上各點(diǎn)的自功率譜相等，因此本方法適用于各類風(fēng)譜；若空間域D定義在迎風(fēng)面的鉛直方向上，為保證互功率譜密度函數(shù)特征問題存在解析解，則一般選取與高度無關(guān)的脈動(dòng)風(fēng)速譜，如Davenport譜。

(3) 在傳統(tǒng)的POD法中，往往需要成千上萬個(gè)隨機(jī)變量來實(shí)現(xiàn)對(duì)脈動(dòng)風(fēng)速時(shí)-空隨機(jī)場(chǎng)的模擬，極大地增加了計(jì)算工作量。在本文方法中，隨機(jī)函數(shù)將POD模擬公式中的正交隨機(jī)變量集表達(dá)為基本隨機(jī)變量的正交函數(shù)形式，實(shí)現(xiàn)了僅用兩個(gè)基本隨機(jī)變量即可表達(dá)脈動(dòng)風(fēng)速隨機(jī)場(chǎng)。同時(shí)，生成的脈動(dòng)風(fēng)速代表性時(shí)程可構(gòu)成一個(gè)完備的概率集，這為應(yīng)用概率密度演化理論進(jìn)行結(jié)構(gòu)風(fēng)振響應(yīng)和抗風(fēng)可靠度分析奠定了基礎(chǔ)。

[1] SHINOZUKA M, JAN C B, SEYA H. Stochastic methods in wind engineering[J]. Journal of Wind Engineering and Industrial Aerodynamics, 1990, 36: 829-843.

[2] DEODATIS G. Simulation of ergodic multivariate stochastic processes[J]. Journal of Engineering Mechanics, 1996, 122(8): 778-787.

[3] CAO Y H, XIANG H F, ZHOU Y. Simulation of stochastic wind velocity field on long-span bridges[J]. Journal of Engineering Mechanics, 2000, 126(1): 1-6.

[4] CHEN X, KAREEM A. Proper orthogonal decomposition-based modeling, analysis, and simulation of dynamic wind load effects on structures[J]. Journal of Engineering Mechanics, 2005, 131(4): 325-339.

[5] 李杰, 劉章軍. 隨機(jī)脈動(dòng)風(fēng)場(chǎng)的正交展開方法[J]. 土木工程學(xué)報(bào), 2008, 41(2): 49-53.

LI Jie, LIU Zhangjun. Orthogonal expansion method of random fields of wind velocity fluctuations[J]. China Civil Engineering Journal, 2008, 41(2): 49-53.

[6] LIU Zhangjun, CHEN Jianbing, LI Jie. Orthogonal expansion of Gaussian wind velocity field and PDEM-based vibration analysis of wind-excited structures[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2011, 99(12): 1207-1220.

[7] CHEN L, LETCHFORD C W. Simulation of multivariate stationary Gaussian stochastic processes: Hybrid spectral representation and proper orthogonal decomposition approach[J]. Journal of Engineering Mechanics, 2005, 131(8): 801-808.

[8] CARASSALE L, SOLARI G. Wind modes for structural dynamics: a continuous approach[J]. Probabilistic Engineering Mechanics, 2002, 17: 157-166.

[9] LI Jie, PENG Yongbo, YAN Qi. Modeling and simulation of fluctuating wind speeds using evolutionary phase spectrum[J]. Probabilistic Engineering Mechanics, 2013, 32: 48-55.

[10] 劉章軍, 萬勇, 曾波. 脈動(dòng)風(fēng)速過程模擬的正交展開-隨機(jī)函數(shù)方法[J]. 振動(dòng)與沖擊, 2014, 33(8): 120-124.

LIU Zhangjun, WAN Yong, ZENG Bo. Simulation of fluctuating wind processes with an orthogonal expansion-random function approach[J]. Journal of Vibration and Shock, 2014, 33(8): 120-124.

[11] LI Jie, CHEN Jianbing. Stochastic dynamics of Structures[M]. Singapore: John Wiley & Sons, 2009.

[12] 劉章軍, 陳建兵. 結(jié)構(gòu)動(dòng)力學(xué)[M]. 北京: 中國水利水電出版社, 2012.

[13] LIU Zhangjun, LIU Wei, PENG Yongbo. Random function based spectral representation of stationary and non-stationary stochastic processes[J]. Probabilistic Engineering Mechanics, 2016, 45: 115-126.

[14] PAOLA M D. Digital simulation of wind field velocity[J]. Journal of Wind Engineering and Industrial Aerodynamics, 1998, 74/75/76(1): 91-109.

[15] LI Jie, CHEN Jianbing. The number theoretical method in response analysis of nonlinear stochastic structures[J]. Computational Mechanics, 2007, 39(6): 693-708.