滯后阻尼體系地震反應的中心差分虛初始條件法

祁文睿，潘旦光,2，高永濤，付相球

(1. 北京科技大學土木與資源工程學院，北京 100083；2. 城市地下空間工程北京市重點實驗室，土木與資源工程學院，北京科技大學，北京 100083)

阻尼為系統的固有特征，對結構的動力反應有顯著影響。根據材料的耗能特性，可建立相應的阻尼模型[1]。對于土木工程而言，粘滯阻尼和滯后阻尼是兩個廣泛使用的模型。粘滯阻尼力與速度成正比，具有計算簡便的特點而成為廣泛使用的模型[2]，如Rayleigh 阻尼模型和速度相關型的阻尼器[3?5]。但是，粘滯阻尼的耗能和頻率相關。已有實驗表明，很多材料的耗能與頻率無關，如型鋼混凝土梁[6]、土[7?8]和粘彈性夾層梁和板[9?10]等，此時采用阻尼力與位移成正比相位差π/2 的滯后阻尼模型更符合實驗結果。滯后阻尼將導致復剛度運動方程，而常采用頻域方法進行求解。頻域計算方法以Fourier 變換為基礎，理論上僅適用于線彈性體系，對于非線性系統，常采用等效線性化進行計算[11?13]。為進行真非線性計算，需要在時域中進行直接積分法進行計算。

為在時域中進行復本構運動方程的求解，朱鏡清[14?15]、何鐘怡[16]等根據對偶原則，建立與實荷載對應的對偶虛荷載，完善了滯后阻尼體系的輸入理論。但是，采用時域直接積分法計算滯后阻尼體系動力反應，即使是無條件穩定的直接積分方法也易于出現發散的結果[17?19]。微分方程數值解不穩定的原因有兩種：一種是數值方法的不穩定；另一種是由于方程本身具有發散解。復阻尼運動方程的逐步積分法的發散解是由后一種原因造成的[17,20]。為使復阻尼運動方程的逐步積分法計算結果穩定，一種常用的方法是將滯后阻尼等效為近似的粘滯阻尼[21?22]，但計算結果的誤差較大。

在直接積分法得到滯后阻尼體系穩定解方面，孫攀旭等分別基于激勵插值方法[23]和滯變阻尼時域理論[24]，建立了穩定的直接積分計算方法。事實上，滯后阻尼體系的特征值為互為相反數的復數特征對[25]，必然有一個特征值的實部是正的，由此導致強迫荷載下補解中的一項是沒有物理意義的發散項，在解析解時人工刪除發散項而使計算結果穩定。而直接積分法的計算結果是包含補解的，且滯后阻尼體系直接積分法是收斂到包含發散項補解[17]而導致計算結果不穩定。Pan 等[26]首先提出虛初始條件的概念，然后基于無條件穩定的Newmark 法建立使滯后阻尼直接積分解不出現發散項的計算方法。本文進一步建立了恒載下的虛初始條件，以及基于有條件穩定的中心差分法建立滯后阻尼的中心差分虛初始條件法。在此基礎上，通過算例驗證驗證所提方法的穩定性、計算精度和計算效率。

1 滯后阻尼體系的虛初始條件

1.1 自由振動

1.2 簡諧荷載作用

在簡諧荷載Aeiθt作用下，單自由度體系的運動方程為：

1.3 恒載作用

1.4 地震作用

若已知地面運動加速度為a(t)，將a(t)采用Fourier 級數展開：

2 虛初始條件的中心差分計算方法

若已知tn時刻及以前時刻的反應已知的情況下，根據中心差分法，求解tn+1時刻的位移方程為：

表1 任意荷載下中心差分虛初始條件法計算步驟Table 1 Calculation procedure of central differential virtual initial condition method under arbitrary load

3 算例

為驗證本文算法的有效性，下面以文獻[17]分析過的算例進行本文算法的驗證。已知體系的滯后阻尼系數η=0.1，對無阻尼自振頻率f=0.1 Hz、1 Hz、10 Hz 的三個體系進行地震反應計算。以表2中的3 條地震波分別作為體系的地震輸入。輸入地震波的加速度時程如圖1 所示。在進行中心差分法計算時，根據中心差分法的穩定性條件和輸入地震波的采樣時間間隔，計算時間步長 ?t取為：

式中：Tn為體系的自振周期；Δts為地震波加速度時程的采樣時間間隔。

表2 地震波主要參數Table 2 The main parameters of ground motions

圖1 輸入地震波的加速度時程Fig. 1 Acceleration time histories of input seismic waves

作為對比，對地震作用形成式(31)的方程后，計算各項簡諧荷載的解析解并求和，所得結果作為解析解。當u(0)=u0,u˙(0)=v0時，位移、速度和加速度的解析解為：

利用式(46)進行Fourier 逆變換，可得位移、速度和加速度的頻域解為：

對比式(44)和式(47)可知，解析解包括由初始條件引起的自由振動，由荷載引起的伴生自由振動和荷載所引起的純強迫振動。而頻域解僅包含荷載所引起的純強迫振動。

當u0=v0=0，f=10 Hz 時，三條地震波作用下體系的位移反應如圖2 所示，f=10 Hz 時天津波作用下的速度和加速度反應如圖3 所示，f=1 Hz 和0.1 Hz 時天津波作用下的位移反應如圖4 所示。

圖3 天津波作用下的速度和加速度反應(f=10 Hz, u0=v0=0)Fig. 3 Velocity and acceleration response under Tianjin wave excitation

由圖2～圖4 可知，對于不同的地震輸入和自振頻率體系，中心差分虛初始條件法計算的位移、速度和加速度都和解析解基本一致，計算結果穩定而不存在臨界自振周期的問題[17]。而對于頻域解，當f=10 Hz 時與解析解幾乎重合；當f=0.1 Hz 時，頻域解與解析解存在明顯差別。這是因為頻域解得到的是體系的穩態反應，而解析解是包含穩態和伴生自由振動的瞬態反應，自振頻率越高，瞬態反應衰減越快，因此，對于自振頻率較高的體系，瞬態反應對體系的總反應影響很小而可直接采用穩態解進行描述，但是，對于自振頻率低的體系，瞬態反應衰減需要較長的時間，當自振周期和地震波持時為同一量級時，則瞬態反應在整個地震反應過程中都有顯著影響，此時采用僅包含穩態解的頻域分析方法將產生顯著誤差。本文所提的中心差分虛初始條件法，計算結果包含瞬態反應和穩態反應，對于不同自振周期體系的反應都和解析解吻合的很好。

圖2 三條地震作用下位移反應(f=10 Hz, u0=v0=0)Fig. 2 Displacement responses under three seismic excitations(f=10 Hz, u0=v0=0)

圖4 天津波作用下不同自振頻率體系的位移反應(u0=v0=0)Fig. 4 Displacement response for various natural vibration frequencies under Tianjin wave excitation

當體系初始條件非零時，在天津波作用下不同自振頻率的位移反應如圖5 所示，f=10 Hz 體系的速度與加速度反應如圖6 所示。初始條件非零情況與零初始條件下的計算結果對比可知，由于初始位移和初始速度瞬態振動的進一步影響，導致頻域解在初始階段誤差增大。對于自振頻率較高的10 Hz 體系，由于瞬態振動很快消失，此時忽略初始條件對體系總反應的影響較小。但是對于自振周期與輸入地震波持時為同一量級的0.1 Hz 體系，初始位移和初始速度引起的瞬態振動將對總反應產生顯著影響，此時忽略初始條件將引起顯著誤差。

圖5 天津波作用下不同自振頻率體系的位移反應Fig. 5 Displacement response of various natural vibration frequencies under Tianjin wave excitation

為定量研究不同方法的精確性，不同方法相對解析解的峰值相對誤差為：

式中： |r?(t)|max表示精確解的峰值； |r(t)|max表示近似解的峰值。

圖6 天津波作用下體系的速度和加速度反應(f=10 Hz)Fig. 6 Velocity and acceleration response under Tianjin wave excitation

零初始條件下本文方法及頻域解的峰值相對誤差如表3 所示。由表3 可知，當體系的自振頻率較高時，頻域解的計算精度高；而當體系的自振頻率較低時，頻域解將產生顯著誤差；而中心差分虛初始條件法對于不同自振周期體系反應的峰值相對誤差都小于5%，顯示了良好的精度。

表4 為u0=v0=0，f=1 Hz 時，天津波作用下三種方法位移反應的計算時間。解析解為輸入地震波的Fourier 變換時間和式(44a)各項直接求和的計算時間，中心差分虛初始條件法的計算時間包括輸入地震波的Fourier 變換時間、計算虛初始條件的時間和每個時刻位移、速度和加速度的遞推計算的時間，頻域法則為Fourier 變換和Fourier 逆變換的總時間。在進行解析解計算時，由于自由振動和伴生自由振動的影響，無法采用快速Fourier變換進行計算而采用逐項求和，因此，計算時間最長。頻域采用快速Fourier 變換，計算時間最短。中心差分虛初始條件法通過前一步的計算結果遞推后一步的計算結果，計算效率僅次于頻域方法，且計算時間顯著小于解析解。結合表3 和表4 的數據分析結果表明，中心差分虛初始條件法兼顧了計算精度和計算效率。

表3 本文方法與頻域解的峰值相對誤差Table 3 Peak relative errors of the proposed method and frequency domain solution

表4 不同算法計算時間比較 /sTable 4 Comparison of calculation time of various methods

4 結論

針對滯后阻尼體系直接積分法收斂到包含發散項而不穩定的問題，提出了中心差分虛初始條件計算方法。由理論分析和數值計算結果可知：

(1)對于不同峰值加速度、卓越頻率的地震輸入，及不同自振頻率體系，本文所提的中心差分虛初始條件法對位移、速度和加速度均能得到穩定的計算結果，不存在臨界自振周期的問題。實初始條件是可觀測和測量部分，虛初始條件是伴隨實初始條件而存在的，可無需人為干涉而使滯后阻尼體系計算結果穩定。由此進行直接積分法計算，即使有條件穩定的中心差分法，依然可以得到穩定的計算結果。

(2)對于自振頻率較高的體系，瞬態反應對體系的總反應影響很小而可直接采用穩態解進行描述，而對于自振頻率低的體系，如自振周期和地震波持時為同一量級時，瞬態反應在整個地震反應過程中都有顯著影響，此時采用僅包含穩態解的頻域分析方法將產生顯著誤差。

(3)中心差分虛初始條件法的計算結果包括瞬態反應和穩態反應，計算誤差和體系的自振頻率無關，可適用于不同的自振頻率體系。

(4)中心差分虛初始條件法計算的峰值相對誤差小于5%。同時，計算時間顯著小于解析解，因此，這種方法兼顧了計算精度和計算效率。