基于單位脈沖響應的人造地震波模擬和結構抗震隨機模擬分析

李雪平 張磊源 魏鵬

(華南理工大學 土木與交通學院∥亞熱帶建筑科學國家重點實驗室，廣東 廣州510640)

地震動具有很強的非平穩性，地震動隨機模型是應用隨機振動理論研究結構隨機地震反應的基礎。20世紀50年代，Housner[1]首次提出利用理想白噪聲模型來模擬地震動地面加速度過程。地震動隨機模型發展至今，已經建立了一些較好的隨機模型，如平穩模型中的Kanai-Tajimi模型[2]、Clough-Penzien模型[3]等，非平穩模型中的強度非平穩模型[4]等。在應用地震動隨機模型時，選擇合理的模型設計參數對于計算分析具有重要意義。歐進萍等[5]、洪峰等[6]分別利用隨機極值理論、強震記錄統計分析的結果來確定模型參數。薛素鐸等[7]、任文杰等[8]根據我國的抗震規范進行研究，并給出了隨機地震動模型參數。

結構動力時程分析的可靠性不僅取決于合理的結構計算模型和精確的計算方法，而且與地震波的輸入選擇有重要關系[9- 11]。賈宏宇等[12]、黃朝光等[13]采用三角級數法合成人造地震波，該法原理簡單,但所得地震波的平均反應譜與規范反應譜有一定差距；江近仁等[14]利用由目標反應譜轉換的功率譜合成人造地震波，為了提高精度，需要多次迭代，計算效率較慢；丁佳偉等[15]提出一種能同時滿足功率譜與反應譜精度的時頻域雙控的合成地震波方法；王向超等[16]將數值模擬方法與地震動預測方程相結合，為重大水電工程抗震分析提供了一種地震波合成方法。受限于已有地震動隨機功率譜模型的固定形式，通過其模擬的地震波與規范的反應譜存在很大的差異。蘇成等[17]直接采用現行《建筑抗震設計規范》(GB 50011—2010)中的規范反應譜，采用諧波合成法[18]模擬了與規范反應譜相一致的人造地震波。

振型分解反應譜法是我國現行規范規定的抗震計算方法，時程分析法是對振型分解反應譜法計算結果的補充計算?，F行規范規定，當采用時程分析法時，多條時程曲線計算所得結構底部剪力的平均值不得小于振型分解反應譜法計算結果的80%。

顯式時域隨機模擬分析是一種針對大型復雜結構的高效非平穩隨機振動分析方法[19- 23]，其從工程技術人員較為熟悉的時域分析方法入手，更加簡單、高效。以往的顯式時域法(ETDM)都是基于精細積分法、Newmark-β法等數值積分方法推導而來，而本文基于單位脈沖的解析響應推導了單自由度和多自由度動力系統的顯式時域表達式；然后將單自由度動力系統的顯式時域表達式應用于求解地震動隨機模型過濾白噪聲獲得人造地震波，利用其能快速求解響應的特點，通過迭代，確定與規范反應譜相符的譜強度因子，進而高效地模擬指定強度和與規范反應譜相符的人造地震波；最后將本文方法模擬的人造地震波和多自由度動力系統的顯式時域表達式應用于結構抗震隨機模擬分析，并通過算例進行驗證。

1 基于單位脈沖的顯式時域表達

1.1 單自由度系統

受動力荷載u(t)作用下的單自由度體系的運動微分方程為：

(1)

將荷載u(t)離散為時刻t0、t1、…、tl處的激勵值u0、u1、…、ul，時間步長為Δt。通過數值積分的方法可以把式(1)的解寫為離散的遞推表達式：

(2)

(3)

當u0=1且ui=0(1≤i≤l),即系統(1)在t0時刻作用單位脈沖激勵時，各時刻響應為：

(4)

Zi=TiRi, 1≤i≤l

(5)

式中：

(6)

(7)

e-ζωnti-1(ωDcos(ωDti-1)+ζωnsin(ωDti-1))]/

(8)

則ti時刻的速度表達式為：

(9)

1.2 多自由度系統

(10)

文中采用Rayleigh阻尼模型，對式(10)應用振型分解法，取振型截斷數為Nd，則Y(t)可表示為：

(11)

式中：Φ為振型矩陣，{q(t)}為廣義坐標向量，且

(12)

將式(12)代入式(11)，并在方程兩邊左乘ΦT，則式(10)可被解耦成Nd個單自由度體系動力方程：

(13)

根據單自由度結構基于單位脈沖的顯式時域表達，第k個單自由度體系響應的顯式時域表達為：

(14)

式中，ak,i(1≤i≤l)為相應的二維列向量系數。

將系數ak,i寫成：

(15)

式中,bk,i和ck,i分別為ak,i的第1個元素和第2個元素。

那么，第k個單自由度體系在ti時刻的位移響應和速度響應的顯式時域表達式可分別表示：

(16)

將式(16)按式(11)進行振型疊加，可以得到多自由度線性結構ti時刻響應的顯式時域表達式為：

(17)

(18)

若令

(19)

則式(17)可進一步表示為：

Qi=HiRi, 1≤i≤l

(20)

式(20)即為多自由度線性結構基于單位脈沖與振型分解法的顯式時域表達式。

(21)

相比式(20)，采用式(21)計算動力響應可以實現多自由度系統的降維計算，同時也可以減少系數矩陣的儲存量,在進行大規模結構動力時程分析時可以顯著減少計算量，大幅提升計算效率。

2 地震動隨機模型及設計參數

工程中常用的非平穩地震模型為強度非平穩模型[4]，它是將平穩隨機過程乘上均勻調制函數來模擬非平穩隨機過程。該非平穩模型地震加速度可表示為：

(22)

本文采用Kanai-Tajimi過濾白噪聲模型(以下簡稱為K-T模型)進行研究，該模型的地面加速度功率譜密度為：

(23)

式中：ωg、ζg為濾波器的參數[25]；S0為白噪聲的譜強度因子，它是直接與地震動強度有關的參數。

圖1 Kanai-Tajimi模型

該模型的覆蓋土層運動方程為：

(24)

那么，平穩地面加速度為：

(25)

式(22)中的均勻調制函數采用工程中常用的均勻調制函數[26]，其表達式為：

(26)

式中,t1、t2為主震平穩段的起始時間,c為衰減系數,這些參數取值與震級和場地條件等有關。

3 與規范反應譜相符的人造地震波

3.1 指定地震強度的人造地震波的模擬

采用基于單位脈沖的顯式時域法求解運動方程(24)，有：

(27)

式中：

將式(27)代入式(25)可得平穩地面加速度：

(28)

將式(28)乘上均勻調制函數g(t)，可得人造地震波的表達式為：

(29)

參考文獻[2]，在表1中給出了當基于K-T模型模擬指定強度的地震波時，在多遇地震、設計地震第一組的抗震設計條件下譜強度因子S0的取值。

表1 K-T模型的譜強度因子S0的取值Table 1 Spectral intensity factor S0 reference value of K-T model 10-4 m2/s3

采用基于單位脈沖的顯式時域法,求解K-T模型模擬樣本數為10 000條的人造地震波樣本。限于篇幅，本文只列出抗震設計條件為8度多遇地震、設計地震分組第一組、場地類別分別為Ⅰ0、Ⅱ、Ⅲ時的計算結果。地震動持時依次取為15、15和20 s，時間步長取為0.02 s。

圖2給出了在Ⅱ類場地下的模擬人造地震波樣本里的單條地震波時程曲線。圖3為在Ⅱ類場地下模擬的10 000條人造地震波的樣本方差和K-T模型目標方差的比較，其中目標方差可根據K-T譜的自相關函數[27]求得。為了驗證本文方法(單位脈沖+ETDM)的正確性，同時給出了Newmark-β直接法和基于Newmark-β法的顯式時域法[28]模擬人造地震波的計算結果。各種方法的計算耗時和最大相對誤差列于表2中。由表2可以看出：基于單位脈沖的顯式時域法在模擬人造地震波樣本時，計算效率較高于Newmark-β直接積分法；基于本文方法得到的結果樣本方差與目標方差的符合程度較高，計算結果的最大相對誤差均在容許誤差范圍內。以上說明基于單位脈沖的顯式時域法在求解K-T模型模擬人造地震波時具有良好的計算效率和計算精度。事實上，在求解其他平穩模型時，如Clough-Penzien模型等，采用基于單位脈沖的顯式時域法也能達到同樣的目的，限于篇幅，本文只給出了基于K-T模型模擬人造地震波的結果。

圖2 單條地震波時程曲線

圖3 樣本方差與目標方差的對比

Newmark-β直接積分法和基于Newmark-β法的顯式時域法的結果與目標方差的偏差較大，主要原因在于Newmark-β法屬于逐步積分法，當取時間步長為0.02 s時，輸入的白噪聲并非理想白噪聲，而是限帶白噪聲，其帶寬為0～50 Hz。

表2 各種方法的計算時間和最大相對誤差Table 2 Calculation time and maximum relative error of various methods

3.2 與規范反應譜相符的人造地震波的模擬

我國現行的抗震規范規定，采用時程分析法時使用的地震波的平均反應譜曲線應與規范反應譜曲線在統計意義上相符[29]。對3.1節中利用基于單位脈沖的顯式時域法模擬的人造地震波樣本進行大量抽樣，獲得平均反應譜,將其與規范反應譜比較，如圖4所示。其中，α為地震影響系數，T為周期。由圖可見，根據單一譜強度因子S0模擬的人造地震波樣本與規范反應譜在形狀上是相似的，但在多數周期點上兩者的反應譜數值誤差較大。這說明地震動隨機模型無法根據單一譜強度因子S0過濾白噪聲，模擬與規范反應譜在整個周期內都吻合的人造地震波。

圖4 兩種反應譜的比較

為模擬與現行規范反應譜相符的人造地震波，采用鄧軍等[30]在研究時程分析法中選擇輸入地震記錄時，對[0 s,6 s]頻段進行劃分控制的思想，考慮利用基于單位脈沖的顯式時域法能快速求解響應的特點，通過迭代，先確定當結構的基本周期為T時與規范反應譜相符的譜強度因子S0(T)的統計值，進而用于模擬大量與規范反應譜相符的人造地震波。其主要步驟如下：

(1)根據譜強度因子與地震烈度、場地類別有關的特點，確定S0,L和S0,R的初始值，致使[S0,L,S0,R]包含S0(T)所有的可能取值，這里取S0,L=10-4，S0,R=1，利用二分法取S0(T),即S0(T)=(S0,L+S0,R)/2；

(3)計算結構基本周期為T時的平均反應譜值與相應規范反應譜值的相對誤差E(α)；

(5)反復進行上述步驟(1)至(4)，直至S0,L-S0,R和E(α)的絕對值滿足迭代容許誤差ε為止；

(6)得到與規范反應譜相符的譜強度因子S0(T)；

(7)模擬大量與規范反應譜相符的人造地震波。

采用MATLAB編程實現上述步驟。為說明上述修正譜強度因子S0(T)方法的有效性，取[0 s,6 s]內有限個周期點，考慮8度多遇地震、設計地震第一組、Ⅱ類場地、體系阻尼比0.05，利用上述方法得到T點上平均反應譜值與規范反應譜值的對比情況，如圖5所示。由圖5可見，修正譜強度因子S0后獲得的平均反應譜與規范反應譜的吻合度非常高，其他抗震設計條件下的結果也類似，這說明了上述修正譜強度因子S0(T)方法的有效性。

為便于應用，表3中給出了采用K-T模型時，在7度、8度和9度多遇地震、設計地震第一組、Ⅱ類場地、體系阻尼比0.05的抗震設計條件下，按本文方法得到的與規范反應譜相符的譜強度因子S0(T)的參考取值。

圖5 修正譜強度因子后平均反應譜與規范反應譜的對比

表3 采用K-T模型時部分譜強度因子S0(T)的參考取值Table 3 Partial spectral intensity factor S0(T) reference value when using K-T model 10-4 m2/s3

4 數值算例

分別采用Newmark-β直接積分法、基于單位脈沖與振型疊加法的顯式時域法計算A點水平自由度的位移響應，兩種方法的計算時間分別為145.9 s和2.1 s。圖7給出了單條地震波作用下求得的位移響應，可以看出兩種方法的計算結果基本吻合。這說明基于單位脈沖與振型疊加法的顯式時域法在進行結構動力時程分析時，具有良好的計算精度和較高的計算效率。

圖6 平面框架結構模型

圖7 A點的位移響應時程

為驗證基于單位脈沖與振型疊加法的顯式時域法作為時程分析法時的合理性，采用該法作為顯式時域隨機模擬分析法，得到結構受隨機地震作用時的抗震隨機模擬分析結果，并與振型分解反應譜法的計算結果進行對比。層間剪力和質心位移的計算結果如圖8所示，層間剪力最大差異為16.54%，質心位移最大差異為13.52%?；趩挝幻}沖與振型疊加法的顯式時域法考慮了地震荷載的隨機性，計算結果具有統計意義，兩種方法的層間剪力最大差異滿足現行規范對采用時程法分析時的要求，說明基于單位脈沖與振型疊加法的顯式時域法可用于建筑結構的抗震計算。

圖8 層間剪力和質心位移

5 結論

本文基于單位脈沖解析響應推導了線性動力系統響應的顯式時域表達式，并將其應用于模擬指定強度和與現行規范反應譜相符的人造地震波，以及結構抗震隨機模擬分析，主要得到以下結論：

(1)無論選用何種地震動隨機模型，利用基于單位脈沖的顯式時域法能快速求解地震動隨機模型的地面運動方程的特點，通過迭代，可以高效地模擬大量指定強度的地震波和與現行規范反應譜相符的人造地震波；

(2)基于單位脈沖與振型疊加法的顯式時域法在進行大型結構抗震分析時，可以進一步實現降維計算，具有相當高的計算效率，非常適用于大型結構的隨機響應時程分析。