在地震反应过程中，阻尼是影响反应结果的重要因素^[1-2].为进行强震下结构的弹塑性时程反应分析，需要构建相应的阻尼矩阵.在各种阻尼模型中，Rayleigh阻尼由于计算简便而得到广泛应用.

采用矩阵表达，Rayleigh阻尼矩阵C可表达为质量矩阵M和刚度矩阵K的线性组合，即

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

(1)

式中α和β分别为质量和刚度比例阻尼系数.α和β这两个系数可以通过指定两个参考频率(ω_i和ω_j)及其阻尼比(ζ_i^*和ζ_j^*)进行计算.显然，Rayleigh阻尼除两阶参考频率的阻尼比外，其他各阶模态的阻尼比有一定的误差存在.因此，所构造的Rayleigh阻尼是否合理依赖于选取的两个参考频率是否合理.在结构地震反应分析时，Rayleigh阻尼参考频率的选取有三类常用方法.

第一类方法是直接从结构的动力特性角度选择参考频率.Idriss等^[3]直接利用体系基频ω₁建立阻尼矩阵，此时所构造的Rayleigh阻尼高估了除基频以外的所有阻尼，从而导致结构动力反应偏小^[4].对于简单的建筑结构，常用两个低阶自振频率作为参考频率^[5].但对复杂结构而言，对结构有显著贡献的模态数目多达几百阶^[6].此时，如果用两个低阶模态建立Rayleigh阻尼矩阵将使结构高阶模态的阻尼比偏大，从而导致低估结构的地震反应.Chopra^[7]定性的指出，选择的两阶参考频率应使对结构反应有显著贡献模态的阻尼比取值合理.Clough等^[8]建议ω_i=ω₁，ω_j从对结构动力反应有显著贡献的高阶振型中选取.Youssef等^[9]的计算结果表明，最优参考频率的阶数随土层深度变化而变化.

第二类方法是令ω_i=ω₁，第二个参考频率ω_j由地震波的频谱特性进行确定.Hudson等^[10]令ω_j为地震波卓越频率ω_e和结构基频的奇数倍.即

$ {{\omega }_{\mathit{j}}}\rm{=}\mathit{n}{{\omega }_{\rm{1}}}, $

(2)

式中n为大于ω_e/ω₁的奇数.楼梦麟等^[11]建议以输入地震波加速度反应谱的峰值频率作为ω_j进行深覆盖土层的地震反应计算.

第三类方法采用优化算法直接得到Rayleigh阻尼系数.杨大彬等^[12]以振型的峰值应变能系数作为权重函数，建立了加权最小二乘法.潘旦光等^[13-14]分别以结构峰值位移误差和基底剪力误差最小为目标函数，提出了一种求解Rayleigh阻尼系数的优化求解方法.

为比较上述三类Rayleigh阻尼系数计算方法的特点及对结构动力反应的影响.本文将以大跨网壳结构水平方向和竖向地震反应为例，讨论不同方法所得阻尼矩阵对网壳顶点位移及基底剪力的影响.除此以外，本文还从Rayleigh阻尼构建角度讨论了地震波主要频率成分及其选取范围的问题.

1 大跨网壳有限元模型及输入地震波 1.1 有限元模型

结构计算模型如图 1所示，上部为85.2 m长，61.8 m宽单层网壳结构，下部由8组支架支撑整个结构.结构平面关于x轴对称，关于y轴不对称，z为竖向坐标.结构主要构件如下：网壳内部杆件和环梁分别为300 mm×50 mm×8 mm×12 mm和1 000 mm×600 mm×30 mm×30 mm的空心方管，下部斜撑和底柱分别为300 mm×20 mm和800 mm×35 mm的空心圆管.所有杆件材料为Q235，采用梁单元进行建模.模型的模态特征如表 1所示.表中r_y和r_z分别表示y和z方向的振型参与质量比，s_y和s_z分别表示y和z方向的累积振型参与质量比.部分典型模态图如图 2所示.并设各阶模态的阻尼比为0.02.

1.2 输入地震波

为比较不同类型地震波对Rayleigh阻尼系数计算的影响，选用表 2中的3条不同场地类型地震波分别作为柱根部的水平和竖向地震输入.输入地震波的加速度时程如图 3所示，加速度时程的幅值统一调整为0.35 m/s².

2 Rayleigh阻尼系数计算方法

由模态分析可知，单层网壳结构在y方向(水平)地震输入时，前45阶模态的振型参与质量即超过90%，第一阶模态的反应贡献具有绝对统治地位.对于z方向(竖向)地震输入时，前205阶模态的振型参与质量超过90%，第三阶模态的反应贡献最大，但不具有统治地位.对于水平地震反应，结构的显著贡献模态少；而竖向地震反应，结构的显著贡献模态多，且没有具有绝对统治地位的模态.同时，在水平地震反应时，结构的第一个显著贡献模态的频率小于大部分地震波的卓越频率，而竖向地震反应时，结构的第一个显著贡献模态的频率高于部分地震波的卓越频率，因此，讨论y方向和z方向的地震反应，相当于讨论了两种结构类型下Rayleigh阻尼系数的计算问题.不同方法计算阻尼系数的差别在于两个参考频率选取的不同，为此，下面讨论单层网壳结构y方向和z方向地震反应时三类方法所得的参考频率及对结构地震反应的影响.

2.1 第一类方法

对于y方向地震输入，一种是根据经验选择ω_i=ω₁和ω_j=ω₁₉(记为i=1 & j=19)，另一种是直接令ω_i=ω₁和ω_j=ω₂(i=1 & j=2).对于z方向地震输入，分别考虑两种组合ω_i=ω₃和ω_j=ω₄₅(i=3 & j=45)及ω_i=ω₁和ω_j=ω₂(i=1 & j=2).并在后面的讨论中将这种方法称为传统方法.

2.2 第二类方法

在第二类方法中，y方向和z方向地震输入的第一个参考频率分别选为ω₁和ω₃.第二个参考频率可基于反应谱和Fourier谱进行选择.图 4为3条地震波的位移、拟速度和拟加速度反应谱.

在利用反应谱选择Rayleigh阻尼第二个参考频率时，常直接选取反应谱峰值频率^[15].但这种方法存在的一个问题是：当结构的基频高于反应谱峰值频率时，构建的Rayleigh阻尼将使所有的高阶模态阻尼比大于真实的阻尼比，从而使计算结果偏小，这一点在位移反应谱中尤其明显，因为位移反应谱的峰值常出现在长周期.譬如，自振周期在0~6 s范围内，El Centro、Northridge和Tianjin波的位移反应谱峰值频率分别为0.352、0.410和0.568 Hz，这些峰值频率都小于基频，此时，都无法构造出合理的Rayleigh阻尼.事实上，反应谱反映的是在一个确定地震波作用下不同自振周期单自由度体系的最大反应，然后，以一定组合规则得到结构的反应.结构高阶模态的周期T_i都是小于基频周期T₁的，即T_i＜T₁ (i=1, 2, …).因此，[0, T₁]区间的反应谱才参与结构反应的计算.从这个角度看，以[0, T₁]区间内反应谱的峰值频率作为第二个参考频率更合理.表 4和表 5中列出的拟加速度反应谱峰值频率f_Ra、拟速度反应谱峰值频率f_Rv和位移反应谱峰值频率f_Rd都是指[0, T₁]区间内的峰值频率.

考虑到多自由度体系是一系列单自由度体系反应的叠加，因此，部分学者将地震波的主要频率定义为反应谱形心的频率^[15].反应谱的形心频率f_Rg-i定义如下

$ \frac{1}{{{f}_{{\rm{Rg}}-i}}}=\frac{\int_{0}^{{{T}_{\rm{R}}}}{T{{S}_{i}}\left( \zeta, T \right)\rm{d}\mathit{T}}}{\int_{0}^{{{T}_{\rm{R}}}}{{{S}_{i}}\left( \zeta, T \right)\rm{d}\mathit{T}}}, $

(3)

式中:S_i(ζ, T)为反应谱，i分别取位移、拟速度和拟加速度.显然f_Rg-i是积分范围T_R的函数.图 5为不同地震波反应谱f_Rg-i随积分区间变化的曲线.由图可知，随着积分区间的增加，形心点的频率逐步降低.这是长周期部分反应谱影响的必然结果.由于只有[0, T₁]区间的反应谱参与结构反应的计算，因此，计算形心频率的反应谱范围定义为[0, T₁]更合理.表 4和表 5中列出拟加速度反应谱形心频率f_Rg-a、拟速度反应谱形心频率f_Rg-v和位移反应谱形心频率f_Rg-d都是指[0, T₁]基本区间不同地震波反应谱的形心频率.

图 6为不同地震波Fourier谱.Fourier谱表明地震波中含有哪些频率分量，及哪些频率分量振幅大.因此，将分量振幅最大的频率称为卓越频率，并记为Fourier谱峰值频率f_F.实际地震波的Fourier谱呈剧烈起伏的锯齿状，为避免Fourier谱中个别尖刺的影响，常将平滑化方法所得的峰值频率作为地震波卓越频率.本文将平滑化后Fourier谱的峰值频率记为f_P.图 6中的平滑化曲线是采用矩形脉冲窗的结果，窗的带宽取为1.2 Hz.

在得到地震波的卓越频率后，可采用式(2)计算Rayleigh阻尼系数.当ω_e取为Fourier谱的峰值频率f_F时，称为Ⅰ-1方法；当ω_e取为平滑化Fourier谱的峰值频率f_P时，称为Ⅰ-2方法.且当ω_e＜ω₁时，令n=1.

2.3 第三类方法

对于第三类方法，以文献[14]提出的约束优化解法作为对比.该方法首先建立代数方程组：

$ \left[{\begin{array}{*{20}{c}} \mathit{\boldsymbol{G}}&{{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_s}}\\ {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_s^{\rm{T}}}&0 \end{array}} \right]\left\{ \begin{array}{l} \mathit{\boldsymbol{a}}\\ \lambda \end{array} \right\} = \left\{ \begin{array}{l} \mathit{\boldsymbol{H}}\\ \zeta _s^ * \end{array} \right\}, $

(4)

式中：G= ΩwΩ^T，H= Ωwy，w=[w_ij]，w_ij=ρ_ijM_iM_jS′_a×(ζ_i, ω_i)S′_a(ζ_j, ω_j)，y={ζ₁^*, ζ₂^*, …, ζ_N^*}^T，ζ_i^*(i=1, 2, …, N)是指各阶模态的精确阻尼比.

Ω=[Ω₁, Ω₂, …, Ω_N]，a={α β}^T为Rayleigh阻尼系数组成的向量，Ω_i=0.5{ω_i^-1 ω_i}^T(i=1, 2, …, N)中的ω_i为结构的第i阶固有频率，y方向和z方向激振下N分别取为45和205.S′_a(ζ_n^*, ω_n)为拟加速度反应谱关于阻尼的导数，M_i为第i阶模态的振型参与质量百分比，Ω_s=0.5{ω_s^-1 ω_s}^T是指约束第s阶的模态阻尼比等于ζ_s^*.对于y方向和z方向地震输入，约束模态分别取为ω₁和ω₃.

根据以上讨论，共形成12种Rayleigh阻尼系数的计算方法，其中传统方法在算例比较时，采用两种频率组合.表 3列出了各种计算方法的简化名称和所对应的类型.下面讨论各种方法所得Rayleigh阻尼对单层网壳结构地震反应的影响.

3 数值计算结果 3.1 计算误差公式

为比较不同方法的计算误差，进行结构的线弹性地震反应分析.同时，以前300阶模态振型分解时程分析方法的计算结果为精确解，记为r^*，采用Rayleigh阻尼模型所得的近似解记为r，则Rayleigh阻尼模型计算结果的相对误差为

$ \mathit{e} = \frac{{\left| r \right| - \left| {{r^ * }} \right|}}{{\left| {{r^ * }} \right|}} \times 100\% . $

(5)

在m条地震波作用下，各反应量的平均误差可采用以下两式进行统计：

$ {{\mathit{\bar e}}_{\rm{r}}} = \frac{1}{m}\sum\limits_{i = 1}^m {{e_i}, \;} {{\mathit{\bar e}}_{\rm{a}}} = \frac{1}{m}\sum\limits_{i = 1}^m {\left| {{e_i}} \right|.} $

(6)

式中：e_r用于度量不同方法所得计算结果是偏大还是偏小，而e_a用来表示这个方法的计算精度.

3.2 不同Rayleigh阻尼模型计算误差

El Centro、Northridge和Tianjin三条地震波分别沿y方向(水平)和z方向(竖向)输入下，前述12种方法的参考频率如表 4和表 5所示.在y方向地震输入时，仅统计y方向的地震反应，z方向地震输入时，仅统计z方向的地震反应.图 1中A点的位移计算误差和基底剪力误差如表 6~9所示.El Centro波作用下A点位移反应的精确解及其Fourier谱如图 7所示，基底剪力的精确解及其Fourier谱如图 8所示.表和图中的u_Ay和u_Az分别表示A点的y方向和z方向的位移.F_y和F_z分别表示y方向和z方向的基底剪力.

由计算结果可看出：

1) 除i=1 & j=2组合的传统方法进行竖向地震反应计算外，其余方法所得A点位移u_Ay和u_Az的平均误差基本相同.对于y方向地震输入，观察表 1和图 7的结果可以发现，第一模态的位移反应远大于其他模态的反应.对于z方向地震输入，虽然第三模态的振型参与质量的优势不是特别明显，但是，u_Az的Fourier谱表明，第三模态的位移反应具有绝对主导地位.此时，对于Rayleigh阻尼矩阵而言，当y方向地震输入时，只要ζ₁=ζ₁^*；当z方向地震输入时，只要ζ₃=ζ₃^*；u_Ay和u_Az的计算精度即可满足精度要求.而任意选择前两阶模态组合计算Rayleigh阻尼系数，进行z方向地震反应计算，计算误差大且偏小.

2) 对于基底剪力F_y和F_z，以f_Rv、f_Rd、f_Rg-v、f_Rg-d为第二参考频率的计算结果都偏小.对比反应谱的峰值频率和形心频率可以发现f_Ra＞f_Rv＞f_Rd，以及f_Rg-a＞f_Rg-v＞f_Rg-d.这是由于反应谱在大于2 Hz区域为加速度敏感区，小于0.3 Hz为位移敏感区，中间区段为速度敏感区^[7].而由图 8基底剪力的Fourier谱可知，高阶模态对基底剪力的反应也有显著贡献.譬如f₁₉=2.905 Hz和f₄₅=4.443 Hz模态的反应分别对F_y和F_z有显著影响，而基于速度反应谱的f_Rv、f_Rg-v和位移反应谱的f_Rd、f_Rg-d的第二个参考频都远小于f₁₉和f₄₅，导致这两阶模态的阻尼比偏大而低估了高阶模态的反应.基于加速度反应谱峰值频率和形心频率的计算误差相对较小.

3) f_F作为第二个参考频的计算结果偏小.这是由于Fourier谱和拟速度反应谱类似，因此，f_F和f_Rv的计算结果也类似.以平滑化的峰值频率f_P作为第二个参考频的计算结果要略优于f_F的计算结果，但无法改变计算结果偏小的特点.同时，当f_F或f_P小于结构的基频时，直接用卓越频率作为第二个参考频率的计算误差更大.这是由于这种方法除基频的阻尼比外，其余模态的阻尼比都大于精确解而导致计算结果偏小.因此，对于地震波卓越频率小于结构基频的地震波，采用地震波卓越频率作为参考频率不合理.

4)Ⅰ-2方法的平均误差小于Ⅰ-1方法，这表明采用平滑化Fourier谱所得卓越频率进行Idriss方法^[10]第二个频率的计算更合理.而且，由于f_F易受Fourier谱中个别低频尖刺的影响，而使计算结果离散性大.同时，当地震波的卓越频率小于结构的基频时，此时，Idriss方法实际上就是基于基频建立阻尼矩阵，Rayleigh阻尼高估了高阶模态的阻尼比，由此使结构的地震反应偏小，因此，当地震波的卓越频率小于结构的基频时，Ⅰ-1和Ⅰ-2方法都是不合理的.当地震波的卓越频率大于结构的基频时，可采用Ⅰ-2方法建立Rayleigh阻尼.

5) 对于优化方法，无论是y方向地震输入还是z方向地震输入，由于优化算法中考虑了结构动力特性、地震波频率特性的影响，因此，对于两个方向的地震输入和所有地震波的计算结果误差都小且稳定.而且，优化方法是直接得到Rayleigh阻尼系数，避免了人为选择两阶参考频率的任意性，适用于不同工程结构的Rayleigh阻尼构建.

6) 对于传统方法，本文所选的y方向和z方向的频率组合i=1 & j=19和i=3 & j=45实际上和优化分析方法的频率基本相同，因此，对于富有经验的计算人员，直接选取两阶合理的参考频率用于Rayleigh阻尼计算也是可行的.但是，任意选择前两阶模态进行Rayleigh阻尼计算所得计算结果偏小，且绝对误差最大，是不合理的计算方法.

4 结论

在大型复杂结构的非线性地震反应分析时，常需建立Rayleigh阻尼矩阵进行直接积分法计算.由于参与结构振动的模态多且复杂，因此，如何选取合理的参考频率是一个需要仔细斟酌的事情.本文以一个长85.2 m，宽61.8 m的大跨屋盖为例，对比分析了12种国内外Rayleigh阻尼系数计算方法对计算结果的影响，由理论分析和数值计算可得出以下结论：

1) 对于第一类方法，工程技术人员如果对结构的动力反应有充分认识，可直接指定两阶参考频率.但是，任意选择前两阶模态进行Rayleigh阻尼计算，易造成计算结果偏小，是不合理的计算方法.

2) 对于第二类方法，建议采用[0, T₁]区间的拟加速度反应谱峰值频率或形心频率、Ⅰ-2方法作为Rayleigh阻尼系数计算的第二个参考频率.

3) 优化方法所得Rayleigh阻尼系数是综合考虑结构的动力特性、输入地震波频谱特性的综合结果，直接得到Rayleigh阻尼系数，避免了基于经验指定两阶频率的任意性，且计算误差较小，适合于工程结构的计算与分析.