廣義Maxwell阻尼耗能系統非均勻響應精細積分精確格式

李創第, 王博文, 昌明靜

(廣西科技大學 土木建筑工程學院, 廣西 柳州 545006)

0 引 言

在結構抗風抗震中, 粘彈性阻尼器已得到廣泛應用。以Maxwell模型為基礎[1-3]可進一步得到更加精確的廣義Maxwell阻尼器模型[4-5], 實際工程中可以用精確的廣義Maxwell模型精準地表示線性流體和線性固體粘彈性阻尼器, 且對于分數導數模型[6-7]和Shen and Soong模型[8]來說, 采用廣義Maxwell模型得到的實驗數據精度更高, 效果更好。

在減震結構系統中, 阻尼器一般與支撐串聯安裝[9-10], 我國《建筑抗震設計規范》[11]通過控制支撐的最小剛度, 使得串聯安裝系統達到或接近純阻尼器的效果, 所以結構系統響應分析要考慮支撐的影響[12]。由于地震發生會首先引起支撐、 阻尼器等結構保護系統的破壞, 進而導致結構體系的損傷甚至毀滅, 目前相關規范明確要求耗能減震系統構件在結構設計基準期內應具備足夠的變形、 耗能能力和良好的抗震動力可靠度[11-12], 故結構及結構保護系統響應方法的建立, 對于結構及阻尼器構件抗震動力可靠度乃至抗震設計方法的研究具有非常重要的意義。地震動的非平穩隨機特性具有頻率非平穩和強度非平穩的特點, 因此, 對于非平穩隨機地震響應的分析較平穩隨機地震響應的分析更加符合工程實際意義[13]。 目前， 均勻調制隨機激勵的研究使大多數地震的非平穩隨機響應分析局限于此, 因此, 非均勻非平穩激勵的研究正日益受到廣大科研人員的高度重視[14-16]。林家浩等提出了高效的虛擬激勵法[17-18], 簡諧、 多項式簡諧、 指數簡諧型精細積分格式[17-20], 已應用于無阻尼器結構的均勻非平穩隨機響應高效分析[21], 但僅對于特定形式的激勵效率較高, 對于均勻與非均勻精細積分精確一般格式尚未建立, 且關于設置支撐的粘彈性阻尼耗能結構系統基于該方法的非平穩響應分析尚未建立。本文提出的精細積分精確一般格式可以退化為林家浩等提出的簡諧、 多項式簡諧、 指數衰減型簡諧3種精細積分格式[17], 應用更廣泛， 更具有代表性。

精細積分法最大的特點是不受計算步長的影響, 任意步長計算效果都是精確的, 從而使計算效率大大提高[17-21]。Newmark法是一種單步積分計算方法[17, 19, 21], 每個時間節點計算結果都是近似解, 隨著步長的增大精度下降很快, 步長增大到一定范圍之后計算結果發散, 結果出錯, 所以計算步長只能取一定限值， 否則將導致計算效率降低。

本文方法得出的結果是精細積分一般格式的精確解析解以及8種具體經典調制精細積分格式的精確解析解, 應用于設置支撐的廣義Maxwell阻尼耗能結構系統均勻與非均勻隨機地震響應分析[22]、 設置支撐的Maxwell阻尼耗能結構非平穩地震響應分析[23]中, 在受演變隨機激勵結構響應的精細逐步積分法[21]中, 使用的HPD-L以及HPD-S精細積分格式得到的結果是近似解, 并且不具有代表性。本文方法建立精細積分一般精確格式以及8種經典調幅精細積分精確格式的具體解析解, 使用方便、 應用范圍廣、 具有代表性, 針對設置支撐的廣義Maxwell阻尼耗能系統, 結合虛擬激勵法, 構建了均勻與非均勻非平穩地震響應的數值分析方法。

1 結構運動方程

1.1 廣義Maxwell阻尼器模型的本構關系

設廣義Maxwell阻尼器受力為PQ(t), 如圖1所示, 阻尼器相對于地面位移為xQ, 其中標準Maxwell阻尼器單元的個數為n, 阻尼器的平衡剛度為k0, 阻尼器第i個阻尼單元的剛度和阻尼分別為ki和ci。那么阻尼器受力可表示為[10]

圖1 廣義Maxwell阻尼器模型

(1)

阻尼器第i個阻尼單元的阻尼力、 松弛時間倒數、 松弛函數分別為Pi(t)、μi、hQi(t)，i=1,2,…,n。

(2)

μi=ki/ci;hQi(t)=kie-μit。

(3)

由式(2)和式(3), 可得

Pi(t)+μiPi(t)=kixQ(t)。

(4)

1.2 支撐與阻尼器的關系

工程實際中阻尼器一般與支撐串聯安裝, 以產生更好的減震效果, 如圖2所示, 支撐剛度為kb, 支撐相對于地面位移為xb, 結構相對于地面位移為x。那么支撐位移xb與結構位移x、 阻尼器位移xQ之間可表示為

圖2 設置支撐的廣義Maxwell阻尼器模型

xb=x-xQ。

(5)

設支撐的受力為Pb(t), 由于串聯安裝， 故支撐受力與阻尼器受力PQ(t)相同, 即

PQ(t)=Pb(t)=kbxb。

(6)

1.3 結構運動方程的建立

圖3 結構模型

(7)

將式(5)代入式(6), 同時考慮式(1)、 式(7), 可以寫為

(8)

(9)

將式(9)分別代入式(8)和(4), 最終可得

(10)

(11)

1.4 擴階方程的建立

令

(12)

式(10)～式(12)以擴階的形式表示為

(13)

式中:

(14)

(15)

(16)

2 非均勻響應分析的虛擬激勵法

2.1 非均勻非平穩地震激勵模型

(17)

g(ω,t)=ε0(ω)eα0(ω)t+ε1(ω)teα1(ω)t+

ε2(ω)t2eα2(ω)t。

(18)

2.2 非均勻響應分析的虛擬激勵法

(19)

(20)

其中:

(21)

對于多自由度耗能結構體系同樣也可化為上式, 同樣可得多自由度耗能結構體系非平穩響應解析解。式(20)的通解為齊次解與特解之和, 即

Z(ω,t)=T(τ)(Z(ω,tk)-Zp(ω,tk))+Zp(ω,t),

(22)

式中: 積分步長t∈[tk,tk+1],τ=t-tk。 關于指數矩陣T(τ)的精細計算, 詳見文獻[18]。問題歸結為求特解Zp(ω,t)及精細地計算T(τ)。

3 非均勻精細積分一般精確格式

由式(18)和式(21), 激勵荷載在每一積分步長t∈[tk,tk+1]內可表示為

f0(ω,t)=-B-1f1(ω)g(ω,t)eiωt

=(r0eα0(ω)t+r1teα1(ω)t+r2t2eα2(ω)t)×

(asinωt+bcosωt),

(23)

式中:r0=-B-1f1(ω)ε0(ω);r1=-B-1f1(ω)ε1(ω);r2=-B-1f1(ω)ε2(ω);a=i;b=1。

由式(23)可得方程的特解Zp(ω,t)為

Zp(ω,t)=(a0eα0(ω)t+a1teα1(ω)t+a2t2eα2(ω)t)sinωt+

(b0eα0(ω)t+b1teα1(ω)t+b2t2eα2(ω)t)cosωt,

(24)

式中:a2=((α2(ω)I-H)2+ω2I)-1×((α2(ω)I-H)ar2+ωbr2);b2=((α2(ω)I-H)2+ω2I)-1×((α2(ω)I-H)br2-ωar2);a1=((α1(ω)I-H)2+ω2I)-1×((α1(ω)I-H)(ar1-2eα2(ω)t-α1(ω)ta2)+ω(br1-2eα2(ω)t-α1(ω)tb2));b1=((α1(ω)I-H)2+ω2I)-1×((α1(ω)I-H)(br1-2eα2(ω)t-α1(ω)tb2)-ω(ar1-2eα2(ω)t-α1(ω)ta2));a0=((α0(ω)I-H)2+

ω2I)-1×((α0(ω)I-H)(ar0-eα1(ω)t-α0(ω)ta1)+ω(br0-eα1(ω)t-α0(ω)tb1));b0=((α0(ω)I-H)2+ω2I)-1×((α0(ω)I-H)(br0-eα1(ω)t-α0(ω)tb1)-ω(ar0-eα1(ω)t-α0(ω)ta1))。

由式(22)和式(24)非均勻精細積分一般精確格式可表示為

Z(ω,tk+1)=T(τ)(Z(ω,tk)-Zp(ω,tk))+

Zp(ω,tk+1)。

(25)

式(23)可以退化為林家浩提出的簡諧、 多項式簡諧、 指數衰減型簡諧3種精細積分格式[14]。由于篇幅所限, 僅介紹: ① 當激勵荷載中r1=r2=0,α0=α1=α2=0時, 可退化為簡諧荷載精細積分格式; ② 當激勵荷載中α0=α1=α2=0時, 可退化為多項式簡諧荷載精細積分格式; ③ 當激勵荷載中r1=r2=0時, 可退化為指數衰減型簡諧荷載精細積分格式。

Szz(ω,t)=z*(ω,t)z(ω,t)，

(26)

(27)

式中: “*”表示復共軛。綜上步驟, 設置支撐的廣義Maxwell阻尼耗能系統的非均勻非平穩響應方差均可得到。

4 8種經典調幅精細積分精確格式

由式(25)得, 精細積分精確格式可由特解求出, 為節省篇幅以下計算只給出特解。

4.1 Shinozuka-Sato型調制函數

g(t)=e-λ1t-e-λ2t,

(28)

式中:λ1、λ2為已知常數。

f0(ω,t)=(r0,0eα0,0(ω)t+r0,1eα0,1(ω)t)×

(asinωt+bcosωt),

(29)

式中:r0,0=-B-1f1ε0,0,ε0,0=1,α0,0=-λ1;r0,1=-B-1f1ε0,1,ε0,1=-1,α0,1=-λ2;a=i;b=1。

可求得特解

Zp(ω,t)=(a0,0eα0,0t+a0,1eα0,1t)sinωt+

(b0,0eα0,0t+b0,1eα0,1t)cosωt,

(30)

式中:a0,1=((α0,1I-H)2+ω2I)-1×((α0,1I-H)ar0,1+ωbr0,1);b0,1=((α0,1I-H)2+ω2I)-1×((α0,1I-H)br0,1-ωar0,1);a0,0=((α0,0I-H)2+

ω2I)-1×((α0,0I-H)ar0,0+ωbr0,0);b0,0=((α0,0I-H)2+ω2I)-1×((α0,0I-H)br0,0-ωar0,0)。

4.2 Hsu-Bernard型調制函數

g(t)=εte-λt,

(31)

式中:ε=λe,λ為已知常數。

f0(ω,t)=r1teα1t(asinωt+bcosωt),

(32)

其中:r1=-B-1f1ε1;ε1=λe;α1=-λ;a=i;b=1。

可求得特解

Zp(ω,t)=(a0+a1t)eα1tsinωt+(b0+b1t)eα1tcosωt,

(33)

式中:a1=((α1I-H)2+ω2I)-1((α1I-H)ar1+ωbr1);b1=((α1I-H)2+ω2I)-1((α1I-H)br1-ωar1);a0=((α1I-H)2+ω2I)-1×((α1I-H)(-a1)+ω(-b1));b0=((α1I-H)2+ω2I)-1×((α1I-H)(-b1)-ω(-a1))。

4.3 Goto-Toki型調制函數

(34)

式中:A0、tp為已知常數。

f0(ω,t)=r1teα1t(asinωt+bcosωt),

(35)

可求得特解

Zp(ω,t)=(a0+a1t)eα1tsinωt+(b0+b1t)eα1tcosωt,

(36)

式中:a1=((α1I-H)2+ω2I)-1((α1I-H)ar1+ωbr1);b1=((α1I-H)2+ω2I)-1((α1I-H)br1-ωar1);a0=((α1I-H)2+ω2I)-1×((α1I-H)(-a1)+ω(-b1));b0=((α1I-H)2+ω2I)-1×((α1I-H)(-b1)-ω(-a1))。

4.4 Iyengar型調制函數

g(t)=(c+dt)e-λt,

(37)

式中:c、d、λ為已知常數。

f0(ω,t)=(r0eα0t+r1teα1t)(asinωt+bcosωt),

(38)

式中:r0=-B-1f1ε0,ε0=c,α0=-λ;r1=-B-1f1ε1,ε1=d,α1=-λ;a=i;b=1。

可求得特解

Zp(ω,t)=(a0eα0t+a1teα1t)sinωt+(b0eα0t+

b1teα1t)cosωt,

(39)

式中:a1=((α1I-H)2+ω2I)-1((α1I-H)ar1+ωbr1);b1=((α1I-H)2+ω2I)-1((α1I-H)br1-ωar1);a0=((α0I-H)2+ω2I)-1×((α0I-H)(ar0-a1)+ω(br0-b1));b0=((α0I-H)2+ω2I)-1×((α0I-H)(br0-b1)-ω(ar0-a1))。

4.5 分段型調制函數

(40)

式中:A0、c、t1、t2為已知常數。

當0≤t≤t1時,

f0(ω,t)=r2t2eα2t(asinωt+bcosωt),

(41)

可求得特解

Zp(ω,t)=(a0+a1t+a2t2)sinωt+(b0+b1t+

b2t2)cosωt,

(42)

式中:a2=((α2I-H)2+ω2I)-1((α2I-H)ar2+ωbr2);b2=((α2I-H)2+ω2I)-1((α2I-H)br2-ωar2);a1=((α2I-H)2+ω2I)-1×((α2I-H)(-2a2)+ω(-2b2));b1=((α2I-H)2+ω2I)-1×((α2I-H)(-2b2)-ω(-2a2));a0=((α2I-

H)2+ω2I)-1×((α2I-H)(-a1)+ω(-b1));b0=

((α2I-H)2+ω2I)-1×((α2I-H)(-b1)-ω(-a1))。

當t1≤t≤t2時,

f0(ω,t)=r0eα0t(asinωt+bcosωt),

(43)

式中:r0=-B-1f1ε0;ε0=A0;α0=0;a=i;b=1。

可求得特解

Zp(ω,t)=a0sinωt+b0cosωt,

(44)

式中:a0=(H2+ω2I)-1(-Har0+ωbr0),b0=(H2+ω2I)-1(-Hbr0-ωar0)。

當t≥t2時,

f0(ω,t)=(r0eα0t+r1teα1t)(asinωt+bcosωt),

(45)

式中:r0=-B-1f1ε0,ε0=-A0e-ct2,α0=0;r1=-B-1f1ε1,ε1=A0e-c,α0=0;a=i,b=1。

可求得特解

Zp(ω,t)=(a0+a1t)sinωt+(b0+b1t)cosωt,

(46)

式中:a1=(H2+ω2I)-1(-Har1+ωbr1);b1=(H2+ω2I)-1((-Hbr1)-ωar1);a0=(H2+ω2I)-1((-H(ar0-a1))+ω(br0-b1));b0=(H2+ω2I)-1((-H(br0-b1))-ω(ar0-a1))。

4.6 余弦型調制函數

g(t)=c+dcosθt,

(47)

式中:c、d、θ為已知常數;c≥d。

f0(ω,t)=(r0,0eα0,0t+r0,1eα0,1t+r0,2eα0,2t)×

(asinωt+bcosωt),

(48)

式中:r0,0=-B-1f1ε0,0,ε0,0=c,α0,0=0;r0,1=-B-1f1ε0,1,ε0,1=d/2,α0,1=iθ;r0,1=-B-1f1ε0,2,ε0,2=d/2,α0,2=-iθ;a=i;b=1。

可求得特解

Zp(ω,t)=(a0,0+a0,1eα0,1t+a0,2eα0,2t)sinωt+

(b0,0+b0,1eα0,1t+b0,2eα0,2t)cosωt,

(49)

式中:a0,2=((α0,2I-H)2+ω2I)-1((α0,2I-H)ar0,2+ωbr0,2);b0,2=((α0,2I-H)2+ω2I)-1((α0,2I-H)br0,2-ωar0,2);a0,1=((α0,1I-H)2+ω2I)-1((α0,1I-H)ar0,1+ωbr0,1);b0,1=((α0,1I-H)2+ω2I)-1((α0,1I-H)br0,1-ωar0,1);a0,0=(H2+ω2I)-1((-Har0,0)+ωbr0,0);b0,0=(H2+ω2I)-1((-Hbr0,0)-ωar0,0)。

4.7 正弦型調制函數

g(t)=c+dsinθt,

(50)

式中:c、d、θ為已知常數;c≥d。

f0(ω,t)=(r0,0eα0,0t+r0,1eα0,1t+r0,2eα0,2t)×

(asinωt+bcosωt),

(51)

式中:r0,0=-B-1f1ε0,0,ε0,0=c,α0,0=0;r0,1=-B-1f1ε0,1,ε0,1=d/(2i),α0,1=iθ;r0,1=-B-1f1ε0,2,ε0,2=-d/(2i),α0,2=-iθ;a=i,b=1。

可求得特解

Zp(ω,t)=(a0,0+a0,1eα0,1t+a0,2eα0,2t)sinωt+

(b0,0+b0,1eα0,1t+b0,2eα0,2t)cosωt,

(52)

式中:a0,2=((α0,2I-H)2+ω2I)-1((α0,2I-H)ar0,2+ωbr0,2);b0,2=((α0,2I-H)2+ω2I)-1(α0,2I-H)br0,2-ωar0,2);a0,1=((α0,1I-H)2+ω2I)-1((α0,1I-H)ar0,1+ωbr0,1);b0,1=((α0,1I-H)2+ω2I)-1((α0,1I-H)br0,1-ωar0,1);a0,0=(H2+ω2I)-1((-Har0,0)+ωbr0,0);b0,0=(H2+ω2I)-1((-Hbr0,0)-ωar0,0)。

4.8 Spanos-Solomos型非均勻調制函數

g(ω,t)=ε(ω)te-λ(ω)t,

(53)

式中:ε(ω)、λ(ω)表示以ω為自變量的函數。

f0(ω,t)=r1teα1t(asinωt+bcosωt),

(54)

式中:r1=-B-1f1ε1,ε1=ε(ω);α1=-λ(ω);a=i,b=1。

可求得特解

Zp(ω,t)=(a0+a1t)eα1tsinωt+

(b0+b1t)eα1tcosωt,

(55)

式中:a1=((α1I-H)2+ω2I)-1((α1I-H)ar1+ωbr1);b1=((α1I-H)2+ω2I)-1((α1I-H)br1-ωar1);a0=((α1I-H)2+ω2I)-1×((α1I-H)(-a1)+ω(-b1));b0=((α1I-H)2+ω2I)-1×((α1I-H)(-b1)-ω(-a1))。

5 算 例

如圖4所示, 單自由度設置支撐的五參數Maxwell阻尼器減震系統, 其結構的基本參數為: 質量m=38 600 kg, 剛度k=146.01×105N/m, 阻尼比s0取0.04。Maxwell阻尼器的基本參數為: 平衡剛度k0=0.36×105N/m, 支撐剛度kb=rbk,rb分別為0.5、 1.5、 3、 10、 ∞, Maxwell阻尼器兩分支單元的剛度和阻尼分別為k1=42.08×105N/m,c1=0.83×105N·s/m;k2=6.87×105N/m,c2=2.15×105N·s/m。

圖4 結構計算簡圖

調幅函數分別取為Shinozuka-Sato型[24-25]均勻調幅和Spanos-Solomos型[26]非均勻調幅, 計算參數分別取為

g(ω,t)=g1(ω,t)=e-0.6t-e-t,

取不同工況的支撐剛度, 利用本文均勻與非均勻調制非平穩激勵下, 其中兩種精細積分格式的具體精確解析解得到計算結果。在Shinozuka-Sato和Spanos-Solomos型均勻調幅非平穩地震作用下阻尼耗能結構系統(結構、 阻尼器)的響應分別如圖5、 圖6所示。可見支撐剛度取值不同， 對耗能系統響應有較大的影響, 而對產生響應最大值的時刻影響較小, 隨著支撐剛度取值增大, 阻尼器受力響應方差也隨之增大, 但是結構的位移、 速度響應方差反而減小。為保證結構獲得很好的耗能作用, 在支撐剛度kb≥10k情況下, 可以按kb=∞近似計算; 在kb較小情況下, 可以按kb的實際剛度進行計算。 在均勻與非均勻非平穩激勵下, 耗能系統的響應具有相似波動趨勢的特點, 表現出明顯的非平穩隨機特性, 符合工程實際。

圖5 Shinozuka-Sato型均勻調幅結構位移、 速度與阻尼器受力響應方差

圖6 Spanos-Solomos型非均勻調幅結構位移、 速度與阻尼器受力響應方差

6 結 論

為建立粘彈性耗能結構及其保護系統(結構、 阻尼器)的抗震分析與設計方法, 本文將虛擬激勵法引入粘彈性耗能阻尼系統, 提出了非均勻精細積分一般精確格式, 并成功高效地應用于設置支撐廣義Maxwell阻尼耗能系統的均勻與非均勻非平穩地震響應分析中。

(1)通過算例, 驗證了該方法適應于設置支撐的廣義Maxwell阻尼系統的非均勻非平穩響應分析, 并且可直接應用于粘彈性阻尼耗能系統響應分析。支撐剛度對粘彈性耗能系統有重要影響, 在支撐剛度較耗能系統剛度很大的情況下, 支撐剛度對耗能系統響應的影響效果不再增加, 一般情況下, 應考慮有限支撐剛度對耗能系統響應的影響。

(2)本文非均勻精細積分一般精確格式, 不受計算步長的影響, 任意步長計算效果都是精確的。相同步長而言, 非均勻精細積分一般精確格式比傳統辦法(如Newmark法)精度和效率都有顯著提高。簡諧、 多項式簡諧、 指數簡諧型精細積分格式都可以用本文非均勻精細積分一般精確格式表示, 應用范圍更加廣泛。本文得到8種經典均勻與非均勻調制非平穩隨機地震響應分析的精細積分精確格式具體解析解, 使用方便、 應用范圍廣、 具有代表性。