阻尼是指体系振动过程中能量耗散的特性, 阻尼模型及参数选取对建筑结构、地基地震响应分析至关重要.由于具有形式简单、数学处理方便等优点, 瑞利阻尼在动力问题分析中得到广泛应用^[1-2].

瑞利阻尼模型最早是由Caughey^[3]提出的广义瑞利阻尼模型退化而来, 可以实现基于正交模态解耦的复杂动力方程，在早期研究中得到广泛应用.大量研究围绕瑞利阻尼系数的选取.Idriss等^[4]采用基频确定瑞利阻尼系数, Yoshida等^[5]定义一个结构敏感的频率范围来计算瑞利阻尼系数, 楼梦麟等^[6]利用实际地震波的分析结果, 开展回归分析, 得到一个阻尼系数转换频率取值的经验公式, 邹德高等^[7]提出同时考虑结构频率特性和地震动频谱特性来确定瑞利阻尼系数, 楼梦麟等^[8]建议可以通过比较土层基频与地震波主频来确定瑞利阻尼系数.丁海平等^[9-11]对现有的系数确定方法进行改进或优化.楼梦麟等^[12-15]对不同参数选定方法得到的计算结果进行比较, 指出不合理的瑞利阻尼系数会导致较大的计算误差.

尽管瑞利阻尼能够较好地表征结构振动过程中的能量耗散, 并在地表分析中得到较多的应用, 但模型的物理意义尚不甚明确, 将瑞利阻尼应用于动响应分析可能会隐含各种问题.Nielsen等^[16-17]指出以总位移表述运动方程时包含了与质量比例阻尼有关的外部阻尼力；Chen等^[18]指出Abaqus/Explicit分析中, 忽略刚度比例阻尼会产生高频振荡；Hall^[17]通过算例分析发现, 若瑞利阻尼取值无上限, 则会得到与实际不符的较大阻尼力, 特别是在非线性问题中, 如果采用初始刚度求解瑞利阻尼, 则会使计算结果偏不安全.

由于瑞利阻尼是目前地震响应分析中经常采用的阻尼模型之一, 本文试图从一个视角阐析瑞利阻尼模型的物理意义, 揭示模型参数对动力响应, 特别是应力响应的影响.

1 瑞利阻尼的物理意义

在体系振动过程中, 能量耗散机制十分复杂, 产生阻尼的因素众多, 但总体上可以归结为两大类.一类是由材料本身的黏滞性引起的, 表现为应力不仅与应变有关, 还与应变的速率有关, 本文称这一类阻尼为材料阻尼.另一类为体系振动过程中支座摩擦、空气阻力等外部环境引起的能量耗散, 本文统称为环境阻尼.体系振动过程中的能量耗散源于材料阻尼和环境阻尼共同作用.

目前, 已有大量的材料阻尼本构模型被提出, 由弹簧和黏壶并联得到的Kelvin模型^[19]是其中最简单的模型之一, 可以表述为

$ \mathit{\boldsymbol{\sigma = D\varepsilon + \eta \dot \varepsilon }}. $

(1)

式中：σ和ε分别为应力和应变矢量, ${\mathit{\boldsymbol{\dot \varepsilon }}}$为应变率矢量, D和η分别为弹性和黏性系数矩阵.

将式(1) 代入虚功方程

$ \int_\mathit{\Omega } {{{\left( {{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{\varepsilon }}} \right)}^{\rm{T}}}\mathit{\boldsymbol{\sigma }}{\rm{d}}V} = \int_\mathit{\Omega } {{{\left( {{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{u}}} \right)}^{\rm{T}}}\mathit{\boldsymbol{f}}{\rm{d}}V} + \int_{{\mathit{\Gamma }_{\rm{t}}}} {{{\left( {{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{u}}} \right)}^{\rm{T}}} \cdot \mathit{\boldsymbol{\bar t}}{\rm{d}}\mathit{\Gamma }} $

(2)

的左端, 利用位移插值u=N·d, 可得

$ \int_\mathit{\Omega } {{{\left( {{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{\varepsilon }}} \right)}^{\rm{T}}}\mathit{\boldsymbol{\sigma }}{\rm{d}}V} = {\left( {{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{d}}} \right)^{\rm{T}}}\left( {\mathit{\boldsymbol{K}} \cdot \mathit{\boldsymbol{d + }}{\mathit{\boldsymbol{C}}_{\rm{M}}} \cdot \mathit{\boldsymbol{\dot d}}} \right), $

(3)

$ \mathit{\boldsymbol{K}} = \int_\mathit{\Omega } {{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{DB}}{\rm{d}}V} , $

(4)

$ {\mathit{\boldsymbol{C}}_{\rm{M}}} = \int_\mathit{\Omega } {{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{\eta B}}{\rm{d}}V} . $

(5)

式中：f和t分别为给定的体积力和表面力矢量, d和${\mathit{\boldsymbol{\dot d}}}$分别为节点位移和速度矢量, N为形函数矩阵, B为应变矩阵, Ω为分析域, Γ_t为给定力的边界, K和C_M分别为刚度矩阵和由材料阻尼引起的阻尼矩阵.

当黏性系数矩阵η和弹性系数矩阵间满足η=βD时,

$ {\mathit{\boldsymbol{C}}_{\rm{M}}} = \beta \mathit{\boldsymbol{K}}. $

(6)

外部环境引起的阻尼作用通常与体系振动速度有关, 作用方式非常复杂.为了简化起见, 假定环境阻尼作用可以等效为与速度相关的体积力f_E, 并表示为

$ {\mathit{\boldsymbol{f}}_{\rm{E}}} = - \mu \left( \mathit{\boldsymbol{v}} \right) \cdot \mathit{\boldsymbol{v}}. $

(7)

将式(7) 代入式(2) 右端的第一项, 可以求得环境阻尼作用对应的等效节点力为

$ {\mathit{\boldsymbol{F}}_{\rm{E}}} = - {\mathit{\boldsymbol{C}}_{\rm{E}}}\mathit{\boldsymbol{\dot d}}, $

(8)

$ {\mathit{\boldsymbol{C}}_{\rm{E}}} = \int_V {{\mathit{\boldsymbol{N}}^{\rm{T}}}\mu \left( \mathit{\boldsymbol{v}} \right)\mathit{\boldsymbol{N}}{\rm{d}}V} . $

(9)

式中：C_E为环境阻尼引起的阻尼矩阵.当μ(v)为常数, 且与质量密度ρ成正比, 即μ(v)=αρ时,

$ {\mathit{\boldsymbol{C}}_{\rm{E}}} = \alpha \mathit{\boldsymbol{M}}, $

(10)

$ \mathit{\boldsymbol{M = }}\int_\mathit{\Omega } {{\mathit{\boldsymbol{N}}^{\rm{T}}}\rho \mathit{\boldsymbol{N}}{\rm{d}}V} . $

(11)

式中：M为一致质量阵.

将式(6)、(10) 叠加, 总阻尼阵C可以表示为

$ \mathit{\boldsymbol{C}} = \alpha \mathit{\boldsymbol{M}} + \beta \mathit{\boldsymbol{K}}. $

(12)

式中：C为瑞利阻尼阵, α和β为相应的瑞利阻尼系数.

综上分析可以看出, 瑞利阻尼的刚度相关部分反映了材料对动力响应的阻滞作用, 表现为应力不仅与应变有关, 而且与应变率相关；瑞利阻尼的质量相关部分反映了环境对结构体系振动的阻滞作用.瑞利阻尼模型是一种综合反映材料和环境对结构体系动力响应阻滞作用的最简单模型, 只需2个材料参数.

需要指出的是, 瑞利阻尼模型的提出仅基于离散的结构体系, 与材料的性质无关, 因此在有限元软件中计算应力时, 没有考虑应变率对应力的贡献, 导致计算应力与实际应力存在差异.

2 一维黏弹性振动问题

为了进一步验证瑞利阻尼模型物理意义的合理性, 揭示模型参数对动力响应, 特别是应力响应的影响, 分别采用解析和数值方法对比分析均质等截面杆件一维黏弹性动力响应问题, 给出该问题的描述和解析解.如图 1所示, 杆件长L=10 m, 截面面积A=1 m², 弹性模量E=25 MPa, 质量密度ρ=1 800 kg/m³, 采用瑞利阻尼.在初始时刻, 杆件处于静止状态, 初始速度和初始位移均为0.杆件左端(x=0 m)均匀分布幅值p=50 Pa的随时间t正弦变化的荷载f(t), 杆件右端(x=L)自由.

黏弹性杆件振动的控制方程和相应的定解条件如下：

$ E\frac{{{\partial ^2}u}}{{\partial {x^2}}} + \beta E\frac{{{\partial ^3}u}}{{\partial {x^2}\partial t}} = \rho \frac{{{\partial ^2}u}}{{\partial {t^2}}} + \alpha \rho \frac{{\partial u}}{{\partial t}};0 < x < L,t > 0. $

(13)

$ u = 0,\frac{{\partial u}}{{\partial t}} = 0;0 < x < L,t = 0. $

(14)

$ \frac{{\partial u}}{{\partial x}} + \beta \frac{{{\partial ^2}u}}{{\partial x\partial t}} = - \frac{{f\left( t \right)}}{{EA}};x = 0,t > 0. $

(15)

$ \frac{{\partial u}}{{\partial x}} + \beta \frac{{{\partial ^2}u}}{{\partial x\partial t}} = 0;x = L,t > 0. $

(16)

式中：u为沿x方向的位移.

求解上述定解问题, 可得解析解.限于篇幅, 省略推导过程, 仅给出解的最终形式：

$ u\left( {x,t} \right) = \bar u\left( {x,t} \right) + \sum\limits_{k = 0}^\infty {{T_k}\left( t \right)\cos \left( {{\gamma _k}x} \right)} . $

(17)

式中:

$ {\gamma _k} = k{\rm{\pi }}/L,k = 0,1,2, \cdots , $

(18)

$ {T_k}\left( t \right) = \frac{1}{{\rho \left( {{r_{k2}} - {r_{k1}}} \right)}}\left[ {{K_{k2}}\left( t \right) - {K_{k1}}\left( t \right)} \right]. $

(19)

其中, ${K_{ki}}\left( t \right) = {{\rm{e}}^{^rk{i^t}}}\int_0^t {{h_k}} \left( \tau \right){{\rm{e}}^{^{ - r}k{i^\tau }}}{\rm{d}}\tau ;i = 1,2.$

$ \begin{array}{l} {h_k}\left( \tau \right) = \frac{{2\rho }}{L}\int_0^L {\left[ {\frac{E}{\rho }\left( {\frac{{{\partial ^2}\bar u}}{{\partial {x^2}}} + \beta \frac{{{\partial ^3}\bar u}}{{\partial {x^2}\partial \tau }}} \right) - } \right.} \\ \;\;\;\;\;\;\;\;\;\;\;\left. {\left( {\frac{{{\partial ^2}\bar u}}{{\partial {\tau ^2}}} + \alpha \frac{{\partial \bar u}}{{\partial \tau }}} \right)} \right]\cos \left( {{\gamma _k}x} \right){\rm{d}}x.\\ {r_{k1,k2}} = \\ \;\;\; - \frac{1}{2}\left[ {\left( {\frac{E}{\rho }\beta \omega _k^2 + \alpha } \right) \pm \sqrt {{{\left( {\frac{E}{\rho }\beta \omega _k^2 + \alpha } \right)}^2} - 4\frac{E}{\rho }\omega _k^2} } \right].\\ \bar u\left( {x,t} \right) = \left\{ \begin{array}{l} \frac{{{{\left( {L - x} \right)}^2}}}{{2EAL}} \cdot f\left( t \right),\beta = 0;\\ \frac{{{{\left( {L - x} \right)}^2}}}{{2EAL}} \cdot \int_0^t {\frac{{f\left( \tau \right)}}{\beta }{{\rm{e}}^{\frac{\tau }{\beta }}}{\rm{d}}\tau \cdot {{\rm{e}}^{ - \frac{1}{\beta }}},\beta \ne 0} . \end{array} \right. \end{array} $

(20)

将位移解代入几何关系, 可得杆件应变解：

$ \varepsilon \left( {x,t} \right) = \frac{{\partial \bar u\left( {x,t} \right)}}{{\partial x}} - \sum\limits_{k = 0}^\infty {{\gamma _k}{T_k}\left( t \right)\sin \left( {{\gamma _k}x} \right)} . $

(21)

利用应力和应变关系, 可得杆件应力解：

$ \begin{array}{*{35}{l}} \sigma \left( x,t \right)=E\left[ \frac{\partial \bar{u}\left( x,t \right)}{\partial x}+\beta \frac{{{\partial }^{2}}\bar{u}\left( x,t \right)}{\partial x\partial t} \right]- \\ \ \ \ \ \sum\limits_{k=0}^{\infty }{E{{\gamma }_{k}}\sin \left( {{\gamma }_{k}}x \right)\left[ {{T}_{k}}\left( t \right)+\beta {{{\dot{T}}}_{k}}\left( t \right) \right]}. \\ \end{array} $

(22)

如果忽略应变率对应力的贡献, 可得如下形式的应力解：

$ {\sigma ^ * }\left( {x,t} \right) = E\left[ {\frac{{\partial \bar u\left( {x,t} \right)}}{{\partial x}} - \sum\limits_{k = 0}^\infty {{\gamma _k}{T_k}\left( t \right)\sin \left( {{\gamma _k}x} \right)} } \right]. $

(23)

上述一维黏弹性杆件振动问题可以采用有限单元法进行数值求解.本文采用ABAQUS软件进行数值分析, 选用四节点四边形平面单元, 单元纵向尺寸为0.25 m.采用HHT算法^[20]进行动力隐式时程积分, 时间步长为0.001 s.

为了对比分析α和β对动响应, 尤其是应力响应的影响, 采用Yoshida等^[5]提出的基于体系敏感频率范围的方法确定基准阻尼系数α和β.

$ \left. \begin{array}{l} \xi = \frac{{\bar \alpha }}{{2{\omega _a}}} + \frac{{\bar \beta {\omega _a}}}{2},\\ \xi = \frac{{\bar \alpha }}{{2{\omega _b}}} + \frac{{\bar \beta {\omega _b}}}{2}. \end{array} \right\} $

(24)

式中：ξ为敏感频率上、下限对应的阻尼比, 本文取土体常用值0.05^[21].Clough等^[22]建议ω_a取多自由度体系的基频, ω_b在对动力响应有显著贡献的高阶振型频率中选取.本文分别取ω_a、ω_b为第1、3阶非零固有频率.经分析可知, ω_a=1.19 rad/s、ω_b=18.51 rad/s.由式(24) 可得, α=0.1 s^-1, β=0.005 s.为了方便起见, 定义无量纲瑞利阻尼系数λ_α=α/α, λ_β=β/β.

3 计算结果及分析

分别对比分析α和β对动力响应的影响, 重点分析刚度相关阻尼项对应力计算值的影响.为了简单起见, 取激振频率ω=10π rad/s(5 Hz).

图 2、3分别给出λ_β=1和λ_α=1时杆件中点(x=5 m)处的应变时程曲线.可见：1) 无论是α还是β变化, 应变数值解ε_N给出了与解析解ε_A几乎完全相同的结果, 说明有限元给出了可靠的应变解；2) 随着α或β的增加, 应变振动幅值减少, 响应滞后, 比较而言, 系数β的影响更明显.

图 4给出当λ_β=1和λ_α=1时, 杆件中点(x=5 m)处稳态应变幅值|ε|_max随α和β的变化曲线.可见, 在基准瑞利阻尼系数附近, α的变化对动力响应的影响远小于β的变化, 说明材料阻尼是主要的结构能量耗散源.

图 5、6给出当λ_β=1时, 杆件左端(x=0 m)和中点(x=5 m)处应力时程曲线, 包括解析解σ_A、解析解σ_A^*和数值解σ_N3种结果.图 5中, λ_α=1和λ_α=100时的曲线几乎完全重合.可以看出, 当λ_β=1时, σ_A^*和σ_N的结果一致, 但与σ_A的结果略有不同；α对σ_A和σ_A^*差别的影响不大.图 7、8给出当λ_α=1时, 杆件x=0 m、x=5 m处的应力时程曲线.由图 7、8可以看出, σ_A^*和σ_N的结果一致, σ_A和σ_A^*间的差异随着β的增加而明显增大.尤其值得注意的是, x=0 m处的σ_A^*不能满足边界条件, 而σ_A满足边界条件.说明现有有限元软件在分析瑞利阻尼的动响应问题时所求得的应力从理论角度上是不精确的.

为了进一步验证α和β对动力响应影响的上述规律, 取ω分别为5π rad/s和15π rad/s进行类似分析.

为了方便起见, 仅对比分析稳态响应振幅的变化而不分析相位差别.以杆件左端p为基准值, 对比分析3种稳态应力幅值|σ_A|_max、|σ_A^*|_max和|σ_N|_max, 定义如下形式的应力无量纲误差：

$ \left. \begin{array}{l} {e_{\rm{A}}} = \frac{{{{\left| {{\sigma _{\rm{A}}}} \right|}_{\max }} - {{\left| {\sigma _{\rm{A}}^ * } \right|}_{\max }}}}{{\bar p}} \times 100\% ,\\ {e_{\rm{N}}} = \frac{{{{\left| {{\sigma _{\rm{A}}}} \right|}_{\max }} - {{\left| {{\sigma _{\rm{N}}}} \right|}_{\max }}}}{{\bar p}} \times 100\% . \end{array} \right\} $

(25)

表 1给出不同ω时杆件中点(x=5m)处的稳态应变幅值.表 2、3分别给出不同ω时杆件x=0 m和x=5 m处的稳态应力幅值和应力的无量纲误差.由表 1~3的计算结果可以看出：1) 当ω取不同值时, 随着α或β的变化, 杆件稳态响应的应变和应力幅值均表现出如上文所述的规律；2) 稳态应力幅值的无量纲误差与ω有关.

4 应力计算差异影响因素分析

从上文分析可以看出σ_N和σ_A^*一致, 导致应力计算误差的原因在于有限元分析软件在计算应力时没有考虑阻尼刚度部分对应力的贡献.本节取λ_α=1, 利用σ_A和σ_A^*对比分析β和ω对应力计算误差的影响.此处仅对比分析|σ_A^*|_max相对|σ_A|_max的应力无量纲误差, 且不分析相位差别.

利用杆件非零基频, 定义无量纲激振频率λ_ω=ω/ω_a.图 9、10分别给出当λ_ω=26和λ_β=1时杆件沿纵向长度x的应力无量纲误差分布曲线.由图 9可见, e在杆件不同位置处起伏分布, 杆件右端因应力为0而误差很小.同时, e随着λ_β的增加而不断加大, 当λ_β=5时, 误差可以超过p的1/5.

图 10表明, 随着ω的增加, e明显增大, 且在杆件内起伏分布更加明显, 这是由应力波在杆件内的波动传播特性决定的.当p不变时, 随着ω增加, 应变率增大, 导致e上升；同时, 随着ω的增加, 应力波波长减小, 导致应变率幅值和应力计算误差幅值在杆件内的起伏分布更加明显.

图 11、12分别给出杆件中点(x=5 m)处e随λ_β和λ_ω的变化曲线.图 11表明, 在不同ω下, e均随λ_β的增加而增大.图 12表明, 杆件中点(x=5 m)处e随着λ_ω的增加起伏变化, 且随着基频ω₁的增加而增大, 原因是E增加引起与应变率相关应力项增大以及杆件内应力波的波动传播.

5 结论

(1) 瑞利阻尼的刚度相关部分反映了材料对动力响应的阻滞作用, 表现为应力不仅与应变有关, 而且与应变率相关；瑞利阻尼的质量相关部分反映了环境对结构体系振动的阻滞作用.瑞利阻尼模型是一种综合反映材料和环境对结构体系动力响应阻滞作用的最简单模型.

(2) 随着α或β的增加, 应变振动幅值减少, 响应滞后, 比较而言, β的影响更明显.在基准瑞利阻尼系数附近, α的变化对动力响应的影响远小于β的变化, 材料阻尼是主要的结构能量耗散源.

(3) 现有的有限元分析考虑了材料阻尼作用和环境阻尼作用对动力平衡的影响, 可以给出可靠的位移、应变及其变化率的解；在应力计算中, 没有考虑材料阻尼的作用, 导致计算应力和实际应力存在差异, 并随着β和ω的增加而明显增大.计算应力不能满足力的边界条件, 因此从理论角度上是不正确的.