多點地震動激勵下的高效反應譜方法＊

王君杰郭進

1） 中國上海 200092 同濟大學土木工程學院

2） 中國石家莊 050043 石家莊鐵道大學土木工程學院

引言

地震動是空間變化的，由于波列傳播速度的有限性和相干性損失以及局部場地地質條件不同，結構各支承點的地震激勵出現顯著差異，這種差異可能對平面大尺度結構（例如橋梁、水壩）的地震反應產生重要影響.歐洲規范 8 （European Committee for Standardization，2004）和城市軌道交通結構抗震設計規范GB 50909—2014 （中華人民共和國住房和城鄉建設部，2014）也提到了地震動空間變化對結構的影響.地震動空間變化可以通過以下三種方法來分析：① 時間歷程方法；② 隨機振動方法；③ 反應譜方法.后兩種方法的優點是應用了地震動的統計特征對結構響應進行分析.

對于隨機振動分析，地震動輸入由自功率譜密度矩陣（auto-power spectral density，縮寫為APSD）確定；對于反應譜分析，地震動輸入由地震反應譜確定，模態組合系數通常由結構隨機振動理論推導得出（王君杰，1992；Berrah，Kausel，1992；Kiureghian，Neuenhofer，1992；王君杰等，1995）.一般情況下，隨機振動法的數值積分運算量很大，基于隨機振動理論的反應譜組合系數的計算也非常耗時.為減少積分時間，Harichandran和Vanmarcke （1986）提出了行波激勵的解析公式，Loh和Ku （1995）基于簡化的自功率譜密度函數提出了一種有效的計算方法，但采用了簡單的空間相干函數.孫建梅等（2003）對Kiureghianh和Neuenhofer （1992）提出的多點輸入反應譜方法進行了簡化，即在相關系數的計算中忽略頻率比介于0.7—1.2范圍以外的部分，減少了計算量，但損失了精度.

為提高多點地震動激勵下譜分析方法中相關系數的計算效率，本文擬基于對自功率譜和相關系數積分表達式的近似處理，提出相關系數的近似解析公式.

1 結構地震響應

假定在空間某點r（r＝ 1，2，···，m）處的地震激勵是零均值平穩過程，其位移、速度和加速 度 函 數 分 別 記 為u（rt），u˙ （rt）和u¨（rt）， 結 構 平 穩 響 應 的 方 差 為 （王 君 杰 ，1992；Kiureghian，Neuenhofer，1992）

式中：σk為第k個結構反應量的均方差；σgr和σgs分別為第r個空間點和第s個空間點的地面運動位移的均方差；σir和σjs分別為第r個地面運動激勵下結構第i振型反應的均方差和第s個地面運動激勵下結構第j振型反應的均方差；參數a和b均為與結構有關的系數；m為多點地震動的數目；n為參與計算的結構的振型數目；ρgrgs，ρgrjs和ρirjs定義為

其中

式中，H（iω）為結構自由振動第i階振型的頻響函數，）為第r空間點地震加速度與第s空間點地震加速度的互功率譜密度函數，ω為圓頻率.根據式（1），多點地震激勵下結構反應的組合方法為

式中，SD（rωi，ζi）和SD（sωj，ζj）分別代表第r個地震激勵下i階振型最大位移和第s個地震激勵下j階振型最大位移，即平均位移反應譜.同時，ur，max＝SD（r0，ζi），us，max＝SD（s0，ζj）.

對于式（3），根據復相干函數，互功率譜密度函數Sr（sω）可以寫為

式中，γ為相關因子， ρr（sω，Δrrs）為地震動的空間相干函數， θr（sω，Δrrs）為地震動的相位差，Δrrs為空間距離的向量.

國內外研究人員已經提出了多個 ρr（sω，Δrrs）表達式，例如Harichandran和Vanmarcke（1986）及Qu等（1996）提出的表達式等.Qu等（1996）提出的地震動空間相干函數模型如圖1所示，其經驗表達式為

圖1 空間相干函數模型（Qu et al，1996）（a） ρrs-空間距離關系；（b） ρrs-頻率關系Fig.1 Spatial coherence function model （Qu et al，1996）（a） ρrs-spatial distance relationship；（b） ρrs-frequency relationship

式中：a1＝1.678×10-5，a2＝1.219×10-3；b1＝-5.5×10-3，b2＝0.7674.

Harichandran和Vanmarcke（1986）提出的地震動空間相干函數模型如圖2所示，其經驗表達式為

圖2 空間相干函數模型（Harichandran，Vanmarcke，1986）（a） ρrs-空間距離關系；（b） ρrs-頻率關系Fig.2 Spatial coherence function model （Harichandran，Vanmarcke，1986）（a） ρrs-spatial distance relationship；（b） ρrs-frequency relationship

式中A， α ，K，θ， ω0和b為經驗參數.對于臺灣SMART-1臺網第24次地震，A＝0.736，α＝0.147，K＝5210， ω0＝6.85 rad/s，b＝2.78.

本文中使用到的兩個自功率譜函數經驗模型見圖3.第一個是Hu模型（胡聿賢，周錫元，1962，1965；王君杰，1992；王君杰，江近仁，1997），即

圖3 自功率譜密度函數 S r（rω）的特征（ζg 為場地土的阻尼比，下同）（a） 胡聿賢自功率譜密度模型；（b） 克拉夫-彭津自功率譜密度模型Fig.3 The characteristic for auto-power spectral density function（ζg is the damping ratio of site soil，the same below）（a） Hu’s auto-power spectral density model；（b） Clough-Penzien’s auto-power spectral density model

第二個是克拉夫-彭津模型（Clough，Penzien，1993），即

本文在計算相關系數時，將以上兩個模型分別簡化為

式中，S0為APSD的強度因子， ωg為 場地土層的頻率參數， ζg為場地土層的阻尼參數， ωc為胡聿賢模型在低頻段的濾波參數， ωf和 ζf為Clough-Penzien模型在低頻段的濾波參數.計算中取ωg＝1.5 Hz， ωc＝0.3 Hz， ωf＝0.25 Hz， ζg＝0.6， ζf＝0.4，S0＝1.

2 近似自功率譜密度函數的精度

分別比較相干系數 ρgrgs， ρgrjs和 ρirjs的近似值［ 式 （10）］ 與 精確值［ 式 （8）］ 之 間的差異，結果如圖4—6所示.其中，空間相干函數使用式（6），視波速向量vapp的模取1 000 m/s.

圖4 系數 ρ grgs的精確值與近似值的比較Fig.4 The comparison between exact value and approximate value of coefficient ρgrgs

從圖4可以看出，式（10）給出的系數 ρgrgs幾乎是精確值.圖5表明：當 ωj處于低頻段時，系數 ρgrjs非常大；隨著 ωj的增大及r點與s點之間空間距離的增大， ρgrjs逐漸減小至零，直至可以忽略，說明式（10）可以給出具有高精度的 ρgrjs近似解答.

圖5 |Δ r rs|＝300m （a）和1 000 m （b）時系數 ρ grjs 的精確值與近似值的比較（ζj ＝ 0.05）Fig.5 The comparison between exact value and approximate value of coefficient ρg r js （ζj ＝ 0.05）with |Δ r rs|＝300m （a） and 1 000 m （b）

從圖6a-d可以看出：當 ωi和 ωj位于低頻段時，由式（10）計算得出的相干系數 ρirjs近乎精確值， 而 近 似的 ρirjs在 ωi和 ωj較 大時會 出 現 更 大的誤差；在高頻段內，地震動自功率譜的值比低頻段的值小得多，相比之下重要性不大.同時，系數 ρirjs的值隨著 ωi和 ωj的增加而減小.以上結果表明，相干系數 ρirjs的近似值在工程實踐中具有很好的準確性.

圖6 系數 ρ irjs （ζi = ζj=0.05）精確值與近似值的比較Fig.6 The comparison between exact value and approximate value of coefficient ρi rj s （ζi = ζj=0.05）（a） |Δ r rs|＝300m，ωi＝0.1 Hz；（b） |Δ r rs|＝1000m，ωi＝0.1 Hz；（c） |Δ r rs|＝300m，ωi＝1.0 Hz；（d） |Δ r rs|＝1000m，ωi＝1.0 Hz；（e） |Δ r rs|＝300m ，ωi＝4.0 Hz；（f） |Δ r rs|＝1000m，ωi＝4.0 Hz

3 相關系數的近似公式

3.1 近似空間相干函數 ρr（sω，Δrrs）

如前所述，由于空間相干函數的經驗表達式 ρr（sω，Δrrs）的復雜性，需要進行數值積分來計算系數 ρgrgs， ρgrjs和 ρirjs.為了得到這三個系數的近似解析表達式，引入近似的 ρr（sω，Δrrs）.

組合系數 ρgrgs與擬靜態響應有關.擬靜態響應是由位移貢獻的，它由地震動中非常低的頻率成分控制.因此，計算系數 ρgrgs的相干函數近似表達式為

式中， ωd，max是地震動中位移的APSD達到最大值時對應的圓頻率.

組合系數 ρgrjs與式（1）或式（4）中的第二項耦合響應項有關，此耦合項由地震動中低頻成分控制.因此，計算系數 ρgrjs的相干函數近似表達式被定義為

式中， ωgrjs取或（ ωj＋ωd，max）/2.

組合系數 ρirjs與式（1）或式（4）中的第三項動力響應項有關，動力響應項由地震動中靠近以及介于 ωi和 ωj之間的頻率成分控制.因此，計算系數 ρirjs的相干函數近似表達式為

從圖7—9所示結果可以看出，使用式（6）和式（8）計算得出的上述近似相干函數 ρr（sω，drs）具有很好的準確性.由此可知，本文所提方法在ρr（sω，Δrrs）的近似表達上具有良好的精度.

圖7 系數 ρ grgs的精確值與近似值的比較Fig.7 The comparison between exact value and approximate value of coefficient ρgrgs

圖8 |Δ rr s|＝300m （a）和1 000 m （b）時系數 ρ grjs （ζj ＝ 0.05）精確值與近似值的比較Fig.8 The comparison between exact value and approximate value of coefficient ρg r js （ζj ＝ 0.05）with |Δ rr s|＝300m （a） and 1 000 m （b）

經過上述簡化后，三個相關系數 ρgrgs， ρgrjs和 ρirjs可以簡化為

對于胡聿賢自功率譜密度模型有

圖9 系數 ρ irjs （ζi ＝ ζj＝0.05）精確值與近似值的比較Fig.9 The comparison between exact value and approximate value of coefficient ρi rj s （ζi ＝ ζj＝0.05）（a） |Δ r rs|＝300m，ωi＝0.1 Hz；（b） |Δ r rs|＝1000m，ωi＝0.1 Hz；（c） |Δ r rs|＝300m，ωi＝1.0 Hz；（d） |Δ r rs|＝1000m，ωi＝1.0 Hz；（e） |Δ r rs|＝300m，ωi＝4.0 Hz；（f） |Δ r rs|＝1000m，ωi＝4.0 Hz

式中： R e 為復數的實部； τ為地震波從r點傳播到s點的時間，，當r與s在相同位置時， τ等于零.

對于克拉夫-彭津的自功率譜密度模型則有

3.2 胡聿賢自功率譜密度模型下相關系數的解析表達式

根據復變函數的基本理論，式（16）的積分形式為

式中， R e［s（fzp）］ 為復變函數的留數，zp是函數f（z）的極點，Imzp是zp的虛部.

根據式（18）得到

如果 τ ≥0 ，則

如果 τ ＜0，則

如果 τ ≥0 ，則

如果 τ ＜0，則

3.3 Penzien-Clough 自功率譜密度下相關系數的解析表達式

式（16）的積分形式為

如果 τ ≥0 ，則

如果 τ ＜0 ，則

如果 τ ≥0 ，則

其中，z1＝ωjeiθj，z2＝ - ωje-iθj，z3＝ωfe-iθf，z4＝ - ωfeiθf；z5＝ωfeiθf，z6＝ - ωfe-iθf，z7＝ωie-iθi，z8＝-ωieiθi；這里取下半平面的單極點留數 R e［s（fz3）］ ， R e［s（fz4）］ ， R e［s（fz7）］ 和 R e［s（fz8）］.

如果 τ＜ 0，則

以上諸式中， θi＝arcsinζi， θj＝arcsinζj， θf＝arcsinζf.

4 精度和計算效率的討論

4.1 相關系數解析表達式的精度

本節以一座斜拉橋為例進行計算，其中跨240 m，邊跨117.5 m，其有限元模型示意圖如圖10所示.假設地震波沿縱向傳播，視波速取1 000 m/s.用式（9）或式（11）的自功率譜密度函數以及式（6）的地震動空間相干模型，計算該橋梁在多點地震激勵下的反應.

圖10 算例橋梁的有限元模型Fig.10 Finite element model of the example bridge

計算中取前50階振型.表1—4分別列出了橋梁不同構件中某些位置處的絕對位移、相對位移、軸力、剪力和彎矩.表中：Ln（n＝1—6）代表橋梁左側的反應，Rn代表右側的反應；n＝1時代表塔頂的絕對位移；n＝2時代表塔頂和塔底之間的相對位移；n＝3時代表邊墩頂和梁端之間的相對位移，這項參數在落梁問題中很重要；n＝4時代表塔底的軸力；n＝5時代表塔底的剪力；n＝6時代表塔底的彎矩.

表1 Qu 等（1996）相干模型下的地震位移響應及相對誤差Table 1 Seismic displacement response and its relative error under the coherence model of Qu et al （1996）

定義兩個相對誤差：

式中：A表示基于精確APSD模型和精確空間相干模型的響應的數值積分解；B表示基于精確APSD模型和近似空間相干模型的響應的數值積分解；C為3.2和3.3節中提出的解析解；eBA和eCA是相對誤差.

表2 Harichandran 和 Vanmarcke （1986）相干模型下的地震位移響應及相對誤差Table 2 Seismic displacement response and its relative error under the coherence model of Harichandran and Vanmarcke （1986）

在方法A和B中，使用了變步長的辛普森積分法.在積分極限從1 000，500，200到100 rad/s的情況下，將基于胡聿賢的功率譜密度函數與屈鐵軍、王君杰的空間相干函數的系數 ρgrgs， ρgrjs， ρirjs進行了比較，最終方法A采用200 rad/s.當相對積分精度分別為1 0-8，1 0-5，10-4和1 0-3時，分別計算出系數 ρgrgs， ρgrjs， ρirjs，并采用了精度為1 0-4時的積分步長.

表3 Qu 等（1996）相干模型下的地震力或彎矩響應及相對誤差Table 3 Seismic force or moment response and its relative error under the coherence model of Qu et al （1996）

表4 Harichandran 和 Vanmarcke （1986）相干模型下的地震力或彎矩響應及相對誤差Table 4 Seismic force or moment response and its relative error under the coherence model of Harichandran and Vanmarcke （1986）

從表1—4可以看出，相對誤差非常小，這說明近似空間相干函數對結構的反應影響非常小.從3.1節中的圖7—9也可以得出這個結論.

除了塔頂的絕對位移（即L1和R1）外，A和C方法的反應相對誤差都小于1%.同時在地震分析中，絕對位移并不重要，而相對位移卻很重要.這表明解析表達式可以與數值方法很好地吻合.

4.2 計算效率的比較

表5列出了數值方法A和解析方法C對該例橋梁的計算時間，計算中取前50階振型.表6列出了計算中不同階振型的計算時間.在本例中，只計算了以下反應：① 四個塔頂的絕對位移；② 四個塔頂與塔底之間的相對位移；③ 兩個邊墩頂部與梁端之間的相對位移；④ 四個塔底的軸力、剪力和彎矩.在本例中，使用英特爾Visual Fortran 11.1.067編寫有限元程序，前50階的特征值計算耗時50 s.

從表5可以發現，解析法求解的時間消耗只有數值求解的1/52左右.就絕對時間而言，解析法需要93—111 s，這在實際設計中是完全可以接受的.同時，數值解決方案需要5 833—6 273 s，比解析法要多耗時近98%.從表6中也可以得出同樣的結論，計算時間隨著所取振型階數的增加而增加.對于大型橋梁，在計算中可能需要取數百階振型.與用數值法相比，采用解析法計算大型結構響應可大大縮短計算用時.

表5 計算時間比較（單位：s）Table 5 Comparison of time consumption of computation （Unit：s）

表6 計算時間比較（單位：s）Table 6 Comparison of time consumption of computation （Unit：s）

5 結論

本文在提出近似空間相干函數的基礎上，建立了相關系數 ρgrgs， ρgrjs和 ρirjs的解析表達式，一方面保證了計算結果的精度，另一方面采用解析法可以更快速地獲得三個相關系數的積分值.同時使用了僅需要一至兩個地震動參數的APSD近似表達式.通過對數值計算結果的分析，本文所提相關系數 ρgrgs， ρgrjs和 ρirjs的解析式有足夠的工程精度和較高的計算效率，可為平面大尺度結構多點激勵下的地震響應計算提供高效方法.