徐軍 李哲帆 張洋 李揚(yáng)

(湖南大學(xué) 土木工程學(xué)院，湖南 長沙 410082)

現(xiàn)有結(jié)構(gòu)設(shè)計(jì)規(guī)范通常忽視風(fēng)荷載的動(dòng)力特征，而將其等效為靜力荷載。對(duì)于高層建筑、大跨網(wǎng)殼結(jié)構(gòu)、大跨度橋梁結(jié)構(gòu)等風(fēng)敏感性較強(qiáng)的結(jié)構(gòu)，風(fēng)荷載的動(dòng)力效應(yīng)往往不可忽視。風(fēng)荷載下結(jié)構(gòu)動(dòng)力響應(yīng)最主要來源于脈動(dòng)風(fēng)速場的作用[1]，而脈動(dòng)風(fēng)速場具有顯著的時(shí)空分布特征和隨機(jī)性，通常被視為在不同空間位置相關(guān)的平穩(wěn)隨機(jī)向量過程，或稱為時(shí)-空隨機(jī)場[2]。合理地描述和模擬隨機(jī)脈動(dòng)風(fēng)速場的時(shí)域動(dòng)力特征，對(duì)于高層建筑和大跨結(jié)構(gòu)等復(fù)雜工程結(jié)構(gòu)的抗風(fēng)分析和保障其風(fēng)災(zāi)安全性具有重要意義。

目前，模擬平穩(wěn)隨機(jī)向量過程的主要方法有：譜表達(dá)法[3]、線性濾波法、本征正交分解法[4]以及小波變換方法[5]等。其中，由于譜表達(dá)法僅涉及一系列諧波疊加，其算法相對(duì)簡單且具有嚴(yán)謹(jǐn)?shù)睦碚摶A(chǔ)，在工程實(shí)際中被廣泛采用。譜表達(dá)法的起源可追溯到20世紀(jì)40年代對(duì)一維單變量白噪聲的描述[6]。此后，Shinozuka[7]將其引入模擬一維多變量平穩(wěn)隨機(jī)向量過程；Yang[8]通過引入FFT算法進(jìn)一步提高了平穩(wěn)隨機(jī)向量過程的模擬效率；Shinozuka等[9]從高斯性、無偏性和遍歷性的角度進(jìn)一步詳細(xì)闡述譜表達(dá)法的原理[9]；為生成遍歷性更好的樣本函數(shù)，Deodatis等[10]引入了雙索引頻率，使得譜表達(dá)法成為模擬隨機(jī)向量過程的經(jīng)典方法，被廣泛地應(yīng)用到風(fēng)場、地震、波浪荷載的模擬中。

平穩(wěn)隨機(jī)脈動(dòng)風(fēng)速場模擬的譜表達(dá)法中往往涉及數(shù)目龐大的隨機(jī)變量，通常采用蒙特卡洛模擬的方式產(chǎn)生基本隨機(jī)變量樣本點(diǎn)來得到多個(gè)空間點(diǎn)上的脈動(dòng)風(fēng)速過程。一方面，為精確地捕捉平穩(wěn)隨機(jī)脈動(dòng)風(fēng)速場的統(tǒng)計(jì)特性，所需的蒙特卡洛模擬樣本數(shù)目往往較大，造成后續(xù)結(jié)構(gòu)隨機(jī)風(fēng)振響應(yīng)和可靠度分析的計(jì)算量難以令人接受；另一方面，蒙特卡洛模擬結(jié)果本質(zhì)上并非完備的概率集合，因此不能在概率密度層次精細(xì)地描述風(fēng)場的概率特性[11]。為克服此問題，Chen等[12]建議用一類隨機(jī)諧波函數(shù)的譜表達(dá)法對(duì)隨機(jī)風(fēng)場進(jìn)行模擬，僅需要少量的隨機(jī)變量便可較好地模擬平穩(wěn)隨機(jī)脈動(dòng)風(fēng)速場特性；李杰等[13]從物理隨機(jī)系統(tǒng)角度出發(fā)，提出了隨機(jī)脈動(dòng)風(fēng)速場的物理模型，將若干關(guān)鍵物理參數(shù)描述為隨機(jī)變量；劉章軍等[14]提出基于隨機(jī)函數(shù)的降維模擬方法，結(jié)合譜表達(dá)方法僅用1-2個(gè)基本隨機(jī)變量即可在二階統(tǒng)計(jì)意義上對(duì)平穩(wěn)隨機(jī)脈動(dòng)風(fēng)速場進(jìn)行精確描述。在上述研究基礎(chǔ)上，本文擬提出一類高維概率空間選點(diǎn)法，實(shí)現(xiàn)基于譜表達(dá)法的平穩(wěn)脈動(dòng)風(fēng)速場精確與高效模擬，為工程結(jié)構(gòu)隨機(jī)風(fēng)振響應(yīng)與可靠度分析奠定基礎(chǔ)。

1 平穩(wěn)隨機(jī)脈動(dòng)風(fēng)速場的譜表達(dá)

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

(1)

式中，n為空間點(diǎn)數(shù)目，ω為圓頻率，Sii(ω)為自譜密度函數(shù)，Sij(ω)為空間點(diǎn)i處的風(fēng)速分量和空間點(diǎn)j處的風(fēng)速分量互譜密度函數(shù)。一般地，SX(ω)為正定Hermitian矩陣，滿足以下關(guān)系：

Sii(ω)=Sii(-ω)，i=1,2,…,n

(2)

(3)

(4)

為避免直接對(duì)譜密度函數(shù)矩陣進(jìn)行喬列斯基分解而造成存儲(chǔ)困難,文獻(xiàn)[15]中引入如下空間相干函數(shù)矩陣：

(5)

β(ω)=B(ω)B*T(ω)

(6)

那么，互譜密度函數(shù)Sij(ω)可進(jìn)一步表示為

(7)

根據(jù)式(2)-(4)及式(7),功率譜密度函數(shù)矩SX(ω)可進(jìn)一步分解為[16]：

SX(ω)=K(ω)B(ω)B*T(ω)K*T(ω)

(8)

式中,K(ω)為如下對(duì)角矩陣,

(9)

由Fourier-Stieltjes積分可知：

(10)

式中，Z(ω)為增量正交零均值復(fù)向量過程，且滿足

E[dZ(ω)]=0n×1,dZ(-ω)=dZ*(ω)

(11)

E[dZ(ω)dZ*T(ω′)]=SX(ω)dωdω′δ(ω-ω′)

(12)

式中，E表示期望算子，“*”是共軛符號(hào)，δ為狄拉克函數(shù)。

那么,第i個(gè)復(fù)過程增量dZi(ω),i=1,2,…,n,表示為

(13)

式中，Δω=ωu/N為頻率間隔，ωu為頻率上限，N為頻率的截?cái)囗?xiàng)數(shù)；Pgh=Ugh-iVgh表示零均值標(biāo)準(zhǔn)化的復(fù)隨機(jī)變量，Ugh和Vgh為實(shí)正交隨機(jī)變量，且滿足

E[Ugh]=E[Vgh]=0,E[UghVrs]=0;

E[UghUrs]=E[VghVrs]=(δgrδhs)/2;

g,r=1,2,…,n;h,s=1,2,…,N

(14)

為各態(tài)歷經(jīng)性，選用了如下雙索引頻率：

(15)

令Big(ωgh)=ξig(ωgh)-iΨig(ωgh)，由式(15)可知-ωgh=ωg(-h)，式(10)可進(jìn)一步表示為

{ξig(ωgh)(Ughcosωght+Vghsinωght)+

Ψig(ωgh)(Ughsinωght-Vghcosωght)}

(16)

(i=1,2,…,n)

空間相干矩陣β(ω)一般是實(shí)對(duì)稱矩陣，其下三角矩陣β(ω)亦為實(shí)值矩陣，式(16)可進(jìn)一步簡化為

[Ughcosωght+Vghsinωght]

(17)

由此可見：平穩(wěn)脈動(dòng)風(fēng)速隨機(jī)向量過程的模擬即可轉(zhuǎn)化為一系列正交隨機(jī)變量調(diào)制的諧波函數(shù)疊加，進(jìn)而可獲得多個(gè)空間點(diǎn)處的脈動(dòng)風(fēng)速時(shí)程序列。一般來說，式(16)和(17)中涉及的正交隨機(jī)變量(即Ugh和Vgh)一般取為標(biāo)準(zhǔn)正態(tài)分布隨機(jī)變量。由此可見：平穩(wěn)脈動(dòng)風(fēng)速隨機(jī)向量過程的譜表達(dá)中涉及成千上萬個(gè)標(biāo)準(zhǔn)正態(tài)分布隨機(jī)變量(2×n×N個(gè))，為典型高維隨機(jī)問題。傳統(tǒng)的蒙特卡洛模擬所需樣本數(shù)量往往過大，模擬效率低下，亦會(huì)對(duì)后續(xù)結(jié)構(gòu)隨機(jī)動(dòng)力響應(yīng)分析效率造成嚴(yán)重負(fù)擔(dān)。為結(jié)合概率密度演化方法進(jìn)行結(jié)構(gòu)隨機(jī)風(fēng)振響應(yīng)與可靠度分析，需產(chǎn)生脈動(dòng)風(fēng)速場典型代表樣本及其賦得概率，有必要對(duì)高維正態(tài)空間離散代表點(diǎn)的選取問題展開專門研究，以期產(chǎn)生少量的代表點(diǎn)數(shù)目，構(gòu)成一個(gè)完備的概率集，提高后續(xù)結(jié)構(gòu)隨機(jī)動(dòng)力響應(yīng)的計(jì)算效率。

2 高維正態(tài)空間離散代表點(diǎn)選取法

借鑒部分分層抽樣的思想[17]，本文提出一類新的高維正態(tài)空間離散代表點(diǎn)選取方法來解決這一問題，由如下3步構(gòu)成：

步驟1 首先，將高維標(biāo)準(zhǔn)正態(tài)空間剖分為一系列2維正交子空間；

步驟2 其次，在每個(gè)正交子空間中應(yīng)用現(xiàn)有低維概率空間選點(diǎn)方法進(jìn)行子空間離散代表點(diǎn)選取和賦得概率的計(jì)算；

步驟3 最后，將子空間中的離散代表點(diǎn)和賦得概率通過隨機(jī)映射的方式獲得原始高維標(biāo)準(zhǔn)正態(tài)空間的離散代表點(diǎn)及其賦得概率。將上述高維標(biāo)準(zhǔn)正態(tài)空間離散代表點(diǎn)及其賦得概率代入至式(17)中，則可生成代表性風(fēng)速過程樣本進(jìn)行有限元計(jì)算。

2.1 高維概率空間分解

由式(17)可見：平穩(wěn)隨機(jī)脈動(dòng)風(fēng)速場譜表達(dá)模擬中共包括2Q(Q=n×N)個(gè)獨(dú)立標(biāo)準(zhǔn)正態(tài)隨機(jī)變量(即Ugh和Vgh)，涉及高維概率空間問題。鑒于直接處理高維隨機(jī)問題的困難，本文首先將原始高維正態(tài)空間分解為一系列正交的低維子空間。根據(jù)文獻(xiàn)[18]中對(duì)正交子空間的定義，可以將高維正態(tài)空間Ωr劃分為Q=n×N個(gè)二維正交子空間，即：

(18)

式中，Ωi,2表示第i個(gè)2維標(biāo)準(zhǔn)正態(tài)子空間，例如，1個(gè)10維標(biāo)準(zhǔn)空間可以剖分為5個(gè)2維正交子空間。

2.2 正交子空間中離散代表點(diǎn)的選取

在獲得一系列2維正交子空間后，對(duì)正交子空間采用數(shù)論選點(diǎn)方法確定其離散代表點(diǎn)及賦得概率。首先，采用華羅庚-王元的數(shù)論方法獲得二維子空間上的均勻分布點(diǎn)集uk,k=1,2,…，nsel[19]:

(19)

(k=1,…，F(xiàn)m)

式中，uk=[u1,k,u2,k],Fm=Fm-1+Fm-2,m=2,3,…;F0=F1=1為斐波那契數(shù)列，int(·)表示取數(shù)值的整數(shù)。本文取m=13,點(diǎn)集總數(shù)目nsel=F13=377，二維均勻分布點(diǎn)集如圖1所示。

圖1 二維均勻分布點(diǎn)集Fig.1 Uniform points in the 2-dimensional unit square

(20)

(j=1,2;k=1,2,…,nsel)

式中，Φ-1指的是標(biāo)準(zhǔn)正態(tài)分布的逆累積分布函數(shù)。

圖2所示為圖1中點(diǎn)集轉(zhuǎn)化至標(biāo)準(zhǔn)正態(tài)空間后所得到的子空間初始離散代表點(diǎn)結(jié)果。

圖2 二維標(biāo)準(zhǔn)正態(tài)空間點(diǎn)集Fig.2 Points in the 2-dimensional standard normal space

(21)

(22)

圖3所示為子空間中離散代表點(diǎn)所對(duì)應(yīng)的Voronoi域剖分，計(jì)算所得的子空間離散代表點(diǎn)賦得概率如圖4所示。

圖3 子空間中Voronoi域剖分Fig.3 Voronoi dissection in subspaces

圖4 子空間離散代表點(diǎn)的賦得概率Fig.4 Assigned probabilities for points set in subspaces

(23)

(j=1,2)

若A成立，I{A}=1；否則I{A}=0。圖5所示為初始離散代表點(diǎn)和重整后離散代表點(diǎn)的對(duì)比圖，可見：重整后點(diǎn)集的位置與初始點(diǎn)集位置略有不同。

圖5 初始點(diǎn)集與重整后點(diǎn)集對(duì)比Fig.5 Comparison between the initial points set and the one after rearrangement

2.3 高維離散代表點(diǎn)及賦得概率

(24)

圖6 點(diǎn)集隨機(jī)配對(duì)Fig.6 Random permutation of the point set

由于低維子空間均正交，可認(rèn)為高維離散代表點(diǎn)的賦得概率為上述隨機(jī)配對(duì)的低維正交子空間離散代表點(diǎn)賦得概率的乘積：

(25)

進(jìn)一步,由歸一化條件得到高維離散代表點(diǎn)賦得概率為:

(26)

3 框架-剪力墻結(jié)構(gòu)平穩(wěn)隨機(jī)脈動(dòng)風(fēng)場模擬

文中采用某三十層框架-剪力墻結(jié)構(gòu)的平穩(wěn)隨機(jī)脈動(dòng)風(fēng)場來驗(yàn)證本文提出的方法的正確性。該框架-剪力墻結(jié)構(gòu)立面圖和平面圖如圖7所示(黑色區(qū)域代表剪力墻和框架柱)。模型底部柱高為6 m，標(biāo)準(zhǔn)層高度為3 m，迎風(fēng)寬度15 m，坐標(biāo)系xoz與風(fēng)向平行。將該結(jié)構(gòu)上的隨機(jī)脈動(dòng)風(fēng)場近似為45個(gè)空間點(diǎn)處的隨機(jī)脈動(dòng)風(fēng)速矢量過程。其空間分布如圖7所示。表1所示為該模型中用于ANSYS建模的單元、尺寸參數(shù)和材料參數(shù)信息。

表1 模型相關(guān)參數(shù)Table 1 Related parameters of the model

圖7 結(jié)構(gòu)模型立面圖與平面圖Fig.7 Elevation and plan layout of structural model

本例中，采用Davenport譜作為雙邊自功率譜密度函數(shù)，定義如下[21]：

(27)

(i=1,2,…,n)

采用本文所建議的高維標(biāo)準(zhǔn)正態(tài)空間選點(diǎn)方法，可以模擬結(jié)構(gòu)上45個(gè)空間點(diǎn)(n=45)的脈動(dòng)風(fēng)速矢量代表過程。表2所示為模擬隨機(jī)脈動(dòng)風(fēng)速矢量過程所需的參數(shù)。

表2 隨機(jī)脈動(dòng)風(fēng)速場模擬所需相關(guān)參數(shù)Table 2 Related parameters required to simulate a random fluctuating wind speed vector

選擇高度為6、42、78 m的空間點(diǎn)1、19和37的仿真結(jié)果來說明本文方法的正確性和效率。圖8所示為3個(gè)空間點(diǎn)的平穩(wěn)風(fēng)速時(shí)程代表性樣本過程，可以清楚地看到：3個(gè)風(fēng)速時(shí)程過程樣本的值域范圍大致相同，且樣本在該值域范圍內(nèi)波動(dòng)，呈現(xiàn)平穩(wěn)特性。

圖8 風(fēng)速時(shí)程曲線代表性樣本Fig.8 Representative samples of fluctuating wind speeds

由國內(nèi)外相關(guān)文獻(xiàn)可知:一般選擇均值誤差、標(biāo)準(zhǔn)差誤差、功率譜密度函數(shù)平均值誤差3項(xiàng)指標(biāo)來考察平穩(wěn)風(fēng)速過程的統(tǒng)計(jì)特性,從而進(jìn)一步驗(yàn)證本文方法的準(zhǔn)確性。那么，第q點(diǎn)模擬樣本過程的均值、標(biāo)準(zhǔn)差和功率譜密度函數(shù)平均值可表示為

(28)

(29)

(30)

(31)

(32)

圖9 均值模擬值與目標(biāo)值的比較Fig.9 Comparison between simulated mean values and target values

圖10 標(biāo)準(zhǔn)差模擬值與目標(biāo)值比較Fig.10 Comparison between simulated standard deviation values and target values

圖11 功率譜密度函數(shù)平均值與目標(biāo)值比較Fig.11 Comparison between the mean of simulated PSD values and target values

從圖中不難發(fā)現(xiàn)：采用本文方法模擬的平穩(wěn)隨機(jī)脈動(dòng)風(fēng)場的均值、標(biāo)準(zhǔn)差及功率譜密度函數(shù)平均值與目標(biāo)值均吻合良好，表明本文方法具有較高的精度。此外，由于本文方法僅采用帶有賦得概率的377個(gè)高維標(biāo)準(zhǔn)正態(tài)空間離散代表點(diǎn)來模擬平穩(wěn)隨機(jī)脈動(dòng)風(fēng)場，亦表明本文方法具有較高的計(jì)算效率。進(jìn)一步地，采用平均相對(duì)誤差進(jìn)行定量分析來表明本文方法的準(zhǔn)確性。那么，第q個(gè)空間點(diǎn)脈動(dòng)風(fēng)速隨機(jī)過程均值、標(biāo)準(zhǔn)差和功率譜密度平均值的平均相對(duì)誤差可定義為[22]：

(33)

(34)

(35)

式中，Nt表示時(shí)間步數(shù)，Nω表示頻率步數(shù)。同時(shí)，采用相同樣本數(shù)目的蒙特卡洛模擬和文獻(xiàn)[3]中的降維譜表達(dá)法，與本文方法進(jìn)行對(duì)比，結(jié)果見表3。

表3 3種方法平均相對(duì)誤差比較Table 3 Comparison of average relative errors among three methods %

從表中可以發(fā)現(xiàn)：相較于這兩種方法，本文方法計(jì)算得到的均值、標(biāo)準(zhǔn)差及功率譜密度平均值的精度均較高，相對(duì)誤差更小。特別地，均值的精度可以提高1-2個(gè)量級(jí)，這意味著：可利用較少數(shù)目樣本模擬隨機(jī)脈動(dòng)風(fēng)場，滿足精度要求，從而進(jìn)一步降低后續(xù)結(jié)構(gòu)隨機(jī)風(fēng)振響應(yīng)與可靠度分析的計(jì)算量。

將上述平穩(wěn)隨機(jī)脈動(dòng)風(fēng)場導(dǎo)入至通用有限元軟件(ANSYS)，對(duì)結(jié)構(gòu)模型進(jìn)行確定性風(fēng)振響應(yīng)分析，可以獲得每個(gè)風(fēng)速矢量樣本相對(duì)應(yīng)的結(jié)構(gòu)動(dòng)力響應(yīng)?？疾旖Y(jié)構(gòu)的層間位移角

(36)

式中，xl(t)表示在第l層在x方向上的水平位移，hl表示第l層的高度。那么，可方便地計(jì)算得到層間位移角的均值與標(biāo)準(zhǔn)差。圖12所示為頂層層間位移角的均值與標(biāo)準(zhǔn)差。進(jìn)一步地，將這些帶有概率信息的響應(yīng)過程代入至概率密度演化方法中，即可方便地獲得結(jié)構(gòu)隨機(jī)風(fēng)振響應(yīng)概率密度演化全過程及動(dòng)力可靠度信息。

圖12 頂層層間位移角響應(yīng)的均值和標(biāo)準(zhǔn)差Fig.12 Mean and standard deviation of the top inter-storey drift

考慮到結(jié)構(gòu)形式的多樣化，筆者還選取了K8單層球面網(wǎng)殼進(jìn)一步驗(yàn)證了本文方法的有效性，限于篇幅，本文不再贅述。

4 結(jié)語

采用譜表達(dá)方法模擬平穩(wěn)隨機(jī)脈動(dòng)風(fēng)場隨機(jī)向量過程往往涉及高維標(biāo)準(zhǔn)正態(tài)隨機(jī)變量，對(duì)數(shù)值模擬的精度和效率提出挑戰(zhàn)。鑒于此，本文借鑒部分分層抽樣思想，提出了一種新的高維正態(tài)空間選點(diǎn)方法,這一方法采用了Voronoi域剖分、點(diǎn)集重整和隨機(jī)配對(duì)等技術(shù)手段獲得高維正態(tài)空間離散代表點(diǎn)及其賦得概率，將其代入至平穩(wěn)隨機(jī)脈動(dòng)風(fēng)場譜表達(dá)模型中，即可得到帶有完備概率信息的風(fēng)場代表性樣本過程。

采用在30層框架-剪力墻結(jié)構(gòu)的平穩(wěn)隨機(jī)風(fēng)速隨機(jī)向量過程模擬驗(yàn)證了本文方法的準(zhǔn)確性與高效性，研究結(jié)果表明：本文方法能夠準(zhǔn)確地獲得平穩(wěn)脈動(dòng)風(fēng)速隨機(jī)向量過程的統(tǒng)計(jì)信息。值得注意的是：相較于蒙特卡洛模擬方法，本文方法能用較少的樣本數(shù)構(gòu)成完備的概率集，可大幅減少后續(xù)的結(jié)構(gòu)隨機(jī)風(fēng)振響應(yīng)與可靠度分析所需的有限元分析計(jì)算量，展現(xiàn)了這一方法在工程實(shí)際運(yùn)用中的優(yōu)越性。

若考慮隨機(jī)脈動(dòng)風(fēng)場的非高斯性與非平穩(wěn)性，擬進(jìn)一步拓展這一方法的適用性。