隨機風場模擬的連續(xù)本征正交分解-隨機函數(shù)方法

劉章軍， 劉增輝

(1. 防災減災湖北省重點實驗室(三峽大學)， 湖北 宜昌 443002； 2. 三峽大學 土木與建筑學院， 湖北 宜昌 443002)

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

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

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

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

(1)

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

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

(2a)

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

(2b)

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

(3)

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

(4)

λi(ω)δij

(5)

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

(6)

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

(7)

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

(8)

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

(9)

(10)

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

f0(x,t)=

(11)

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

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

f(x,t)=

(12)

式中：f(x,t)為模擬的時-空隨機場；ωm=mΔω，ωu=NΔω為截斷圓頻率；N為頻率截斷項數(shù)；Pkm=Wk(ωm)為一組零均值的正交復隨機變量，滿足如下的正交性

(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為零均值的實正交隨機變量，滿足如下的基本條件

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

(16)

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

Rkmcosωmt+Ikmsinωmt

(17)

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

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

(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分別為基本隨機變量Θ1和Θ2的定義區(qū)間，pΘ1(θ1)和pΘ2(θ2)分別為基本隨機變量Θ1和Θ2的概率密度函數(shù)。

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

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

(19)

3 脈動風場的互功率譜密度函數(shù)及其特征問題的封閉解

一般地，假定脈動風場是一個零均值的非完全均勻的連續(xù)隨機場v(x,t)，定義在0≤x≤L和0

(20)

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

x,x′∈[0,L]

(21)

結合式(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,ω)關于空間域x∈[0,L]的中點x=L/2具有對稱性，即：當tan(μk/2)=α/μk時，特征函數(shù)關于中點x=L/2是正對稱的；當tan(μk/2)=-μk/α時，特征函數(shù)關于中點x=L/2是反對稱的，這與結構動力學中的結構振型具有類似的物理意義。事實上，特征函數(shù)表征了隨機風場關于空間域的各階模態(tài)形狀[14]。

4 數(shù)值算例

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

(27)

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

(28)

式中：z0為地面粗糙長度，本文取z0=0.03 m。

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

(29)

為簡便之，僅以水平坐標x=50 m，80 m和120 m三點處的脈動風速進行分析。

在隨機風場的數(shù)值模擬中，取截斷頻率ωu=2π rad/s，頻率步長Δω=0.01 rad/s，則頻率截斷項數(shù)N=628；持時T=600 s，時間步長Δt=0.1 s。同時，根據(jù)式(7)可確定POD模擬公式(17)中的展開項數(shù)n，為保證計算精度，本文要求展開精度Π(n)≥90%。

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

(a) 前五階特征值

(b) 前20階特征值對應的展開精度

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

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

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

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

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

圖4分別給出了x=50 m與x=80 m處、x=50 m與x=120 m處以及x=80 m與x=120 m處610條脈動風速代表性時程的集合互功率譜與目標互功率譜的比較。從圖4可知，模擬的互功率譜與目標互功率譜擬合較好，進一步證明本方法的有效性。

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

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

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

5 結 論

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

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

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

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

[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]. 土木工程學報, 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] 劉章軍, 萬勇, 曾波. 脈動風速過程模擬的正交展開-隨機函數(shù)方法[J]. 振動與沖擊, 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] 劉章軍, 陳建兵. 結構動力學[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.