广义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种经典均匀与非均匀调制非平稳随机地震响应分析的精细积分精确格式具体解析解, 使用方便、 应用范围广、 具有代表性。