土层地震反应显式计算中阻尼矩阵系数的选取

楼梦麟，邵新刚

（同济大学土木工程防灾国家重点实验室，上海200092）

在进行土层地震反应的确定性分析时，通常有两种分析方法，即时域分析和频域分析.在进行时域分析时，一般采用粘滞阻尼假定，即阻尼力与质点的运动速度成正比，通过求解实系数微分方程组，可求得各自由度反应的时程.而在频域分析时，通常采用Fourier变换，通过求解复代数方程组，可求得土层地震反应的稳态解［1］，此时一般采用滞后阻尼假定.这两种分析方法各有不同的适用范围，采用频域法时，可以方便地对线性体系的地震反应进行计算，但对于非线性地震反应的分析只能采用等效线性化方法［2］，无法进行真非线性分析.而时域分析方法能够适应非线性体系地震反应分析的不同要求，既可采用等效线性化的分析方法，也可采用真非线性分析方法［3］.除了简单结构如梁、柱等外，一般工程动力体系地震反应的时域计算多采用直接积分法［4］，直接积分法一般又分为隐式算法和显式算法［5－6］.例如：对于无阻尼体系或阻尼矩阵为一对角矩阵的情况，中心差分方法所给出的格式为显式.但是对于一般非对角线的阻尼矩阵，中心差分方法所给出的格式为隐式.在LS－DYNA［7］程序中对土层非线性动力分析时具有良好的非线性本构模型，且采用时域显式算法，具有较高计算效率，得到较为广泛的应用.然而，程序中土层计算体系的阻尼矩阵采用了仅与质量矩阵成线性比例的假定.由于集中质量矩阵为对角阵，阻尼矩阵也为对角阵，满足了显式算法的要求.

已有研究工作［8］表明，对于深覆盖土层地震反应计算采用土动力学中常用的与刚度矩阵成线性比例的阻尼矩阵将会产生低估土层地震反应的问题.本文通过线性条件下的算例分析，讨论LS－DYNA计算程序中基于质量比例阻尼矩阵显式算法所存在的问题，并进一步以频域反应的解为准确解，通过对比不同工况的时域解，讨论确定质量比例阻尼系数的取值问题.

1 质量比例阻尼矩阵存在的问题

对于实际结构，在进行动力反应分析时，采用有限元方法对计算土层体系进行离散化处理后，可以得到如下动力方程式［9］：

式中：M，C 和K 分别为土层的质量矩阵、阻尼矩阵和刚度矩阵分别为对应于土层下卧基岩面的相对位移向量，速度向量和加速度向量；p（t）为作用于土层体系上的地震等效动力荷载向量，其 中：ej为位移变换向量，aj（t）是作用在基岩上的加速度历程.

在有限元分析中，一般取质量矩阵M 为对角矩阵.当C＝a0M 时（a0为比例系数），阻尼矩阵也为对角矩阵，因此，式（1）中可采用显式算法求解.DYNA程序中建议：

式中：ω1为土层的基频；ξ为土层的振型阻尼比；一般假定土层各阶振型阻尼比相同.以ω1来确定比例系数a0符合人们的常规的认识，即在结构体系地震反应中，第1阶振型频率贡献最大.

设：

式中：Φ 为计算土层的振型模态矩阵；q（t）为广义坐标向量.应用式（3）对式（1）进行坐标变换，并利用振型正交性可得：

令

式中：ηix，ηiy，ηiz分别称为i振型在土层x，y和z方向地震反应中的参与系数（或贡献系数），简称i振型j方向地震参与系数，φi为振型向量.一般来说低阶振型的地震参与系数的绝对值大于高阶振型的地震参与系数.根据式（2），可得第i阶振型的计算阻尼比为

很显然，当假定体系各阶振型阻尼比都为ξ时，按质量比例阻尼矩阵的假定，则高阶振型阻尼比的计算值小于假定值ξ，折减比例系数为ω1／ωi，振型阶序越高，振型阻尼比计算值越小，使得该振型地震反应的计算值就越大.

当土层基频接近于或高于输入地震的主要频率分量时，显然由于土层第1阶振型的振型参与系数大，其对体系地震反应的贡献远大于其他高阶振型，因此当第1阶振型阻尼比与实际阻尼比相等时，在第1阶振型地震反应计算正确得到保证的前提下，其他高阶振型阻尼比偏小时使得土层地震反应的计算误差非常有限.而当计算体系的基频远低于输入地震波的主要频率分量时，很显然与地震波中主要频率分量相近的土层高阶振型将存在共振效应.众所周知，共振放大因子为Di＝1／（2ξi），低估振型阻尼比，将高估共振放大因子，从而使得该振型的地震反应计算值增大.如当ω1／ωi＝1／2时，则Di增大1倍，依此类推，此时高阶振型的共振效应不可忽略.下面通过具体算例，讨论质量比例阻尼矩阵给土层地震反应计算带来的影响.由于适应LS－DYNA 的程序时域显式算法中C＝a0M的规定，本文建议质量比例阻尼系数a0表示如下：

式中：ωs为控制圆频率；fs为控制频率.一般来说，当土层基频接近于或高于输入地震波的主要分量的频率时，按常规取ωs＝ω1，而本文的重点讨论当土层基频远低于输入地震波主要分量的迫振频率时ωs（或fs）的合理取值.

2 数值分析

2.1 土层计算模型

文中土层模型采用港珠澳隧道场地模型的纵剖面为分析对象，水平范围为1 016m，底部为基岩面，左端厚度为85m，右端厚度为112m，假定为均匀土层，土介质的计算参数选取为：弹性模量391 MPa、密度2 000kg·m－3、泊松比0.35，这一土层模型称为土层模型1；为了考察土层厚度对ωs值取值的影响，在这一土层模型的基础上，采用比例按同一比例减小土层左右两端厚度另外形成两个土层模型，分别称为土层模型2和土层模型3.此外，另建立一个4层土介质组成的水平均匀分层土层模型（简称土层模型4）与均质斜坡土层相对应，自表层向下各层土介质的计算参数为：弹性模量185.1，391，561.3和2 492.4 MPa，质量密度1 870，2 000，2 220和2 100 kg·m－3，泊松比0.4，0.35，0.35和0.35.土层模型的有限元网格如图1所示.

图1 土层有限元网格图Fig.1 Finite element mesh of soil layer

4个土层计算模型的几何尺寸和土层前5阶自振频率见表1，从表中可以看出，土层模型1和土层模型4的自振频率较为接近.

表1 土层计算模型特征长度与自振频率Tab.1 Characteristic length and natural frequency of the soil layer calculation model

2.2 基岩面输入地震波

为了比较全面地考察基岩面输入地震波的影响，在土层地震反应计算中共选用了5条地震波，包括该工程场地地震安全性评估报告给出的基岩人工波和El Centro、Kobe、Taft、汶川基岩波等4条地震实测波，它们的时程、Fourier幅值谱和反应谱分别如图2—6所示.

图2 基岩人工波Fig.2 Bedrock artificial wave

为了考察地震波频谱分布对质量比例阻尼系数取值的影响，选取了一下几个频率参数作为地震波的特征频率参数，分别为：①反应谱峰值频率fR；②反应谱谱曲线重心频率fRg；③Fourier幅值谱峰值频率fF；此外还引入考虑结构基频和地震波特征频率综合影响的频率参数fa和fb；④fa＝（f1＋fR）／2，f1＝w1／2π；⑤fb＝（f1＋fF）／2；⑥文献［5］中滞后阻尼转换频率fc.计算中输入地震波的加速度峰值都调整为0.1g（g为重力加速度），土层各阶振型频率阻尼比取为0.08，对应于5条输入地震波的不同特征频率值见表2—4.

分析表中5条地震波的特征频率参数fR，fF和fRg可以看出所选的汶川基岩地震波由于频谱成分分布很广，所以fR、fF值很高，是比较特殊的一条地震波，表明高频成分很丰富，同样基岩人工波中高频成分也较为丰富.

表2 对应于土层模型的特征频率值Tab.2 Characteristic frequency value corresponding to the model of the soil layer Hz

表3 土层模型1中地表中点加速度与位移反应幅值Tab.3 Acceleration and displacement response amplitude of the surface midpoint on Soil Model 1

表4 土层模型1中地表加速度与位移计算误差Tab.4 Acceleration and displacement errors of calculation of the surface midpoint on Soil Model 1 %

2.3 土层地震反应分析

对于同一土层有限元模型，分别采用时域和频域分析方法，频域计算时采用了滞后阻尼假定，因此以频域计算所得土层反应为准确值，作为不同时域分析结果的比较基准.在时域分析中，不同工况计算所得的反应结果与频域的计算结果相比较，时域近似分析结果的相对误差可表示为

式中：a和a＊为土层同一节点处地震位移或加速度反应峰值，其中a为时域分析结果，a＊为频域分析结果.

表3和表4中列出了土层模型1在不同地震波激励下的土层地表中点加速度和位移反应幅值及相应的时域计算结果的误差.

从表4中可以看出：由于土层模型1 的基频远低于表2中5条地震波的特征频率fR、fF以及fRg，所以当取ωs＝ω1形成土层地震反应计算中的阻尼矩阵时，显然高估了土层地表的加速度峰值以及位移峰值.相比之下取其他特征频率时，土层地震反应峰值计算结果要好于ωs＝ω1时的结果.进一步对比表2和表4中的数据可以看出，总体上ωs＝ωa＝2πfa时，计算精度要好一些，且基本能保证地震加速度峰值为计算误差为正，即计算值略大于实际值.这里主要考虑到：一般来说土层地震反应计算的目的是为地面结构地震反应分析提供结构基底地震动加速度时程.

图7和图8中比较了在El Centro波和汶川波激励下的土层模型1 地表中点加速度反应时程曲线，从图中可以看出，采用ωs＝ω1形成阻尼矩阵后，计算所得的土层地表中点加速度反应时程与频域计算所得的时程相差较大，而采用ωs＝ωa形成阻尼矩阵后（ωa＝2πfa），计算所得的土层地表中点加速度反应时程得到了较好改善.其他地震波激励下的土层地表中点加速度反应时程曲线与El Centro波和汶川波激励下的情况类似，限于篇幅，不一一列出.

表5中列出了土层模型2在不同地震波激励下土层地表中点加速度和位移反应的时域计算结果的相对误差.

表5 土层模型2中地表加速度与位移计算误差Tab.5 Acceleration and displacement errors of calculation of the surface midpoint on Soil Model 2 %

从表5可以看出：由于土层模型2的基频总体 上略低于表2中5条地震波的特征频率，因此误差变化不及表4 中显著，但取ωs＝ω1时，依然高估了土层的地震反应值，从加速度峰值误差分布来看，取ωs＝ωa时误差水平最低.采用ωs＝ω1及ωs＝ωa形成阻尼矩阵后，El Centro波和汶川波激励下的土层地表中点加速度反应时程在相位与波形上非常接近，区别体现在峰值上，且在土层模型2中的误差不如模型1显著了.限于篇幅时程图不再一一列出.

表6中列出了土层模型3在不同地震波激励下的土层地表中点加速度和位移反应的时域计算结果的误差.

表6 土层模型3中地表加速度与位移计算误差Tab.6 Acceleration and displacement errors of calculation of the surface midpoint on Soil Model 3 %

从表6可以看出，由于土层基频总体上高于或远高于（汶川地震波除外）表2中5条地震波的特征频率，因此当取ωs＝ω1时总体上看计算误差最小.这与目前LS－DYNA 程序中约定的方法相一致.采用ωs＝ω1及ωs＝ωa形成阻尼矩阵后，El Centro波和汶川波激励下的土层地表中点加速度时程在波形上同样接近，区别依然体现在峰值上.且在土层模型3中取基频的计算结果的误差水平已经和ωs＝ωa时相当.同样，限于篇幅时程图不再一一列出.

表7中列出了土层模型4在不同地震波激励下的土层地表中点加速度和位移反应的时域计算结果的误差.

表7 土层模型4中地表加速度与位移计算误差Tab.7 Acceleration and displacement errors of calculation of the surface midpoint on Soil Model 4 %

从表7可以看出：均匀分层土层模型4所体现出来的规律与模型1大体相同，并没有因为土层的分层及坡度的改变发生太大的变化，总的来说加速度反应误差在10%以内略好于均匀土层1.体现了该方法的实际应用价值.

采用ωs＝ω1及ωs＝ωa形成阻尼矩阵后，El Centro波和汶川波激励下的土层地表中点加速度时程反应的规律与模型1相类似.

此外，从表4—7可以看出，式（8）中ωs取值对加速度反应的影响远大于位移反应的影响，对于3种不同类型的土层，ωs取为土层基频时，土层地表位移峰值的计算误差最小，且误差都为正.其原因是：对应于加速度简谐波a（t）＝A0sinωt（A0为振幅）的位移波为u（t）＝－（A0／ω2）sinωt＝v0sinωt.很显然，随着激励频率ω的升高，位移简谐波幅值按激励频率平方减小.因此，土层地震位移反应更敏感于低频激励，因此准确计算土层第1阶振型位移反应是最为关键，当ωs＝ω1时可保土层位移反应计算的合理性.

3 结论

从本文所选三种不同土层计算模型在不同输入地震波的地震反应数值结果比较可看出：

（1）在土层地震反应时域显示算法中，所形成质量比例阻尼矩阵对土层地震反应有很大影响，当土层基频等于或高于地震波的主要特征频率时，常规做法是取土层基频来形成阻尼矩阵是可行的.

（2）当土层基频略低于地震波的主要特征频率值时，用土层基频来形成比例质量矩阵时，会高估土层的地震加速度反应.本文算例中，地表加速度峰值有时高估达30%和80%不等.

（3）当土层基频远低于地震波的主要特征频率值时，采用常规方法形成质量比例阻尼矩阵式，将更加高估土层地震加速度反应值，本文算例中地震加速度峰值误差有时高达60%和90%不等.

（4）由于地震波为不规则波，频率分量丰富，很难用一个频率来形成质量比例阻尼矩阵以合理反映各阶振型阻尼比的数值，因此应重视在形成质量比例阻尼矩阵中频率的选择.本文建议按式（8）来形成质量比例阻尼系数a0，其中频率ωs将针对土层和输入地震波的特征频率来选择.一般来说，①当土层基频高于输入地震波的特征频率时，可取ωs为土层基频，与目前常规认识一致；②当土层基频接近于特别是远低于地震波的特征频率时，建议ωs取为土层基频和输入地震波反应谱峰值频率的算术平均值ωa，基本可保证计算所得土层加速度峰值略高于真实解.

（5）由于土层地震反应计算的主要目的之一是为地面结构提供结构基底输入地震加速度，为此本文上述结论主要是针对土层地震加速度反应计算，而对于地震位移反应计算时，ωs仍取为土层基频为好，因为土层第1阶振型的位移反应在土层地震位移反应中起控制作用.

［1］ 廖振鹏.近场波动问题的有限元解法［J］.地震工程与工程振动，1984，4（2）：1.LIAO Zhenpeng.A finite element method for near－field wave motion in heterogeneous materials［J］.Earthquake Engineering and Engineering Vibration，1984，4（2）：1.

［2］ Martin P P，Seed H B.One－dimensional dynamic ground response analyses ［J］.Journal of Geotechnical Engineering Division，ASCE，1982，108（7）：935.

［3］ 李小军，廖振鹏，关慧敏.粘弹性场地地形对地震动影响分析的显式有限元－有限差分方法［J］.地震学报，1995，17（3）：362.LI Xiaojun，LIAO Zhenpeng，GUAN Huimin.An explicit finite element－finite difference method for analyzing the effect of visco－elastic Local topograph on the earthquake motion［J］.Acta Seismologica Sinica，1995，17（3）：362.

［4］ 张汝清，殷学纲，董明.计算结构动力学［M］.重庆：重庆大学出版社，1987.ZHANG Ruqing，YIN Xuegang，DONG Ming.Computational structural dynamics［M］.Chongqing：Chongqing University Press，1987.

［5］ 李小军，廖振鹏，杜修力.有阻尼体系动力问题的一种显式差分解法［J］.地震工程与工程振动，1992，12（4）：75.LI Xiaojun，LIAO Zhenpeng，DU Xiuli.An explicit finite difference method for viscoelastic dynamic problem ［J］.Earthquake Engineering and Engineering Vibration，1992，12（4）：75.

［6］ 杨柏坡，陈庆彬，王琦萍.地震波动研究中的显式有限元法［C］／／第三届全国地震工程会议论文集.大连：大连理工大学出版社，1990：211－216.YANG Baipo，CHEN Qingbin，WANG Qiping.An explicit finite element method of the seismic wave study［C］／／Proceedings of the Third National Conference on Earthquake Engineering.Dalian：Dalian University of Technology Press，1990：211－216.

［7］ 尚晓江，苏建宇，王化峰.ANSYS／LS－DYNA 动力分析方法与工程实例［M］.北京：中国水利水电出版社，2008.SHANG Xiaojiang，SU Jianyu，WANG Huafeng.ANSYS／LSDYNA dynamic analysis method and projects［M］.Beijing：China Water ＆Power Press，2008.

［8］ 楼梦麟，潘旦光.滞后阻尼在土层时域分析中的应用［J］.同济大学学报：自然科学版，2004，32（3）：281.LOU Menglin，PAN Danguang.Hysteretic damping application in time domain analysis of soil layer［J］.Journal of Tongji University：Natural Science，2004，32（3）：281.

［9］ Clough R W，Penzien J.Dynamics of structures［M］.New York：Mc－Graw Hill Inc，1993.