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

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

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

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

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

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

1.1 自由振動(dòng)

1.2 簡(jiǎn)諧荷載作用

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

1.3 恒載作用

1.4 地震作用

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

2 虛初始條件的中心差分計(jì)算方法

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

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

3 算例

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

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

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

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

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

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

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

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

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

由圖2～圖4 可知，對(duì)于不同的地震輸入和自振頻率體系，中心差分虛初始條件法計(jì)算的位移、速度和加速度都和解析解基本一致，計(jì)算結(jié)果穩(wěn)定而不存在臨界自振周期的問(wèn)題[17]。而對(duì)于頻域解，當(dāng)f=10 Hz 時(shí)與解析解幾乎重合；當(dāng)f=0.1 Hz 時(shí)，頻域解與解析解存在明顯差別。這是因?yàn)轭l域解得到的是體系的穩(wěn)態(tài)反應(yīng)，而解析解是包含穩(wěn)態(tài)和伴生自由振動(dòng)的瞬態(tài)反應(yīng)，自振頻率越高，瞬態(tài)反應(yīng)衰減越快，因此，對(duì)于自振頻率較高的體系，瞬態(tài)反應(yīng)對(duì)體系的總反應(yīng)影響很小而可直接采用穩(wěn)態(tài)解進(jìn)行描述，但是，對(duì)于自振頻率低的體系，瞬態(tài)反應(yīng)衰減需要較長(zhǎng)的時(shí)間，當(dāng)自振周期和地震波持時(shí)為同一量級(jí)時(shí)，則瞬態(tài)反應(yīng)在整個(gè)地震反應(yīng)過(guò)程中都有顯著影響，此時(shí)采用僅包含穩(wěn)態(tài)解的頻域分析方法將產(chǎn)生顯著誤差。本文所提的中心差分虛初始條件法，計(jì)算結(jié)果包含瞬態(tài)反應(yīng)和穩(wěn)態(tài)反應(yīng)，對(duì)于不同自振周期體系的反應(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)體系初始條件非零時(shí)，在天津波作用下不同自振頻率的位移反應(yīng)如圖5 所示，f=10 Hz 體系的速度與加速度反應(yīng)如圖6 所示。初始條件非零情況與零初始條件下的計(jì)算結(jié)果對(duì)比可知，由于初始位移和初始速度瞬態(tài)振動(dòng)的進(jìn)一步影響，導(dǎo)致頻域解在初始階段誤差增大。對(duì)于自振頻率較高的10 Hz 體系，由于瞬態(tài)振動(dòng)很快消失，此時(shí)忽略初始條件對(duì)體系總反應(yīng)的影響較小。但是對(duì)于自振周期與輸入地震波持時(shí)為同一量級(jí)的0.1 Hz 體系，初始位移和初始速度引起的瞬態(tài)振動(dòng)將對(duì)總反應(yīng)產(chǎn)生顯著影響，此時(shí)忽略初始條件將引起顯著誤差。

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

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

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

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

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

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

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

表4 不同算法計(jì)算時(shí)間比較 /sTable 4 Comparison of calculation time of various methods

4 結(jié)論

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

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

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

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

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