绝对位移直接求解虚拟激励法采用附加振型求解的理论研究

第一作者李永华男，博士，副教授，1972年生

通信作者桂国庆男，博士，教授，1966年生

绝对位移直接求解虚拟激励法采用附加振型求解的理论研究

李永华1，桂国庆2,1，廖宇1

(1.南昌大学建筑工程学院，南昌330031； 2.井冈山大学建筑工程学院，江西吉安343009)

摘要：为提高已有绝对位移直接求解虚拟激励法的求解效率及精度，对结构支座节点附加大质量块并释放支座约束，采用基于附加振型的振型分解法进行求解，附加振型来源于附加大质量块。对附加振型、常规振型的振型特性和振型贡献及计算效率进行了理论分析，分析表明：通过构建少量附加振型就可精确捕获结构拟静位移，且只需采用与相对运动法同等数量的常规振型就可精确捕获结构动态相对位移，与Wilson位移输入模式的绝对位移振型分解法相比，在保证计算精度的条件下可大大减少绝对位移求解所需的振型数。此外，对附加振型绝对位移求解法阻尼误差进行了理论分析，指出了附加振型法不存在由于阻尼假定不同而引起的误差，计算精度及效率均优于基于完全法的绝对位移直接求解法。

关键词：附加振型；虚拟激励法；绝对位移直接求解；振型分解；误差分析

基金项目：国家自然科学

收稿日期：2013-10-09修改稿收到日期：2014-04-30

中图分类号:TU311.3文献标志码：A

Pseudo excitation method based on directly solving absolute diaplacements with addtional modes

LIYong-hua1,GUIGuo-qing1,2,LIAOYu1(1.School of Architecture and Civil Engineering, Nanchang University, Nanchang 330031, China; 2. School of Architecture and Civil Engineering, Jinggangshan University, Ji’an 343009, China)

Abstract:Structures were solved using the pseudo excitation method with additional modes through adding big masses on support nodes in order to improve the calculation accuracy and efficiency of the old pseudo excitation method based on directly solving absolute displacements. The additional modes were caused by the added big masses. The contribution and characterisitics of additional modes and the conventional ones as well as the calculation efficiency of the new algorithm were analyzed. It was shown that constructing only few additional modes can capture the pseudo-static displacement responses of structures accurately, and adopting the same number of conventional modes obtained with the relative motion method can capture the same accurate dynamic relative displacement responses of structures; compared with the displacement input model presented by Wilson, the proposed method can lead to a significant reduction of the required modes for solving absolute displacements with the same accuracy; the calculation errors caused by different damping assumings don’t exist in the new algorithm, so it is superior to the old pseudo excitation method based on directly solving absolute displacements in accuracy and efficiency.

Key words:additional modes; pseudo excitation method; directly solving absolute displacement; mode decomposition; error analysis

李永华[1]提出了一种随机振动绝对位移直接求解的虚拟激励法，该方法通过在基底附加大质量块直接对绝对位移进行一次性求解，避免了在构建虚拟激励的过程中需提取整体结构质量和刚度矩阵的繁琐，且可在通用有限元软件中直接实现,大大方便了多点激励随机振动工程应用[2-4]，该算法已得到较好的应用和发展。Wilson[5]采用位移输入模式进行绝对位移求解时，发现振型分解法需要大量甚至全部高阶振型组合才能达到所需的计算精度，受此影响李永华等在研究过程中，没有采用振型分解法进行求解，而是采用完整矩阵法(以下简称完全法)进行谐响应求解。

绝对位移直接求解的虚拟激励法在采用完全法进行求解时面临两个问题：①计算时间：随机振动虚拟激励法需对大量离散频率点(上百或上千个)构建虚拟谐响应求解，当结构自由度较大时，采用完全法谐响应求解时间较长；②误差问题：误差问题分为两方面，一方面是由于瑞利阻尼与振型阻尼之间的转换所引起的误差:采用完整矩阵法求解时，需形成整体结构阻尼矩阵(通常采用瑞利阻尼模型构建)，而瑞利阻尼系数只能由两阶自振频率和对应的阻尼比来确定，如指定ξi=ξi+1=0.05后，则更高振型(>i+1)或更低振型(

赵银庆[3]在绝对位移直接求解过程中，尝试了采用振型分解法进行求解，发现采用全部振型进行组合时，绝对位移直接求解的虚拟激励算法与传统虚拟激励算法计算结果基本一致。但采用少量振型进行组合时两种算法计算结果是否一致并没有进行比较,也缺少相关的理论分析。如果必须采用大量或全部振型进行组合才能达到所需的计算精度，采用振型分解绝对位移直接求解法就没有任何优势。

张林等[9]采用大质量法，对FR-4印制板进行算例分析时,发现是否选取与质量点的刚体运动对应的振型决定了结果是绝对解,还是对质量点的相对解。

本文在上述文献研究的基础上，拟进一步从理论上分析绝对位移直接求解的虚拟激励法，当采用振型分解进行求解时面临的两个核心问题：①达到工程应用所需的精度，需要的振型数量;②由于振型分解算法通常采用振型阻尼，在振型阻尼模型下，绝对位移直接求解的虚拟激励法与传统的虚拟激励算法由于阻尼假定不一致会引起一定误差。

1基本思路

多点激励地震响应求解过程中，动力方程求解方法可分为相对运动法和绝对位移直接求解法。相对运动法将结构内部节点绝对位移分为拟静位移和动态位移两部分分别进行求解，这里拟静位移为支座移动引起的内部节点位移，动态位移为支座移动导致的内部节点产生的惯性力所引起的内部节点位移[10-11]。

众所周知，地震作用物理上是由于地面的振动，带动上部结构发生振动，无论是一致激励还是多点激励，地面实际上是运动的，图1为上述两种算法之间的物理模型与动力学计算模型示意图。

图1　物理模型与动力学计算模型示意图 Fig.1 Comparison of physical and dynamic model

一致地震动作用下，通常不采用绝对位移直接求解法进行求解，而是采用相对运动法进行求解,如图1(a)所示，即固定支座，将地震动加速度作用施加于结构内部节点上，以内部节点与支座节点之间的相对位移(动态位移)为变量建立动力方程。

相对运动法的最大优点是：由于动力方程中不含拟静位移，只有动态相对位移，而动态相对位移在物理上与振型密切相关，多数建筑结构在地震作用下的动力响应，动态相对位移主要由低阶振型确定，高阶振型对结构响应的贡献不显著。因此，通常只需采用较少数量的低阶振型的贡献就可得到工程所需的计算精度(90%以上的振型参与质量)，大大提高了计算效率[12]。

绝对位移直接求解法有3种加载形式，如图1(b)所示，其中:a模式为常见的Wilson位移输入模式，即把地震作用转化为首层力和首层弯矩进行施加，详细理论见文献[5]第22章；b模式为位移直接输入模式，多数商用有限元软件都不支持位移直接输入模式的振型分解求解法，相关算法理论及误差分析详见文献[8]；c模式为本文重点研究的大质量法输入模式。

Wilson位移输入模式必须采用大量甚至全部高阶振型才能达到绝对位移响应所需的计算精度,下面分别从物理和数学上来进行分析：

(1)由于绝对位移中包含拟静位移，而拟静位移在物理上与振型无关，只与结构刚度矩阵有关，与结构质量矩阵并没有任何关系，因而物理上无法保证采用少数低阶振型的组合就可确保拟静位移计算结果的精度。

(2)从数学上讲，任何n个线性无关的n维向量都可以是向量空间Rn的一个基[13]，也就是说，如果想用常规振型来表示空间结构的拟静位移，空间结构有n个自由度，就必须采用n个振型去组合，数学上才能保证计算精度，这就解释了Wilson教授的算例中，即使是一个很简单的结构模型，采用Wilson位移输入模式时，需要大量甚至全部振型才能达到绝对位移所需的计算精度。

从上面分析可以看出，采用Wilson位移输入模式，想通过少量振型的组合来获取结构的绝对位移是行不通的，必须进行一定的修正[14]。

大质量法输入模式的振型与Wilson位移输入模式的振型是不同的，由于大质量块的引入，相应的支座约束被取消，增加了部分振型(本文称这部分增加的振型为附加振型，支座固定模式的结构振型为常规振型)，故振型参与的理论基础不一样。本文的基本思路为：

(1)能否通过释放结构支座自由度约束并附加大质量单元构建少量附加振型，且这些附加振型就能精确捕获结构的全部拟静位移，其余的常规振型只需贡献结构的动态相对位移。

(2)由于常规振型只需贡献动态位移，这就类似于相对运动法中的动态位移计算，能否采用与相对运动法相同数量的常规振型组合就能得到相同精度的动态相对位移。

2附加振型的特性分析

2.1　单自由度结构的振型特性分析

先以一单自由度结构为例进行振型特性分析，其计算模型如图2所示。

图2　弹簧-质量模型图 Fig.2 Spring-mass model

图2中，模型A为支座固定的单自由度弹簧-质量体系；模型B取消了支座约束，在支座处附加大质量块M(取m、k的105以上[9])，同时释放支座水平方向的自由度约束；模型C在模型B的基础上取消了质量块m的质量。

三种模型下弹簧-质量模型的自振频率分别为：

(2)模型B的质量和刚度矩阵分别为：

式中，频率ω下标a表示附加频率，余同。

(3)模型C的质量和刚度矩阵分别为：

通过对上述三种模型的自振频率分析可以看出：

(1)模型B由于在支座处增加了两个大质量块，故多出了两阶振型，本文称多出的两阶振型为附加振型。当大质量块M质量大于内部质量块m质量的105倍时，模型B的常规振型频率ω1与模型A的频率ω1基本一致，误差小于万分之一，这个误差在工程应用中是完全可以接受的。

(2)模型C由于取消了内部节点质量块，只有2阶附加振型，且附加振型频率与模型B的附加振型频率完全一致，均为低阶频率。

对于钢筋混凝土结构，阻尼比约为5%，钢结构约为2%～3%，由于结构的阻尼比较小，因此计算结构的自振频率时，可以不考虑阻尼比的影响。

2.2　多自由度结构的振型特性分析

假设一个空间多自由度结构体系，内部质点数为s,自由度为p，支座质点数为t,自由度为q，如图3所示，采用上述3种计算模型进行分析。

图3　多自由度结构模型 Fig.3 Structural model with multi-degree of freedom

图3中，模型A为常规结构体系模型，支座节点固定,内部质点质量为mi,自由度为p；模型B取消支座约束，在支座节点附加大质量块Mi,大质量块质点自由度为q，整个结构自由度为p+q；模型C在模型B的基础上取消了内部节点质量，只留下支座节点质量Mi，整个结构质点自由度为q。

对于具有n个质点自由度的结构，共有n阶振型，且振型间有以下关系：

(1)

(2)

式中，γj为第j阶振型的振型参与系数，xji为j振型i质点水平相对位移，mi为第i个质点质量，r为质点数。

式(1)是抗震设计中常见的基本公式[15]，式(2)为振型参与质量计算公式[11]，是式(1)的另一种表现形式，但式(2)具有明显物理意义。对于n个自由度结构体系，当选择所有振型(n个振型)参与组合时，振型参与质量为100%。

下面对模型A、B、C的振型特性进行分析：

(1)先对模型A进行分析：

A模型内部质点数为s，自由度数为p，结构共有p阶常规振型，由式(1)可知:

(3)

(2)再对模型C进行分析：

模型C内部质点数为0，支座质点数为t，结构虽然有p+q个节点自由度，但只有q个质点自由度。由结构动力学相关知识可知:q个质点自由度的结构只有q阶振型。为了与后述模型B分析结果进行比较，振型编号为p+1,p+2,…，p+q，支座质点编号为要s+1,s+2,…，s+t。

由式(1)可知，对于模型C：

(4)

(3)最后对B模型进行分析：

B模型内部质点数为s，内部质点自由度为p，支座质点数为t，支座质点自由度为q，故结构体系共有p+q阶振型，由式(1)可得：

(5)

将上式振型参与系数具体化并拆分后可写成：

(6)

当附加质量Mi远大于内部节点质量mi时，模型B内部节点的质量mi可近似忽略，模型B与模型C的附加振型形状应基本相同，所以式(6)可写为：

(7)

将式(4)代入式(7)右边第二项可得：

(8)

式中，

将式(8)移项整理可得：

(9)

要使式(9)恒成立，只有：

(10)

从式(10)可以看出：模型B前p阶常规振型中，支座节点处位移恒为零，这与支座固定模型A的常规振型一致。

从上面分析可以得出如下结论：当大质量块的质量远大于内部质量点的质量时，模型B的附加振型与模型C的附加振型基本一致；模型B的常规振型与模型A的常规振型基本一致。

2.3　振型参与系数分析

文献[7]在式(18)、式(23)已经推导了随机地震动激励下虚拟动态位移的计算公式，从文献[7]中上述两式可以看出：与振型分解反应谱法一样，结构虚拟动态位移主要由结构的振型参与系数和振型形状确定。

前面已经得出了模型B的常规振型中支座处位移为零，与模型A的常规振型基本一致；模型B的附加振型与模型C的附加振型基本一致，现对三种模型的振型参与系数进行对比分析。

从以上分析可以看出，模型B的常规振型与模型A的常规振型的振型参与系数基本一致；模型B的附加振型与模型C的附加振型的振型参与系数基本一致。

3基于附加振型的虚拟激励求解方法

绝对位移直接求解的虚拟激励法采用附加振型进行求解时，虚拟激励的构建与文献[1]基本一致，故不再赘述，现重点阐述附加振型的求解原理。

3.1　附加振型求解过程分析

为便于理解，现把图3中空间结构模型C抽象成图4所示等效模型图，由于模型C中内部节点没有质量，故可将模型C等效成大质量块间采用广义弹簧连接的结构体系，如图4所示。

图4　模型C等效模型图 Fig.4 Equivalent model of model C

(11)

由于模型C中不考虑内部节点质量,所以附加振型数量只需由支座节点附加大质量块的自由度数来确定。虽然也需要全部附加振型参与计算才能保证拟静位移的计算精度，但对于大型实际工程来说，附加大质量块数远小于结构内部质量点数(比如大跨度桥梁，支座附加大质量块自由度数远远小于内部质点自由度数)，相应的附加振型的数量远小于常规振型的数量，所以采用附加振型来计算拟静位移，计算效率远高于Wilson位移输入模式。

3.2　振型贡献

基于附加振型的绝对位移直接求解虚拟激励法采用模型B进行求解，模型B与模型A的常规振型形状、振型参与系数基本一致；模型B与模型C的附加振型形状、振型参与系数基本一致。

由于模型A的常规振型主要贡献结构动态相对位移，模型C的附加振型主要贡献结构拟静态位移，所以模型B的常规、附加振型组合在一起可贡献结构的绝对位移。

相对运动法采用常规振型来贡献动态相对位移，采用静力法求解拟静位移；Wilson位移输入模式采用常规振型来贡献绝对位移；而基于附加振型的绝对位移直接求解法采用常规振型来贡献动态相对位移，采用附加振型来贡献拟静位移，这是上述算法本质区别。

3.3　奇异刚度矩阵求解方法

基于附加振型的振型分解法在支座增加了大质量块，且释放了支座约束，故整体结构振型分析时多出了一些附加振型，且这些附加振型中有一些是刚体运动振型，这种情况在旋转轴的扭转振动，飞机或导弹的总体振动中很常见。在这一类结构振动分析时，总体刚度矩阵是奇异的或病态的，这给数值计算带来困难，解决上述问题的最有效的方法就是特征值平移法，动力学教材中大多有该方法介绍，可参见文献[16]。新版本ANSYS软件在稀疏矩阵求解器中嵌入了该算法，能精确求解含病态矩阵的结构特征值和振型。

4基于附加振型求解的误差分析

(12)

由于振型分解法是将多自由度体系结构解耦成互不耦联的单自由度体系进行求解，对于单自由度体系，Ksb=-Kss，Csb=-Css，式中，Css=2ξωnM，所以式(12)可写为：

(13)

(14)

(15)

从式(15)可以看出，基于附加振型绝对位移求解法与相对运动法之间不存在由于阻尼假定不同而引起的阻尼误差。

5算例论证

算例采用一大跨度弦支穹顶，结构模型如图5所示，跨度为120 m，采用ANSYS建模，上部网壳如图5(a)所示，环向一共16圈杆件，下部预应力拉索由外到内隔圈布置，一共布置4圈环向预应力拉索(图5(b))。网壳杆件单元采用Beam188模拟，撑杆与径向拉杆采用Link8单元模拟，拉索采用Link10单元模拟，更多相关参数详文献[17]。分别采用计算模型A(固定支座模型)、计算模型B(支座附加大质量块模型)、计算模型C(支座附加大质量块，取消内部节点质量模型)三种计算模型进行分析。

图5　大跨度弦支穹顶模型 Fig.5 Models of large-span suspend-dome

本算例为多支座结构，但行波地震激励分析时并不需要在每个支座节点都附加大质量块，而是将相邻的支座节点先分为16组，每组节点自由度进行耦合，然后在每组增加一个可水平双向运动的大质量块，故模型B、模型C的附加振型都为32阶，三种力学计算模型下的结构自振频率见表1，表1给出了三种计算模型下的结构前8阶附加振型频率和前10阶常规振型频率，式中fa表示附加频率，从表1可以看出，模型A只有常规振型频率，模型B同时拥有附加振型与常规振型频率，且模型A的常规振型频率与模型B的常规振型频率基本一致；模型C由于取消了内部节点质量，故只有附加振型，且模型B的附加振型频率与模型C的附加振型频率基本一致，这与前面的理论分析完全吻合，证明了前面理论分析的正确性。

表1　弦支穹顶自振频率

图6为3种计算模型下的弦支穹顶结构振型对比，从图6可以看出，模型A的前4阶常规振型与模型B的前四阶常规振型基本一致，模型B的前4阶附加振型与模型C的前4阶附加振型基本一致，模型C由于内部节点无质量，故在振型图上没有显示出灰色的小质量块点。

采用文中论述的基于附加振型的虚拟激励法，对模型B进行行波地震激励分析，视波速取600m/s，结构阻尼比取5%；然后采用ANSYS内嵌的传统随机振动分析模块对模型A进行行波激励分析，最后将上述两种方法计算结果进行对比分析，需说明的是，ANSYS内嵌的传统随机振动分析模块采用相对运动法求解，其计算精度等同于林家浩教授提出的传统虚拟激励法，只不过是计算效率很低。

图6　3种模型下的弦支穹顶结构振型对比 Fig.6 Mode comparison of suspend dome with three dynamic models

图7　白噪声行波激励下单元剪力、位移响应功率谱密度对比 Fig.7 Comparison of shear and displacement power spectrum density under traveling-wave excitation based on white noise model

在对A模型采用传统随机振动进行计算时，内力、位移响应取前38阶常规振型的组合(38阶常规振型参与质量能达95%以上)；传统随机振动方法求解过程在文献[1]已进行了详细的阐述，本文不再赘述。在采用基于附加振型虚拟激励法对B模型进行计算时，拟静力内力、位移响应取前32阶附加振型的组合；动态内力、位移响应取前38阶常规振型的组合，绝对内力、位移响应取前32阶附加振型和前38阶常规振型的组合，两种算法下的单元内力和节点位移计算结果比较见图7。

随机抽取了15根杆件内力和节点位移进行了对比，发现两种算法计算结果基本一致，图7给出了2根杆件内力和2个节点的位移对比，图7(a)为E1单元动态轴力响应功率谱密度对比；图7(b)为E2单元绝对轴力响应功率谱密度对比；图7(c)为N1节点绝对位移响应功率谱密度对比；图7(d)为N2节点动态位移响应功率谱密度对比。从图7可以看出，除图7(b)中E2单元的绝对轴力响应功率谱在0.5Hz附近略有差别外，其余图形的结构响应在两种计算方法下的计算结果基本一致。从图7还可以看出，大跨度弦支穹顶的结构内力或位移响应并不是完全由低频控制，有时候高阶振型的贡献更大，不可忽略。

6结论

(1)绝对位移直接求解的虚拟激励法可以采用完全法进行求解，也可采用基于附加振型的振型分解法进行求解，采用完全法求解通常采用瑞利阻尼，采用振型分解法时通常采用振型阻尼。

(2)采用完全法进行求解时，计算效率相对较低，且瑞利阻尼模型下，当结构阻尼较大,或当激励频率较低时，与传统虚拟激励法相比，两种算法存在一定的质量阻尼误差[7-8]。

(3)采用基于附加振型的虚拟激励法求解时，构建少量附加振型组合就可精确获得拟静位移，选择与传统虚拟激励法相同数量的常规振型组合就能获得与之相同精度的动态位移。

(4)与传统虚拟激励法相比，基于附加振型的虚拟激励法大大方便了工程应用，且只需求解一次，不存在阻尼误差，计算效率高于传统的虚拟激励算法，计算精度与传统虚拟激励法相当。

参考文献

[1]李永华，李思明. 绝对位移直接求解的虚拟激励法[J]. 振动与冲击, 2009, 28(10): 185-190.

LI Yong-hua, LI Si-ming. Pseudo excitation method based on solving absolute displacement [J]. Journal of Vibration and Shock, 2009, 28(10): 185-190.

[2]张婷, 陈清军, 杨永胜. 两种改进虚拟激励法比较及大跨结构的行波效应分析[J]. 力学季刊，2012, 33(4):612-619.

ZHANG Ting, CHEN Qing-jun, YANG Yong-sheng.Comparison between two improved pseudo-excitation methods and analysis of long-span structural traveling waves effects. [J]. Chinese Quarterly of Mechanics，2012, 33(4): 612-619.

[3]赵银庆. 基于虚拟激励的车体结构随机疲劳问题研究[D].大连:大连理工大学, 2010.

[4]贾宏宇，郑史雄. 直接求解多维多点地震动方程的虚拟激励法[J]. 工程力学, 2013, 30(3):341-346.

JIA Hong-yu, ZHENG Shi-xiong. Pseudo excitation method of direct solving ground motion equation of multi-dimensional and multi-support excitation. [J]. Engineering Mechanics, 2013, 30(3):341-346.

[5]Wilson E L. Static and dynamic analysis of structures: A physical approach with emphasis on earthquake engineering [M]. Berkley, California: Computers and Structures, Inc, 2004.

[6]胡聿贤. 地震工程学[M]. 北京:地震出版社, 1988.

[7]李永华，桂国庆.相对运动法与绝对位移直接求解法算法误差分析[J]. 土木建筑与环境工程, 2011, 33(5):83-89.

LI Yong-hua, GUI Guo-qing. Error analysis between relative motion solving method and absolute displacement direct solving method[J]. Journal of Civil, Architectural & Environmental Engineering, 2011, 33(5):83-89.

[8]李永华, 卓平山, 胡远辉. 多点激励下绝对位移直接求解算法的误差频域分析[J]. 振动与冲击, 2012, 31(12): 112-119.

LI Yong-hua, ZHUO Ping-shan, HU Yuan-hue. Error analysis of seismic response obtained bydirectly solving absolute displacement in frequency domain under multi-support earthquake excitations[J]. Journal of Vibration and Shock, 2012, 31(12): 112-119.

[9]张林, 高翠琢, 吴景峰. 基于大质量法的谐响应有限元分析[J]. 半导体技术, 2011, 36(5): 402-405.

ZHANG Lin, GAO Cui-zhuo, WU Jing-feng. Harmonic finite element analysis based on large mass method[J]. Semiconductor Technology, 2011, 36(5): 402-405.

[10]林家浩, 张亚辉. 随机振动的虚拟激励法[M]. 北京:科学出版社, 2004.

[11]徐培福，傅学怡，王翠坤. 复杂高层建筑结构设计[M]. 北京：中国建筑工业出版社，2005.

[12]俞载道. 结构动力学基础[M]. 上海:同济大学出版社, 1987.

[13]同济大学数学教研室. 线性代数[M]. 上海:高等教育出版社, 1996.

[14]石永久,江洋,王元清.直接求解法在结构多点输入地震响应计算中的应用与改进[J].工程力学, 2011, 28(1): 75-81.

SHI Yong-jiu, JIANG Yang, WANG Yuan-qing. Application and improvement of direct solving method in seismic response analysis of structures under multi-support excitations [J]. Engineering Mechanics, 2011, 28(1): 75-81.

[15]丰定国，王社良. 抗震结构设计[M]. 武汉:武汉工业大学出版社, 2001.

[16]赵光恒. 结构动力学[M]. 北京:中国水利水电出版社, 1993.

[17]周海峰. 行波地震作用下大跨度弦支穹顶结构动力稳定性分析[D]. 南昌: 南昌大学, 2012.