橋梁顫振臨界風(fēng)速的概率密度演化計(jì)算

姜保宋, 周志勇, 唐 峰

(土木工程防災(zāi)國(guó)家重點(diǎn)實(shí)驗(yàn)室(同濟(jì)大學(xué))， 上海 200092)

顫振臨界風(fēng)速是衡量橋梁顫振穩(wěn)定性的重要指標(biāo)[1]. 然而，由于橋梁結(jié)構(gòu)自身及外部環(huán)境的隨機(jī)性，如剛度、質(zhì)量、阻尼比、氣動(dòng)導(dǎo)數(shù)等的不確定，從而造成橋梁的顫振臨界風(fēng)速也是不確定的[2]. 所以橋梁顫振失效概率采用可靠度分析方法更為合理[3].

目前在橋梁顫振臨界風(fēng)速概率評(píng)價(jià)方面，主要還是采用矩方法[4]及替代函數(shù)法[5]進(jìn)行估計(jì). 基于現(xiàn)有的結(jié)構(gòu)可靠度理論，一些學(xué)者提出橋梁顫振穩(wěn)定失效概率的計(jì)算方法：2001年，葛耀君教授運(yùn)用中心點(diǎn)法和驗(yàn)算點(diǎn)法，基于實(shí)驗(yàn)與理論相結(jié)合進(jìn)行橋梁顫振穩(wěn)定的概率分析，并對(duì)上海楊浦大橋和江陰長(zhǎng)江大橋進(jìn)行了顫振穩(wěn)定的概率性評(píng)價(jià)[6]. 2007年，周崢等提出的隨機(jī)有限元方法，綜合考慮了結(jié)構(gòu)剛度、質(zhì)量、阻尼比以及顫振導(dǎo)數(shù)等隨機(jī)因素對(duì)顫振臨界風(fēng)速的影響，并對(duì)東海大橋顆珠山斜拉橋進(jìn)行了顫振可靠性分析[7].

目前，對(duì)于線性系統(tǒng)平穩(wěn)反應(yīng)的分析理論基本趨于完備，并已經(jīng)開(kāi)始進(jìn)入工程實(shí)用階段[8]. 然而，對(duì)于非線性系統(tǒng)，人們?nèi)匀幻媾R巨大的困難[9]. 近幾年，李杰和陳建兵提出的廣義概率密度演化方法[10]，為求解非線性系統(tǒng)帶來(lái)了巨大的方便. 應(yīng)用數(shù)值求解技術(shù)可獲得廣義概率密度演化方程的解，從而可以得到結(jié)構(gòu)響應(yīng)的概率密度函數(shù)(PDF)曲線，進(jìn)而獲得結(jié)構(gòu)的可靠度. 概率密度演化方法一方面包含結(jié)構(gòu)中所有的概率信息，另一方面不受概率密度分布情況的影響，對(duì)于復(fù)雜形式的狀態(tài)函數(shù)也適用，為結(jié)構(gòu)可靠度分析提供了極大的方便.

本文采用概率密度演化方法與橋梁顫振多模態(tài)耦合分析方法相結(jié)合，以江陰長(zhǎng)江大橋?yàn)槔x取結(jié)構(gòu)質(zhì)量、結(jié)構(gòu)剛度、結(jié)構(gòu)阻尼比、氣動(dòng)導(dǎo)數(shù)作為隨機(jī)變量，對(duì)結(jié)構(gòu)模態(tài)阻尼比及頻率的概率密度演化過(guò)程進(jìn)行分析，并給出基于概率密度演化方法得到的顫振臨界風(fēng)速. 論文首先概述了目前結(jié)構(gòu)可靠度研究的進(jìn)展以及目前橋梁顫振臨界風(fēng)速可靠度計(jì)算的方法，第1節(jié)介紹概率密度演化方法的基本理論以及對(duì)橋梁進(jìn)行顫振臨界風(fēng)速概率密度演化分析的主要步驟，第2節(jié)對(duì)典型橋梁的基本參數(shù)進(jìn)行了介紹，第3節(jié)為概率演化計(jì)算結(jié)果及分析，第4節(jié)對(duì)研究結(jié)果進(jìn)行了總結(jié).

1 基本理論

1.1 廣義概率密度演化方程

一般的結(jié)構(gòu)動(dòng)力系統(tǒng)運(yùn)動(dòng)方程可以表述為

(1)

結(jié)構(gòu)的基本參數(shù)變異性較大[11-12]. 可以將這些參數(shù)表示為一系列基本隨機(jī)變量的函數(shù). 記為θ1=(θ1，θ2，…,θs1)，其中s1為結(jié)構(gòu)隨機(jī)變量的總個(gè)數(shù)，這樣，方程中的質(zhì)量矩陣、剛度矩陣和阻尼比矩陣都是隨機(jī)向量θ1的函數(shù).

同時(shí)，結(jié)構(gòu)外部激勵(lì)也具有顯著的隨機(jī)性，可將上述隨機(jī)過(guò)程表述為θ2=(θs1+1，θs1+2，…，θs1+s2)，其中s2為隨機(jī)激勵(lì)中基本隨機(jī)變量的個(gè)數(shù).

上述全部隨機(jī)變量記為θ=(θ1,θ2)=(θ1,θ2，…，θs)，其中s=s1+s2為隨機(jī)變量的總個(gè)數(shù). 其中θ的聯(lián)合概率密度函數(shù)pθ(θ)已知，其分布域?yàn)棣甫龋谑牵匠?1)可表示為

(2)

記Z=Z(Z1,Z2,…,Zm)為所要考察的物理量，令

(3)

其中φ(.)是從狀態(tài)向量向所考察的物理量的轉(zhuǎn)換函數(shù)，對(duì)應(yīng)的概率密度函數(shù)Pzθ(z,θ,t) 滿足方程：

(4)

一般的，方程(4)的邊界條件可采用：

Pzθ(z,θ,t)|zj→±=0,j=1,2,…,m.

(5)

而初始條件則為

Pzθ(z,θ,t)|t=t0=δ(z-z0)Pθ(θ).

(6)

z0為初始值，通過(guò)求解廣義概率密度演化方程(4)，可最終得到Z(t)的概率密度函數(shù)：

(7)

對(duì)于一般的隨機(jī)動(dòng)力系統(tǒng)，廣義概率密度演化方程的求解需要借助于有限差分法：一種是單邊差分格式，另外一種是Lax-Wendroff格式. 李杰和陳建兵采用具有TVD性質(zhì)的有限差分格式可以獲得較為理想解答[13]. 最近，Papadopoulos[14]和Kalogeris采用Petrov-Galerkin有限元方法也成功的應(yīng)用于概率密度演化方程的求解.

1.2 廣義概率密度演化方程的求解

(8)

圖1 廣義概率密度演化方程求解流程

Fig.1 Flow chart for the solution of the generalized probability density evolution equation (GDEE)

1.3 橋梁顫振的概率密度演化計(jì)算

1.3.1 顫振臨界風(fēng)速的計(jì)算

本文基于結(jié)構(gòu)有限元模型，進(jìn)行多模態(tài)顫振耦合分析[15]. 橋梁結(jié)構(gòu)進(jìn)行顫振分析時(shí)，自激力作用下的結(jié)構(gòu)運(yùn)動(dòng)方程一般形式[16]可表示為

(9)

(10)

假設(shè)結(jié)構(gòu)的特征運(yùn)動(dòng)形式可描述為

X=φeλt.

(11)

(12)

1.3.2 概率密度演化計(jì)算

按照式(10)，物理方程可建立為

(13)

令

則有

(14)

至此，通過(guò)求解G的特征值，即可求得各參與模態(tài)在不同風(fēng)速下的特征頻率及阻尼比，方程(2)的物理方程變?yōu)?/p>

eig(G(U,θ)).

(15)

方程中U代表風(fēng)速. 從上述分析可以看出，橋梁顫振臨界風(fēng)速的概率密度演化計(jì)算過(guò)程只需將對(duì)應(yīng)的物理方程換成方程(15)，將風(fēng)速看成廣義的時(shí)間變量，即可求解，求解流程如圖1所示.

2 橋梁基本參數(shù)

2.1 結(jié)構(gòu)參數(shù)

本文選取江陰長(zhǎng)江大橋作為本次計(jì)算的實(shí)例. 大橋主橋跨徑布置為336.5 m+1 385 m+309.34 m. 主梁采用扁平狀閉口鋼箱梁，主梁寬36.9 m，主梁高3.0 m. 主纜橫橋向中心間距為32.5 m，主纜矢跨比為1/10.5，吊桿縱橋向間距為16 m. 橋塔為門(mén)式框架結(jié)構(gòu)，南北橋塔高分別為187、184 m. 橋梁總體立面布置圖及主梁斷面布置圖參見(jiàn)圖2、3.

圖2 橋梁總體立面布置(m)

圖3 橋梁主梁斷面(m)

結(jié)構(gòu)的主要振型、頻率及等效質(zhì)量列于表1，結(jié)構(gòu)一階正對(duì)稱豎彎頻率為0.133 706 Hz，結(jié)構(gòu)一階正對(duì)稱扭轉(zhuǎn)頻率為0.272 963 Hz.

表1 結(jié)構(gòu)主要振型及等效質(zhì)量

2.2 氣動(dòng)導(dǎo)數(shù)

2.3 隨機(jī)參數(shù)選取及取值

表2 結(jié)構(gòu)氣動(dòng)導(dǎo)數(shù)

3 橋梁顫振的概率密度演化分析

3.1 代表點(diǎn)的選取及檢驗(yàn)

代表點(diǎn)選取采用基于概率空間剖分的兩步選點(diǎn)法[10,20]，主要步驟如下：1)獲得一個(gè)在Ωe分布空間均勻散布的點(diǎn)集;2)根據(jù)Ωe概率分布特性對(duì)此點(diǎn)集進(jìn)行變換，獲得F-偏差及一階和二階偏差均較小的點(diǎn)集作為最終的離散代表點(diǎn)集. 本文選取離散代表點(diǎn)2 000個(gè)，對(duì)應(yīng)的各隨機(jī)變量的分布及相關(guān)系數(shù)如圖4所示，各隨機(jī)變量的分布和初始假定分布一致，且各隨機(jī)變量之間的相關(guān)系數(shù)均小于3%，選點(diǎn)符合要求. 圖5、6為初始風(fēng)速下結(jié)構(gòu)頻率及阻尼比與各隨機(jī)變量的分布及相關(guān)性.

圖4 各隨機(jī)變量分布及相互之間的相關(guān)性

圖5 第15階頻率分布及其與其他隨機(jī)變量相關(guān)性

Fig.5 Fifteenth modal frequency distribution and its correlation with other random variables

圖6 第15階阻尼比分布及其與其他隨機(jī)變量相關(guān)性

Fig.6 Fifteenth modal damping ratio distribution and its correlation with other random variables

3.2 概率密度演化過(guò)程

概率演化過(guò)程的計(jì)算參見(jiàn)第1節(jié)的介紹，本文顫振分析采用多模態(tài)耦合顫振分析方法進(jìn)行計(jì)算，選取結(jié)構(gòu)前30階模態(tài)參與顫振分析及概率密度演化分析，給出各個(gè)模態(tài)頻率及阻尼比的概率密度演化過(guò)程. 為方便敘述，本文僅選取具有代表性的模態(tài)對(duì)其計(jì)算結(jié)果進(jìn)行闡述.

3.2.1 阻尼比

圖7為結(jié)構(gòu)一階對(duì)稱豎彎阻尼比概率演化過(guò)程. 7(a)對(duì)應(yīng)的概率密度隨風(fēng)速的演進(jìn)過(guò)程，7(b)對(duì)應(yīng)的是特定風(fēng)速下的概率密度分布. 在低風(fēng)速下，一階對(duì)稱豎彎阻尼比的概率密度分布可近似認(rèn)為服從正態(tài)或?qū)?shù)正態(tài)分布，但方差或者變異系數(shù)會(huì)隨著風(fēng)速的加大而不斷增大，且變異度不斷加劇，而在高風(fēng)速下(v=64 m/s),阻尼比的概率分布已經(jīng)無(wú)法用單一的概率分布來(lái)描述，概率分布會(huì)出現(xiàn)雙峰甚至多峰的情況.

圖8為結(jié)構(gòu)一階對(duì)稱扭轉(zhuǎn)阻尼比概率演化過(guò)程. 8(a)對(duì)應(yīng)的概率密度隨風(fēng)速的演進(jìn)過(guò)程，8(b)對(duì)應(yīng)的是特定風(fēng)速下的概率密度分布. 低風(fēng)速下，一階對(duì)稱扭轉(zhuǎn)阻尼比的概率密度分布可近似認(rèn)為服從正態(tài)或?qū)?shù)正態(tài)分布，方差或者變異系數(shù)會(huì)隨著風(fēng)速的加大而不斷增大，且變異度不斷加劇，其阻尼比均值則呈現(xiàn)先減小后增大的趨勢(shì)，這與按照確定性方法得到阻尼比隨風(fēng)速的變化曲線一致(在3.3節(jié)中論述). 高風(fēng)速下的概率密度則無(wú)法用單一的概率分布來(lái)描述，概率分布會(huì)出現(xiàn)雙峰甚至多峰的情況.

(a)概率密度演化過(guò)程

(b)特定風(fēng)速下的概率分布

Fig.7 Probability density evolution of modal damping ratio for the 1st-S-V mode

(a)概率密度演化過(guò)程

(b)特定風(fēng)速下的概率分布

Fig.8 Probability density evolution of modal damping ratio for the 1st-S-T mode

3.2.2 頻率

圖9為結(jié)構(gòu)一階對(duì)稱豎彎頻率概率演化過(guò)程. 9(a)對(duì)應(yīng)的是概率密度隨風(fēng)速的演進(jìn)過(guò)程，9(b)對(duì)應(yīng)的是特定風(fēng)速下的概率密度分布. 低風(fēng)速下，一階對(duì)稱豎彎頻率的概率密度分布可近似認(rèn)為服從正態(tài)或?qū)?shù)正態(tài)分布，但方差或者變異系數(shù)會(huì)隨著風(fēng)速的加大而增大，但增大率較小. 而在高風(fēng)速下,頻率的概率分布已經(jīng)無(wú)法用單一的概率分布來(lái)描述，概率分布會(huì)出現(xiàn)雙峰甚至多峰的情況. 豎彎頻率均值隨著風(fēng)速加大而加大，與按照確定性方法得到豎彎頻率隨風(fēng)速的變化曲線一致(在3.3中論述).

(a)概率密度演化過(guò)程

(b)特定風(fēng)速下的概率分布

Fig.9 Probability density evolution of modal frequency for the 1st-S-V mode

圖10為結(jié)構(gòu)一階對(duì)稱扭轉(zhuǎn)頻率概率演化過(guò)程. 10(a)對(duì)應(yīng)的是概率密度隨風(fēng)速的演進(jìn)過(guò)程，10(b)對(duì)應(yīng)的是特定風(fēng)速下的概率密度分布. 一階對(duì)稱扭轉(zhuǎn)頻率的概率密度分布可近似認(rèn)為服從正態(tài)或?qū)?shù)正態(tài)分布，但方差或者變異系數(shù)會(huì)隨著風(fēng)速的加大而增大. 頻率均值隨風(fēng)速的增大而減小，這與按照確定性方法得到頻率隨風(fēng)速的變化曲線一致(在3.3中論述).

(a)概率密度演化過(guò)程

(b)特定風(fēng)速下的概率分布

Fig.10 Probability density evolution of modal frequency for the 1st-S-T mode

3.3 顫振臨界風(fēng)速評(píng)價(jià)

3.3.1 按確定性方法的計(jì)算結(jié)果

確定性方法采用多模態(tài)耦合分析方法確定顫振臨界風(fēng)速，選取結(jié)構(gòu)前30階模態(tài)參與顫振臨界風(fēng)速多模態(tài)計(jì)算. 判別顫振臨界風(fēng)速標(biāo)準(zhǔn)為任一參與多模態(tài)計(jì)算的模態(tài)阻尼比大于零，模態(tài)初始阻尼比均取為0.005. 計(jì)算方法參見(jiàn)第2節(jié)中的論述. 結(jié)構(gòu)主要模態(tài)頻率及阻尼比隨風(fēng)速的變化曲線如圖11所示，以阻尼比大于零作為顫振臨界風(fēng)速的判別標(biāo)準(zhǔn)，得到顫振臨界風(fēng)速vcr=68.7 m/s. 與文獻(xiàn)[17]中的vcr=67.69 m/s較為一致，故下文采用vcr=68.7 m/s.

3.3.2 阻尼比的分布

為了采用概率密度演化方法對(duì)顫振臨界風(fēng)速進(jìn)行評(píng)價(jià)，需要首先確定阻尼比在不同風(fēng)速下的概率分布情況，進(jìn)而確定標(biāo)準(zhǔn)差及置信區(qū)間. 圖12～14為結(jié)構(gòu)主要模態(tài)阻尼比不同風(fēng)速下的均值、概率分布及置信線.

(a)頻率

(b)阻尼比

Fig.11 Variation of modal damping ratio and frequency with wind speed

圖12 結(jié)構(gòu)一階對(duì)稱豎彎阻尼比不同風(fēng)速下的均值、概率分布及置信線

Fig.12 Mean, probability density function(PDF), and confidence curves of modal damping ratio for the 1st-S-V mode at different wind speeds

通過(guò)對(duì)比可以發(fā)現(xiàn)，低風(fēng)速下各主要模態(tài)阻尼比在不同風(fēng)速下的概率分布可近似認(rèn)為服從正態(tài)分布或?qū)?shù)正態(tài)分布. 標(biāo)準(zhǔn)差隨風(fēng)速的增大有所增大，近似服從指數(shù)分布，其中一階正對(duì)稱豎彎對(duì)應(yīng)的標(biāo)準(zhǔn)差增幅最為明顯，一階側(cè)彎增幅較小，而扭轉(zhuǎn)模態(tài)對(duì)應(yīng)的阻尼比在低風(fēng)速下增幅較小，在高風(fēng)速下增幅較大(參見(jiàn)圖15).

圖13 結(jié)構(gòu)一階對(duì)稱扭轉(zhuǎn)阻尼比不同風(fēng)速下的均值、概率分布及置信線

Fig.13 Mean, PDFs, and confidence curves of modal damping ratio for the 1st-S-T mode at different wind speeds

圖14 結(jié)構(gòu)一階反對(duì)稱扭轉(zhuǎn)阻尼比不同風(fēng)速下的均值、概率分布及置信線

Fig.14 Mean, PDFs, and confidence curves of modal damping ratio for the 1st-A-T mode at different wind speeds

圖15 結(jié)構(gòu)阻尼比標(biāo)準(zhǔn)差隨風(fēng)速的變化曲線

Fig.15 Variation of standard deviation of modal damping ratio with wind speed

3.3.3 顫振臨界風(fēng)速

基于概率密度演化方法的顫振臨界風(fēng)速確定如下：采用99%置信曲線確定各模態(tài)阻尼比零點(diǎn)對(duì)應(yīng)的風(fēng)速，以最小的風(fēng)速作為顫振臨界風(fēng)速. 根據(jù)3.2節(jié)的研究，本部分僅列出一階對(duì)稱扭轉(zhuǎn)以及一階反對(duì)稱扭轉(zhuǎn)求得的對(duì)應(yīng)風(fēng)速點(diǎn). 如圖16所示，一階正對(duì)稱扭轉(zhuǎn)及一階反對(duì)稱扭轉(zhuǎn)阻尼比概率均值變化曲線與按照確定性方法得到曲線基本重合，說(shuō)明概率密度演化方法可靠性較好. 基于概率密度演化方法，按照一階正對(duì)稱扭轉(zhuǎn)得到的顫振臨界風(fēng)速為vcr=64.82 m/s，按照一階反對(duì)稱扭轉(zhuǎn)得到的顫振臨界風(fēng)速為vcr=72.31 m/s，故最終顫振臨界風(fēng)速為vcr=64.82 m/s.

與按確定性方法到的顫振臨界風(fēng)速相比，顫振臨界風(fēng)速降低了68.7 m/s-64.82 m/s=3.88 m/s. 因此，采用傳統(tǒng)的確定性方法得到的顫振臨界風(fēng)速會(huì)較高地估計(jì)結(jié)構(gòu)的安全度，應(yīng)該采用概率評(píng)價(jià)的方法，進(jìn)行顫振臨界風(fēng)速估算.

(a)一階正對(duì)稱扭轉(zhuǎn)

(b)一階反對(duì)稱扭轉(zhuǎn)

Fig.16 Variation of standard deviation of modal damping ratio and confidence curve with wind speed

4 結(jié) 論

1)將概率密度演化方法與橋梁顫振多模態(tài)耦合分析方法相結(jié)合，以江陰長(zhǎng)江大橋作為實(shí)例, 通過(guò)考慮橋梁自身結(jié)構(gòu)的不確定性以及氣動(dòng)導(dǎo)數(shù)的不確定性，給出不同模態(tài)的阻尼比及頻率的概率密度演化過(guò)程. 按照99%保證率曲線，給出該橋的顫振臨界風(fēng)速，與按確定性方法到的顫振臨界風(fēng)速相比，顫振臨界風(fēng)速降低了3.88 m/s.

2) 結(jié)構(gòu)阻尼比和結(jié)構(gòu)質(zhì)量、結(jié)構(gòu)剛度的相關(guān)系數(shù)為0.2左右，與初始阻尼的相關(guān)系數(shù)為0.91. 而結(jié)構(gòu)阻尼比主要與對(duì)應(yīng)的運(yùn)動(dòng)形態(tài)一致的氣動(dòng)導(dǎo)數(shù)相關(guān)度較大，而與其他氣動(dòng)導(dǎo)數(shù)相關(guān)度較小.

3) 結(jié)構(gòu)模態(tài)阻尼比及頻率在低風(fēng)速情況下的概率密度分布可近似認(rèn)為服從正態(tài)或?qū)?shù)正態(tài)分布，高風(fēng)速下的概率密度則無(wú)法用單一的概率分布來(lái)描述，概率分布會(huì)出現(xiàn)雙峰甚至多峰的情況.

4) 結(jié)構(gòu)模態(tài)阻尼比的標(biāo)準(zhǔn)差隨風(fēng)速的增大有所增大，近似服從指數(shù)增大. 其中一階正對(duì)稱豎彎對(duì)應(yīng)的標(biāo)準(zhǔn)差增幅最為明顯，一階側(cè)彎增幅較小，而扭轉(zhuǎn)模態(tài)對(duì)應(yīng)的阻尼比在低風(fēng)速下增幅較小，在高風(fēng)速下增幅較大.