杨秋伟 张金鑫 李翠红

∗(绍兴文理学院土木工程系，浙江绍兴312000)

†(宁波工程学院建筑与交通工程学院，浙江宁波315211)

∗∗(宁波工程学院浙江省土木工程工业化建造工程技术研究中心，浙江宁波315211)

有限单元法是进行大型复杂结构力学性能分析的有效手段[1-2]。然而，随着有限元模型自由度数目的增大，对其进行静力或动力特性计算所需要的计算成本也呈几何级数增长，对于规模巨大的结构，采用常规的电脑设备甚至完全无法进行相关计算。因此，对结构的有限元模型进行缩聚，以降低计算规模，并获得和原始模型尽量接近的计算结果，是工程中的必然要求。另一方面，工程实践中对结构进行静力或动力数据采集时，往往也只能获得部分自由度处的数据，因此，对结构的原始有限元模型进行缩聚，只保留和那些测试自由度相匹配的模型规模，这也是建立实验分析和理论分析的中间环节。因此，模型缩聚法在结构的静力和动力特性分析中有着广泛的应用。

Guyan[3]和Irons[4]最早研究了结构模型缩聚方法，并将其用于结构静力和动力计算中。他们的方法忽略了所有的惯性量，因此只能获得静力问题的准确解，对于动力问题则误差很大。O’Callahan[5]在Guyan 方法的基础上，在缩聚模型中增加了一阶惯性量，提出了改进的缩聚方法(improved reduced system，IRS)，有效提高了结构振动低价特征对的计算精度。Friswell 等[6-7]进一步改进了IRS 方法，通过迭代运算来逐步提高缩聚模型的精度，提出了迭代的模型缩聚法(iterated IRS，IIRS)。Qu 等[8-9]提出了逆迭代动力缩聚法，比IIRS 方法显著提高了计算效率。Kim 等[10]研究了评估Guyan 方法计算误差的估算方法。近年来，学者们将模型缩聚方法应用于解决结构非线性特征分析问题[11-14]。本文在IRS方法的基础上，提出了一种考虑二阶惯性量的改进IRS 方法，进一步改善了IRS 方法的计算精度。所提方法和迭代的模型缩聚法IIRS 相比，计算量更小且计算精度更高。以桁架结构和框架结构为例，验证了所提方法，并将计算结果与精确值、Guyan 缩聚解、IRS 缩聚解和IIRS 缩聚解进行了比较，结果说明了所提方法计算精度最好，具有良好的工程应用前景。

1 结构模型缩聚IRS 和IIRS 方法

本节对文献[5] 所提IRS 方法和文献[6-7] 中的IIRS 方法做简要回顾。IRS 方法利用级数展开，在模型缩聚方程中保留了一阶惯性量，从而比Guyan方法明显提高了计算精度。IIRS 方法在IRS 方法的基础上进一步通过逐次迭代来逼近精确解。主要方程回顾如下。

设结构有限元模型的刚度矩阵和质量矩阵分别用K和M表示，它们均为n×n维方阵，则结构自由振动的特征值λj和特征向量ϕj可以通过求解如下的广义特征值问题得到

当结构的自由度数目n很大时，直接求解方程(1)需要很大的计算量。为了节省计算花费，可以将结构有限元模型进行缩聚来降低矩阵维数，为此，首先将结构的自由度分为主要自由度(保留自由度)和次要自由度(缩减自由度) 两类，即将方程(1) 分解为

式中，上标m和s分别表示主要自由度和次要自由度，并且，m和s同时也表示主要自由度和次要自由度的数目。由方程(2) 的第二行可得

方程(3) 可以改写为

方程(4) 可进一步通过Neumann 级数展开为

工程中，为了简化计算，有限元模型中质量矩阵通常可以采用对角矩阵，则有

利用方程(6)，方程(5) 可以简化为

由方程(7)，如果忽略所有的惯性量(即含有质量子矩阵Mss的所有项)，则可得Guyan 缩聚模型

根据方程(8)，缩聚后的特征值方程如下

其中自由度转换矩阵T0为

显然，方程(9)中的缩聚刚度矩阵和质量矩阵的维数已降低为m×m维，一般而言m << n，因此和方程(1) 相比，计算量大大降低了，求解缩聚模型(即方程(9)) 可得到原模型(即方程(1)) 特征值和特征向量的近似解，这就是模型缩聚方法的基本原理。但由于忽略了级数展开式(9)中所有的惯性项，因此上述Guyan 缩聚法用于求解动力问题时精度较差。

为了进一步提高缩聚模型的精度，需要在级数展开式(7) 中保留部分惯性项，O’Callaghan 在文献[5] 中提出了IRS 方法，保留了第一阶惯性项，将方程(7) 转化为

方程(13) 中，λjϕmj可以由方程(9) 近似求得，即

将方程(14) 代入方程(13) 中整理可得

根据方程(15)，可得IRS 缩聚模型为

其中自由度转换矩阵T1为

为了进一步提高计算精度，Friswell 等[6-7]提出了IIRS 方法，其基本原理是，在方程(14) 中利用新的缩聚模型(即由方程(17) 和方程(18) 计算所得的新缩聚系统矩阵Kr1和Mr1) 来获得更准确的λjϕmj，重复由方程(14) 到(19) 的过程再次获得改进的缩聚模型，多次进行上述迭代过程，即可逐步逼近特征对的精确解。IIRS 方法自由度转换矩阵Ti的迭代计算公式为

2 模型缩聚改进IRS 方法

文献[15] 中讨论了IIRS 方法收敛的充分条件，从数学理论上证明了，只有当主辅自由度的选择满足收敛的充分条件要求时，才能保证迭代结果收敛到理论上的精确解，文献[15]中的算例表明，只有第一阶特征值满足此收敛条件，而对于其他特征值未能满足，且后期迭代对于精度的改善效果非常有限。根据文献[16] 所提供的计算结果，IIRS 方法并未获得和精确值一样的结果，后继即使增加迭代次数，计算精度并不能得到明显提升。因此，为了进一步提高IRS 方法的精度和计算效率，本文不再采用IIRS 方法那样的迭代操作，而是通过增加二阶惯性量来进一步提高计算精度。具体公式和原理详述如下。

在保留一阶惯性量的基础上，在级数展开式(7)中进一步保留第二阶惯性项，从而可将方程(7) 转化为

方程(22) 中，λjϕmj可以由方程(14) 近似求得，而可以由方程(14) 递推求得

将方程(14) 和方程(23) 代入方程(22) 中整理可得

根据方程(24)，可得二阶IRS 方法缩聚模型如下

其中自由度转换矩阵为

如果继续考虑三阶惯性量，采用类似的推导过程，可得三阶IRS 方法缩聚模型

其中自由度转换矩阵为

需要特别指出的是，由于所提方法中的惯性量是由Guyan 缩聚模型近似得到的，而Guyan 缩聚模型计算所得的特征值一般都是大于或者远大于精确值的，因此，所保留的惯性量阶数并非越大越好，后继的算例中会发现，如果保留了前四阶惯性量，其计算结果精度反而不如只保留二阶惯性量时的计算精度，这是因为随着特征值幂次的增高，惯性量的误差迅速增大，综合比较而言，考虑二阶惯性量的缩聚模型(即二阶IRS) 普适性最好。

3 数值算例

3.1 桁架结构

以图1 所示的桁架结构为例对所提改进IRS 方法进行验证。该结构的基本物理参数为：弹性模量E=200 GPa，密度ρ=7.8×103kg/m3，杆的长度L均为1 m，杆件的横截面面积A=1.759×10−4m2。除去支座约束，该结构共23 个自由度。不失一般性，取节点1、3 和5 所对应的6 个自由度为主要自由度，而其他节点自由度为次要自由度，主要自由度和全部自由度之比约为1/4。采用本文所提方法，分别考虑二、三和四阶惯性量，计算所得的前6 个特征值列于表1 中，为了便于进行比较，表1 中同时给出了精确值、Guyan 缩聚解、IRS 缩聚解和IIRS 缩聚解，其中精确值为采用完整有限元模型计算所得结果。

图1 桁架结构

由表1 可见，和精确值相比，Guyan 方法精度较差，说明Guyan 缩聚法不太适用于解决动力特征值计算问题。采用IRS 和IIRS 方法，可以获得比Guyan 方法更准确的计算结果。而采用本文所提方法，二阶IRS 计算结果精度最高，三阶和四阶IRS 计算结果反而精度变差。和IRS、IIRS 结果相比，也是二阶IRS 精度最好。另外，从计算时间上比较，IRS方法用时0.016 s，IIRS 方法用时0.032 s，本文所提二阶IRS 方法用时0.021 s，本文方法的计算时间介于IRS 和IIRS 之间。因此，综合而言，选择考虑二阶惯性量的缩聚模型(二阶IRS) 最佳。

3.2 框架结构

以图2 所示的框架结构为例对所提二阶IRS 方法再次验证。该结构的基本物理参数为：弹性模量E=69 GPa，密度ρ=2.7×103kg/m3，单元长度L均为0.2 m，单元横截面面积A= 6×10−5m2，惯性矩I=1.25×10−10m4。除去支座约束，该结构共66 个自由度。选取图2 中黑色圆点位置处的7 个自由度为主要自由度(节点2/5/12/16/20/23 均取水平自由度，节点14 取竖直自由度)，剩余其他节点自由度均视为次要自由度，主要自由度和全部自由度之比约为0.1。采用本文所提方法，考虑二阶惯性量，计算所得的前7 个特征值列于表2 中，表2 中同时给出了精确值、Guyan 缩聚解、IRS 缩聚解和IIRS 缩聚解。

表1 前6 阶特征值计算结果(算例1)

图2 框架结构

由表2 可见，和精确值相比，采用IRS 方法，前4 个特征值的计算误差均在10%以内，但第7 个特征值完全失真。而采用IIRS 方法，第7 个特征值的计算误差有明显缩小，但其他特征值的误差和IRS 相比没有明显优势。采用本文所提二阶IRS 方法，前6 个特征值的计算误差均控制在5%以内，第7 个特征值的误差也显著降低为45.77%。从计算时间上比较，IRS 方法用时0.018 s，IIRS 方法用时0.035 s，本文所提二阶IRS 方法用时0.026 s，本文方法的计算时间介于IRS 和IIRS 之间。综合而言，本文方法计算精度更加可靠，实际应用时可将本文方法与IRS方法结合使用，互相对比验证，以增强计算结果的可靠性。

表2 前7 阶特征值计算结果(算例2)

4 结论

模型缩聚方法可以有效降低结构的有限元计算规模，节省计算时间和成本，并能获得和实验测量自由度相匹配的有限元模型。本文在IRS 方法的基础上，进一步考虑二阶惯性量，明显提高了计算精度。本文方法和迭代缩聚IIRS 方法相比，计算量更小且精度更高。以桁架结构和框架结构为例对所提二阶IRS 方法进行了验证，并将计算结果与精确值、Guyan 缩聚解、IRS 缩聚解和IIRS 缩聚解进行了比较，结果表明：综合而言，所提方法计算精度最好，具有良好的工程应用前景。