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

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

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

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

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

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

1.1 自由振動

1.2 簡諧荷載作用

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

1.3 恒載作用

1.4 地震作用

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

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

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

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

3 算例

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

表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 結(jié)論

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

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

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

(3)中心差分虛初始條件法的計算結(jié)果包括瞬態(tài)反應(yīng)和穩(wěn)態(tài)反應(yīng)，計算誤差和體系的自振頻率無關(guān)，可適用于不同的自振頻率體系。

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