劉小璐， 蘇 成， 聶 銘

(1. 廣東電網(wǎng)有限責(zé)任公司電力科學(xué)研究院, 廣州 510080； 2. 華南理工大學(xué) 土木與交通學(xué)院, 廣州 510640)

大跨度橋梁不同支承處的地震激勵存在顯著的非一致性，主要表現(xiàn)為行波效應(yīng)、失相干效應(yīng)和局部場地效應(yīng)，分別反映各支承處地震激勵間的相位差、相干性損失和頻譜差異[1]。同時，地震動具有本質(zhì)上的隨機性，應(yīng)該采用隨機振動方法對大跨度橋梁結(jié)構(gòu)開展非一致地震響應(yīng)分析研究。國內(nèi)外學(xué)者已采用不同隨機振動方法開展了大跨度橋梁非一致地震激勵問題研究。例如，Allam和Datta[2]與Dumanogluid和Soyluk[3]采用功率譜法計算了非一致地震激勵下大跨度斜拉橋的平穩(wěn)隨機響應(yīng)；Lin等[4]、焦?？坪屠類廴篬5]與趙雷和劉寧國[6]均采用虛擬激勵法對非一致地震激勵下的大跨度斜拉橋開展隨機振動分析研究；劉小璐等[7]采用時域顯式法研究了非一致地震激勵下大跨度懸索橋的非平穩(wěn)隨機響應(yīng)。

上述研究主要解決非一致地震激勵下的線性隨機振動問題。強震作用下橋梁局部構(gòu)件會進(jìn)入彈塑性狀態(tài)，此時應(yīng)該采用非線性隨機振動方法進(jìn)行抗震分析。國內(nèi)外學(xué)者發(fā)展了多種非線性隨機振動方法，如FPK方程法[8]、隨機平均法[9]、矩方程法[10]、等效線性化法[11-12]、等效非線性系統(tǒng)法[13]、概率密度演化法[14]和Wiener路徑積分法[15]等。其中，大部分方法僅限于解決小規(guī)模非線性問題，難以應(yīng)用于大跨度橋梁結(jié)構(gòu)[16]。等效線性化法因其較強的適用性在實際工程中得到廣泛的應(yīng)用，例如Colangelo[17]利用等效線性化法計算了多層彈塑性框架在一致平穩(wěn)地震激勵下的隨機響應(yīng)；趙巖等[18]采用虛擬激勵法求解非線性結(jié)構(gòu)的等效線性化運動方程，計算了連續(xù)梁橋在一致平穩(wěn)地震激勵下的彈塑性隨機響應(yīng)。盡管如此，該方法在求解非平穩(wěn)隨機問題時需要在不同時刻建立一系列等效線性系統(tǒng)，對于大規(guī)模工程結(jié)構(gòu)來說，計算量較大。該方法一般只能獲取結(jié)構(gòu)響應(yīng)的一、二統(tǒng)計矩，無法直接獲取非線性結(jié)構(gòu)響應(yīng)的尾部概率分布和峰值等更全面的統(tǒng)計信息。此外，隨機模擬法[19]是一種原理簡單且通用性強的方法，可以準(zhǔn)確獲取非線性結(jié)構(gòu)隨機響應(yīng)的任意階統(tǒng)計矩及概率分布，經(jīng)常用來檢驗其它方法的計算精度。但隨機模擬法需涉及到大量反復(fù)的確定性樣本分析，若采用傳統(tǒng)彈塑性時程分析法進(jìn)行樣本分析將會非常耗時。

在作者前期研究[20]基礎(chǔ)上，針對非一致地震激勵下局部進(jìn)入彈塑性狀態(tài)的大跨度橋梁結(jié)構(gòu)，基于大質(zhì)量法[21]提出了一種時域顯式降維迭代算法，可以快速進(jìn)行非一致地震激勵下的彈塑性時程分析。將該方法與蒙特卡羅法結(jié)合，提出了一種時域顯式降維迭代-隨機模擬法，可高效獲取非一致地震激勵下大跨度橋梁結(jié)構(gòu)彈塑性響應(yīng)的統(tǒng)計矩和平均峰值。以主跨為1 200 m的某大跨度懸索橋為工程算例，采用纖維梁柱單元模擬大橋彈塑性構(gòu)件，開展非一致地震激勵下彈塑性隨機振動分析研究，考察地震動空間效應(yīng)對彈塑性隨機響應(yīng)的影響。

1 非一致地震激勵

考慮地面運動加速度為均勻調(diào)制的非平穩(wěn)零均值隨機過程，則結(jié)構(gòu)支承處的非一致地面運動加速度向量可表示為

X(t)=Ψ(t)x(t)

(1)

式中

(2)

式中:m為非一致地震激勵個數(shù)；Xj(t)(j=1,2,…,m)為第j個非平穩(wěn)地面運動加速度；ψj(t)為反映Xj(t)非平穩(wěn)特性的均勻調(diào)制函數(shù)；xj(t)為Xj(t)的平穩(wěn)部分，它們的互功率譜矩陣Sx(ω)可表達(dá)為

(3)

式中：Sjj(ω)為xj(t)的自功率譜，Sjk(ω)(j≠k)為xj(t)與xk(t)間的互功率譜。

根據(jù)式(3)所示的功率譜矩陣，采用諧波合成法生成非一致非平穩(wěn)地面運動加速度為

(j=1,2,…,m)

(4)

式中:N為頻率分量個數(shù)；ωi為第i個頻率分量的圓頻率；φki為在(0,2π)范圍內(nèi)均勻分布的隨機相位角，且當(dāng)k≠r或i≠s時，φki和φrs相互獨立；ajk和θjk分別為第j與第k個地面運動加速度間的相關(guān)幅值與相關(guān)相位角。

對平穩(wěn)地面運動加速度x(t)的互功率譜矩陣Sx(ω)進(jìn)行Cholesky分解，可得

Sx(ω)=B(ω)B*T(ω)

(5)

式中

(6)

則式(4)中的ajk(ω)和θjk(ω)可分別表示為

(7)

(8)

式中:Δω為圓頻率間隔，函數(shù)Im、Re分別表示取復(fù)數(shù)的虛部、實部。

2 大質(zhì)量法

在一致地震激勵下，結(jié)構(gòu)內(nèi)力只與結(jié)構(gòu)相對地面的位移有關(guān)，因此一般基于結(jié)構(gòu)相對位移建立運動方程，其中地震激勵為結(jié)構(gòu)各質(zhì)點上一致的地面運動加速度。而在非一致地震激勵下，無法建立唯一的結(jié)構(gòu)-地面相對坐標(biāo)系，因此結(jié)構(gòu)運動方程及其地震激勵施加方式與一致地震激勵情況大為不同。常用的處理方法有相對運動法和大質(zhì)量法，其中相對運動法由于需要進(jìn)行疊加運算不適用于非線性情況，而大質(zhì)量法沒有這方面的限制，具有更強的適用性。

大質(zhì)量法是對結(jié)構(gòu)模型進(jìn)行動力等效的一種分析方法，該方法在處理非一致地震激勵問題時需要解除支承點沿地震作用方向的約束，并賦予該支承點大質(zhì)量，此時彈塑性結(jié)構(gòu)的運動方程可寫為

(9)

式中

(10)

(11)

3 時域顯式降維迭代法

對于局部彈塑性結(jié)構(gòu)，將結(jié)構(gòu)非線性恢復(fù)力向量F(U)分解為線性部分和僅與彈塑性單元恢復(fù)力相關(guān)的部分為

F(U)=K0U+LnFn

(12)

式中:K0為結(jié)構(gòu)初始彈性剛度陣；Ln為彈塑性單元恢復(fù)力的指示矩陣；Fn為彈塑性單元恢復(fù)力中的非線性部分可表示為

Fn=R(Un)-Kn,0Un

(13)

式中:R為彈塑性單元的恢復(fù)力向量；Un為彈塑性單元的節(jié)點位移向量；Kn,0為彈塑性單元的初始彈性剛度矩陣，屬于結(jié)構(gòu)初始彈性剛度矩陣K0的一部分。將式(12)代入式(9)，并將非線性部分移至方程右端，可得

(14)

式中

(15)

L=[Lx-Ln]

(16)

Vi=Ai,0P0+Ai,1P1+…+Ai,iPi

(i=1,2,…,nt)

(17)

式中:Vi=V(ti)，Pi=P(ti)，i=1,2,…,nt，其中nt為總時間步數(shù)；Ai,j(j=0,1,…,i)為系數(shù)矩陣，它們只與結(jié)構(gòu)參數(shù)有關(guān)，可表達(dá)為如下閉合公式：

(18)

基于Newmark-β積分方法可獲取T,Q1和Q2為

(19)

式中:I為單位矩陣；取γ=0.5,β=0.25,則Newmark-β積分方法可以實現(xiàn)無條件穩(wěn)定。

對于當(dāng)前時刻ti，假設(shè)已逐步計算出結(jié)構(gòu)響應(yīng)V0,V1,…,Vi-1，則式(17)右端的為已知量P0,P1,…,Pi-1，然而P1中依然含有關(guān)于未知向量Un,i的非線性項，需要進(jìn)行迭代計算。式(17)共有2nd(nd為結(jié)構(gòu)總自由度數(shù))維，對于大規(guī)模結(jié)構(gòu)，直接迭代求解該式的計算量較大。為了提高計算效率，從式(17)中提取Un,i相關(guān)的行，可得

(20)

Ψ(Un,i)=Yi(i=1,2,…,nt)

(21)

式中

(22)

(i=1,2,…,nt)

(23)

采用Newton-Raphson方法求解式(21)，將所獲得的收斂結(jié)果Un,i代入式(17)的其余行中，則得到結(jié)構(gòu)當(dāng)前時刻其它單元的節(jié)點位移向量。例如對于結(jié)構(gòu)某一關(guān)鍵響應(yīng)r，可表示為

(24)

4 時域顯式降維迭代-隨機模擬法

將時域顯式降維迭代法與蒙特卡羅法相結(jié)合，利用時域顯式降維迭代法進(jìn)行樣本分析，高效獲取非一致地震激勵下大跨度橋梁結(jié)構(gòu)的彈塑性隨機響應(yīng)，具體計算步驟如下：

(2)建立結(jié)構(gòu)初始有限元模型，釋放地震激勵方向的支承約束，并在相應(yīng)自由度施加大質(zhì)量。

5 工程實例

5.1 有限元模型

虎門二橋大沙水道橋是一座主跨1 200 m的地錨式懸索結(jié)構(gòu)。該橋跨徑布置為(360+1 200+480)m，矢跨比1∶9.5，如圖 1所示。門式塔高191.1 m，設(shè)上、下兩道橫梁，如圖2所示。主塔為鋼筋混凝土結(jié)構(gòu)，塔底截面及其配筋情況如圖 3所示，其中縱向受力鋼筋強度等級為HRB400，混凝土強度等級為C50。利用ANSYS建立大橋的初始三維有限元模型。其中，主梁、主塔和邊墩均采用梁單元模擬，主纜和吊桿則通過桿單元模擬。主塔和邊墩底部固支，主纜通過大橋的兩側(cè)錨碇固定。主梁側(cè)向和豎向運動受主塔和邊墩的約束，而縱向運動不受主塔和邊墩約束。

圖1 懸索橋立面圖

圖2 主塔立面圖

圖3 主塔塔底截面及其配筋圖

歷史震害數(shù)據(jù)表明，大跨度懸索橋主塔相對其它部位更容易損傷[22]。虎門二橋大沙水道橋抗震專題研究表明，在重現(xiàn)期分別為945年和2450年的中震和大震作用下，整個結(jié)構(gòu)可以依然處于彈性狀態(tài)，然而主塔塔底截面承載能力的冗余度遠(yuǎn)低于其它關(guān)鍵截面。因此，當(dāng)發(fā)生強度高于大震的地震時，主塔塔底截面可能會率先進(jìn)入彈塑性狀態(tài)，應(yīng)該采用彈塑性單元模擬主塔底部，并預(yù)設(shè)塑性發(fā)展范圍為塔底10 m以內(nèi)，在后續(xù)分析中將對該塑性發(fā)展范圍進(jìn)行驗證。在地震和恒載作用下，主塔底部主要受到軸力、剪力和彎矩作用，由于不允許發(fā)生剪切破壞，在此僅考慮軸力和彎矩引起的塑性行為，并采用可以模擬軸力-彎矩耦合作用的纖維彈塑性單元對預(yù)設(shè)的塔底塑性范圍進(jìn)行模擬。分別采用Menegotto-Pinto模型[23]和修正Kent-Park模型[24]模擬鋼筋纖維和混凝土纖維的軸向應(yīng)力應(yīng)變關(guān)系，鋼筋的屈服強度和混凝土的抗壓強度分別取為330 MPa和22.4 MPa。本文通過幾何非線性靜力分析計算恒載引起的大位移和應(yīng)力剛化效應(yīng)，在此基礎(chǔ)上開展地震激勵下的結(jié)構(gòu)動力分析，不考慮地震激勵新引起的幾何非線性效應(yīng)。

5.2 非一致地震激勵

考慮東、西塔底地震激勵的差異性，而同側(cè)塔底與錨碇相距較近，因此可以認(rèn)為同側(cè)塔底與錨碇的地震激勵具有一致性。為此，該結(jié)構(gòu)共有2個非一致地震激勵，式(1)中反映地面運動加速度非平穩(wěn)的均勻調(diào)制函數(shù)可表示為

(25)

式中:ta=3.2 s，tb=17.0 s,c=0.13。地震激勵間的互功率譜可表示為

(26)

圖4 非一致地面運動加速度時程樣本

5.3 彈塑性時程分析

針對圖 4所示的非一致地面運動加速度時程樣本，分別采用本文的時域顯式降維迭代法和傳統(tǒng)非線性Newmark-β法對大橋進(jìn)行彈塑性時程分析。東塔底部截面Ⅰ和下橫梁上截面Ⅳ的彎矩時程分別如圖 5和圖 6所示。從圖中可以看出，兩種方法的結(jié)果完全重合，驗證了時域顯式降維迭代法的準(zhǔn)確性。

圖5 東塔截面Ⅰ的彎矩時程

圖6 東塔截面Ⅳ的彎矩時程

傳統(tǒng)非線性Newmark-β法和時域顯式降維迭代法的計算時間分別為2 518.3 s和127.4 s。顯而易見，本文的時域顯式降維迭代法具有更高的計算效率，這是因為該方法僅僅需要對主塔底部4個彈塑性纖維單元進(jìn)行非線性降維迭代運算，而傳統(tǒng)方法需要將橋梁所有自由度捆綁一起進(jìn)行迭代計算。需要指出的是，本文方法的計算時間主要包括兩部分：一部分用于構(gòu)建如式(17)所示動力響應(yīng)顯式表達(dá)式中的系數(shù)矩陣，為109.9 s；另一部分用于降維迭代運算，僅為17.5 s。其中，系數(shù)矩陣只依賴于結(jié)構(gòu)本身信息，它們只需要計算一次，隨后可重復(fù)用于后續(xù)隨機模擬的所有樣本分析中。

為了考察強震作用下主塔塑性變形的發(fā)展情況，給出如圖 2所示東塔底部截面Ⅰ-Ⅲ的彎矩-曲率關(guān)系，分別如圖7(a)～(c)所示。從圖中可看出，非一致地震激勵下截面Ⅰ的彎矩-曲率曲線呈現(xiàn)出相對飽滿的旗形滯回環(huán)，表明截面Ⅰ表現(xiàn)出顯著的彈塑性行為；截面Ⅱ的滯回環(huán)面積明顯減小，非線性行為相對較弱；而截面Ⅲ的彎矩-曲率關(guān)系幾乎趨近于線彈性。以上結(jié)果表明，長度為10 m的彈塑性纖維單元可以完全捕捉主塔底部的彈塑性變形行為。截面Ⅰ-Ⅲ最外側(cè)鋼筋應(yīng)力-應(yīng)變關(guān)系如圖8(a)～(c)所示，同樣表明在所關(guān)注的地震水準(zhǔn)下，主塔底部彈塑性行為發(fā)展未超出預(yù)期范圍。

圖7 東塔底部截面Ⅰ-Ⅲ的彎矩-曲率關(guān)系

圖8 東塔底部截面Ⅰ-Ⅲ外側(cè)鋼筋應(yīng)力-應(yīng)變關(guān)系

5.4 彈塑性隨機振動分析

時域顯式降維迭代法用于彈塑性時程分析的計算精度和計算效率均得到了驗證，可進(jìn)一步利用時域顯式降維迭代-隨機模擬法開展彈塑性抗震隨機振動分析。為了考察地震激勵非一致性的影響，同時計算懸索橋在一致和非一致順橋向地震激勵下的彈塑性隨機響應(yīng)，其中非一致地震激勵按5.2節(jié)中的參數(shù)確定，一致地震激勵取其東塔底處的激勵，隨機模擬樣本數(shù)取為500。

兩種地震激勵工況下西塔和東塔底部截面內(nèi)力標(biāo)準(zhǔn)差分別如圖9和圖10所示。從圖中可以觀察到，對于西塔底截面，非一致地震激勵下的軸力標(biāo)準(zhǔn)差偏小，而剪力和彎矩標(biāo)準(zhǔn)差偏大；對于東塔底截面，非一致地震激勵下的軸力和彎矩標(biāo)準(zhǔn)差偏小，剪力標(biāo)準(zhǔn)差與一致地震激勵下的結(jié)果非常接近。由此可知，地震激勵空間效應(yīng)會使得該大跨度懸索橋內(nèi)力結(jié)果既可能偏大也可能偏小。

圖9 西塔底截面內(nèi)力標(biāo)準(zhǔn)差

圖10 東塔底截面內(nèi)力標(biāo)準(zhǔn)差

此外，利用本文方法可進(jìn)一步獲取大橋關(guān)鍵響應(yīng)的平均峰值。主塔底截面軸力、剪力和彎矩的平均峰值如表1所示。為了考察地震激勵空間效應(yīng)對結(jié)構(gòu)響應(yīng)的影響，將兩種工況的相對差異也列于表1中。從表中可以看出，地震激勵空間效應(yīng)對響應(yīng)平均峰值的影響規(guī)律與對響應(yīng)標(biāo)準(zhǔn)差的影響規(guī)律基本相同。值得注意的是，當(dāng)考慮地震激勵空間效應(yīng)時，東塔底截面軸力的平均峰值比一致地震激勵下的結(jié)果小了9.4%；西塔底截面剪力和彎矩的平均峰值比一致地震激勵下的結(jié)果分別大了12.4%和11.8%。在該大橋彈塑性抗震設(shè)計中需要考慮地震激勵空間效應(yīng)對結(jié)構(gòu)內(nèi)力的影響。

表1 塔底截面內(nèi)力的平均峰值

時域顯式降維迭代-隨機模擬的計算時間如表2所示，取不同樣本數(shù)，平均每個樣本分析的時間基本相同，均是13 s左右。這表示只需要花13 s的時間就可以完成一次大跨度懸索橋確定性彈塑性時程分析。當(dāng)樣本數(shù)為500時，時域顯式降維迭代-隨機模擬法的計算時間為不到2個小時，而基于傳統(tǒng)非線性Newmark-β方法的隨機模擬法的計算時間則接近200個小時。本文的時域顯式降維迭代-隨機模擬法在計算效率方面具有明顯的優(yōu)越性。以上所有的計算均是在帶有Intel Core i5處理器和4 GB內(nèi)存的臺式機上完成的。

表2 時域顯式降維迭代-隨機模擬的計算時間

6 結(jié) 論

(1)基于大質(zhì)量法推導(dǎo)了非一致地震激勵下結(jié)構(gòu)動力響應(yīng)形式上的時域顯式表達(dá)式，利用該顯式表達(dá)式的降維列式優(yōu)勢，提出僅關(guān)于彈塑性單元節(jié)點自由度的時域顯式降維迭代算法，顯著提高非一致地震激勵下大跨度橋梁彈塑性時程分析的計算效率。并結(jié)合蒙特卡羅法，提出了非一致地震激勵下大跨度橋梁彈塑性隨機振動分析的時域顯式降維迭代-隨機模擬法。

(2)以虎門二橋大沙水道橋為工程實例，開展順橋向非一致地震激勵下的彈塑性隨機振動分析，驗證了本文方法的正確性和高效性。研究結(jié)果表明，對于該大跨度懸索橋，非一致地震激勵下主塔內(nèi)力的標(biāo)準(zhǔn)差和平均峰值比一致地震激勵下的結(jié)果既可能偏大也可能偏小，其中東塔底截面軸力的平均峰值偏小了9.4%，西塔底截面的剪力和彎矩平均峰值分別偏大了12.4%和11.8%。在該大橋彈塑性抗震設(shè)計中需要考慮地震激勵空間效應(yīng)對結(jié)構(gòu)內(nèi)力的影響。