基于有限元Craig缩聚模型的基础激励和载荷识别

朱斯岩

(北京航空航天大学 航空科学与工程学院，北京 100191)

朱礼文

(北京宇航系统工程研究所，北京 100076)

基于有限元Craig缩聚模型的基础激励和载荷识别

朱斯岩

(北京航空航天大学 航空科学与工程学院，北京 100191)

朱礼文

(北京宇航系统工程研究所，北京 100076)

载荷识别是用有限自由度的加速度测量值推断外力.基础激励是在基础自由度上加给定的加速度时间历程激振结构，求整个结构的动响应.利用Craig减缩模型的性质，及界面自由度的已知运动可求出整个结构其余自由度的动力学响应，完成基础激励求解.在此基础上用已知的界面自由运动和求得的非界面自由度的动力学响应，代入结构动力学方程即可实现载荷识别.本方法唯一能给基础激励和载荷识别结果带来的偏差是Craig缩聚模型与原始有限元模型间的差异，对此用频响分析的方式，直接从频响曲线的变化来确定缩聚模型的精度.算例表明这一新的基础激励和载荷识别途径是有效的.

结构动力学;基础激励;载荷识别

基础激励和载荷识别是工程应用中常见的二类结构动力学问题，其共同点是均用有限自由度的加速度时间历程，实现基础激励或载荷识别的目的.基础激励类似响应计算，二者的区别在于，响应计算中动力学方程的外力函数项已知，问题转化为二阶微分方程(组)的数值积分，目前已有很多成熟精确的算法［1－3］;基础激励是用结构少数基础自由度的加速度时间历程，求结构其余自由度的响应［4－5］.文献［4］介绍了 3 种方法:惯性释放法，适用于基础为1个节点6个自由度的情况;大质量法，是一种近似算法;Lagrange乘子法(LMT)，较为繁琐.载荷识别由于种种实际因素限制，只能用少数自由度的加速度响应测量值推断作用于结构的激振力F(t)，而对其余自由度响应所做的不同假设，都将导致载荷识别结果的不同，即通常所说的不适定性［6］，也是研究的难点.本文利用Craig减缩模型的性质，由已知的界面自由度加速度时间历程，可精确简捷地求得全部非界面自由度的动响应，完成基础激励的求解.同时，载荷识别就可以先按上述方式，用测得的界面自由度加速度时间历程，求非界面自由度动力学响应，再将所有界面和非界面自由度的动响应代入结构动力学方程，即可实现载荷识别，获得作用于结构界面的外载荷.

以上方法唯一会给基础激励和载荷识别结果带来的偏差，是Craig缩聚模型与原始有限元模型间的差异.有限元模型的Craig缩聚是固定界面模态综合的基础［7］，最初它与静力学子结构方法一样，只是为了适应当时计算机速度低、容量小的硬件条件.1987年长征系列运载火箭开始卫星发射的国际商业服务时，在星箭耦合载荷计算中运用了这一技术.而此时构造的捆绑火箭结构动力学数学模型，同样采用了固定界面模态综合方法［8］.可见Craig缩聚模型已十分广泛地应用于航空航天结构动力学分析中，因此对Craig缩聚模型与原始有限元模型间的差异进行定量分析是十分必要的.最早王文亮教授在文献［9］中对此作了有效性分析，而普遍采用的是主模态截止频率，认为缩聚模型与原始模型在主模态截止频率以下的模态子空间相当.此外文献［10］运用模态有效质量(modal effectivemass)的概念.本文提出用频响分析的方式，直接从频响曲线的变化来确定缩聚模型的精度.

1 Craig模型缩聚

有限元模型的Craig缩聚，要求把运动参数分为界面自由度和非界面自由度，此时运动方程变为

本文将有外力作用的自由度列为界面自由度，界面模态φc定义为

由式(2)可见，界面自由度的选取应保证矩阵Kmm可逆，因此必要时须选取无外力作用的自由度为界面自由度，并认为其上作用的外力等于0.对于其余的全部自由度均归为非界面自由度，解如下广义特征值问题:

得到主模态φm和主模态特征值λm.

由主模态φm和界面模态φc构造变换矩阵T:

则运动参数的物理空间和模态空间变换关系为

由式(5)可知界面自由度保留了物理空间的性质.

缩聚后式(1)变换为

其中

刚度矩阵的非对角块，与主模态坐标及界面自由度相关的项为0，即

运动方程右侧非界面自由度对应外力为0，坐标变换后对应主模态坐标的广义力仍为0.同时与主模态对应的质量阵和刚度阵均为对角阵，已解耦，对应阻尼项可以用模态阻尼比方式引入.通常界面自由度数较少，可不计相应的阻尼项.故列入阻尼阵后式(6)变为

式(8)可分解为如下二式

模型缩聚是依靠减少主模态数实现的.若算例规模较小，允许计算出与非界面自由度数相等的主模态数，此时通过变换矩阵式(4)生成的结构模型即Craig变换模型(模型大小与原始模型相同).此时，用Craig变换模型计算的特征值、特征矢量，与原始模型计算得到的特征值、特征矢量完全一致.

界面自由度的加速度值，在基础激励问题中是已知的激励项，可据此求整个结构的动响应;而在载荷识别问题中是已知的测量值，需要用它推断作用于结构上的外力.

2 基础激励

基础激励是用加速度时间历程激振结构基础，因此无法像响应计算那样，用外力时间历程通过积分结构动力学方程获得结构的动响应.NASTRAN早期提供了“大质量法”的基础激励计算模块，近期的版本包含了直接用加速度时间历程的计算模块.

本文基础激励问题的求解运用式(9a).由于主模态坐标对应的广义力为0，因此将已知项移至等号右侧并变号，得到响应计算的方程形式如下:

3 载荷识别

载荷识别可依据Craig减缩模型，先按照第2节基础激励的方式，用有限的界面自由度的测量值计算主模态坐标参数的动响应，再根据式(9b)完成载荷识别计算.式中，主模态坐标加速度按基础激励方式求得为已知的加速度测量值;Xc为界面自由度的位移时间历程，是未知项，可由加速度测量值数值积分二次求出，积分时初值可简单取0，同时加速度测量值起点也定为0.初值对计算结果有影响，但因阻尼的存在将只影响计算起始的很小一段时间［11］.将代入式(9b)便能得到通过界面自由度作用于结构的载荷.

上述载荷识别的计算方法同样适用于模拟振动台运动的基础激励分析.由于振动台作为基础只包含1个节点6个自由度，是静定约束.相应的界面模态是对应结构的刚体运动，因此与刚度矩阵相乘为0，故式(9b)中矩阵，可直接根据计算振动台对试件的激振力和夹持力.显然这对振动台反共振问题分析具有重要价值.

4 Craig缩聚模型精度评估

频响分析是在结构的给定自由度上，加简谐形式的单位载荷，计算结构指定自由度的位移(速度、加速度)响应随加载频率变化的规律.计算得到的位移(速度、加速度)频响曲线，是结构固有特性(传递特性).目前所有有限元商用软件都有频响分析模块.

结构动力学方程为

其中，响应U和激励F是简谐函数，频率为ω.

通常响应和激励间存在相位差.将式(12)、

式(13)代入式(11)，得到频响分析运动方程:

将式(14)分成实部、虚部:

通过式(14)或式(15)求解U1和U2，分别由U1，U2对应的自由度处取值，求该自由度的位移响应、相位差，再乘以ω，ω2，即得到速度、加速度响应.

用Craig缩聚模型作频响分析时，式(14)、式(15)中的质量阵、刚度阵和阻尼阵可直接由式(8)得到，但激励自由度必须列为界面自由度，响应自由度建议也列入界面自由度.否则必须将式(14)或式(15)中的 U1和 U2分别取出，用生成Craig缩聚模型的变换阵T，按式(5)恢复到物理空间，再由相应自由度处取值求该自由度的位移响应、相位差，再乘以ω或ω2，获得速度响应、加速度响应.

5 算例

算例是长10m，直径500mm，壁厚3mm的铝质圆管.数学模型包含11个节点，共66个自由度，由10个等长度的铁木辛可梁单元组成，不加约束.每个节点集中质量900 kg，集中转动惯量0.3 kg·m2(绕圆管中心线转动惯量为 0.6 kg·m2).选取 10 个界面自由度，其序号为:1，2，3，4(节点1);25，26，27(节点5);49，50，51(节点9).频响计算以自由度序号50作为加载点，自由度2的加速度频率响应曲线如图1所示.缩聚模型的主模态阶数分别取前8，10，24阶3种情况，界面自由度数为10，因此缩聚模型总自由度数分别为18，20，34.

图1 加速度频率响应曲线

从图1曲线可以看出主模态阶数变化时对应曲线峰值的位置和高度都会发生变化.实际上，该算例当主模态阶数大于等于18后，主模态阶数的增加不再引起频响曲线计算结果的变化.此外，对于根部固支或其他多种约束方式，结构系统的频响曲线也大体类同.

载荷识别算例也采用了同一结构模型，并在界面自由度2，6(即节点1的Y向和节点5的Y向)加激振力，计算结构动响应.其中界面自由度的加速度时间历程作为载荷识别的加速度测量值.同时为方便起见，两个力函数仅在时序上颠倒了一下，其余界面自由度上外力函数皆为0.计算结果显示识别载荷同激振力十分一致.为突出二者差别，在确定界面自由度运动(作为测量值)的响应计算中，结构模型采用主模态阶数等于非界面自由度数的Craig变换模型，而载荷识别时采用主模态阶数等于8的Craig缩聚模型.因图1的频响曲线显示这两个模型的差异主要表现在40Hz以上的频段，因此激振力也应有40Hz以上的频率成分.激振力的幅频特性如图2所示.图3给出了自由度2上识别载荷与激振力的比较，二者十分一致.

图2 原始所加载荷傅里叶幅值谱

图3 原始所加载荷与载荷识别结果比较

6 结论

本文提出利用有限元Craig缩聚模型进行基础激励计算和载荷识别，并根据频响曲线确定模型精度，算例计算表明方法可行有效.

本文基础激励的计算方式在预示振动试验的分析工作中已实际应用，但载荷识别方法由于要求加力点必须同时是加速度测量点，这就要求事先在加速度测点布置中有所安排考虑，因此在现有工作中进行验证有一定困难.此外，外力作用处直接安装加速度测点也有可能存在难度，这时可以将测点设定在外力作用点附近、与着力点间连接刚性较好的结构上，相当于将力移动到测点.若该点处在力作用线上，则对随后的基础激励和载荷识别计算没有影响.若因外力移动出现附加弯矩，则对应的加速度测量不只是线加速度，还要求增加角加速度测量.由于目前尚无法直接测量角加速度，可以在外力作用点附近寻找刚性较大的结构，用多个线加速度测点测量结果换算该点的角加速度.总之应尽可能避免附加弯矩的出现.

References)

［1］钟万勰.结构动力方程的精细时程积分法［J］.大连理工大学学报，1994，34(2):131 －135 Zhong Wanxie.On precise time-integration method for structural dynamics［J］.Journal of Dalian University of Technology，1994，34(2):131－135(in Chinese)

［2］ Bathe K J.Finite element procedures in engineering alalysis［M］.Englewood Cliffs，New Jersey:Prentice-Hall Inc，1982:499－532

［3］邱吉宝，向树红，张正平.计算结构动力学［M］.合肥:中国科学技术大学出版社，2009:260－274 Qiu Jibao，Xiang Shuhong，Zhang Zhengping.Computational structural dynamics［M］.Heifei:University of Science and Technology of China Press，2009:260 －274(in Chinese)

［4］董益亮，郭钢，徐宗俊，等.道路激励作用下的汽车后桥动力响应分析［J］.汽车工程，2002，24(4):339－342 Dong Yiliang，Guo Gang，Xu Zongjun，et al.The dynamic response analysis of the automotive rear axle under road excitation［J］.Automotive Engineering，2002，24(4):339 － 342(in Chinese)

［5］任建峰，仇原鹰，段宝岩，等.一种基础加速度冲击的仿真新方法［J］.现代机械，2006，132(2):21 －23 Ren Jianfeng，Qiu Yuanying，Duan Baoyan，et al.A new idea to access the simulation of a shock from the base［J］.Modern Machinery，2006，132(2):21 －23(in Chinese)

［6］王彦飞.反演问题的计算方法及其应用［M］.北京:高等教育出版社，2007:1－5 Wang Yanfei.Computational methods for inverse problems and their application［M］.Beijing:Higher Education Press，2007:1－5(in Chinese)

［7］ Craig R R Jr，Chang C C.Substructure coupling for analysis and testing［R］.NASA CR-2781，1977

［8］于海昌，贾文成，卫国，等.大型捆绑火箭模态试验/分析的相关性研究［J］.导弹与航天运载技术，1993，202(2):42 －53 Yu Haichang，JiaWencheng，WeiGuo，et al.A correlation investigation ofmodel test/analysis for heavy lift strap-on launch vehicle［J］.Missiles and Space Vehicles，1993，202(2):42 －53(in Chinese)

［9］王文亮，杜作润.结构振动与动态子结构方法［M］.上海:复旦大学出版社，1985:229－426 Wang Wenliang，Du Zuorun.Structural vibration and dynamic substructure methods［M］.Shanghai:Fudan University Press，1985:229－426(in Chinese)

［10］ López-Díez J，Marco-Gómez V，Luengo P.Modal effective mass applied to error localization of finite element models［R］.AIAA-99-1310，1999

［11］朱斯岩，朱礼文.运载火箭动态载荷识别研究［J］.振动工程学报，2008，21(2):135 －139 Zhu Siyan，Zhu Liwen.Dynamic load identification on launch vehicle［J］.Journal of Vibratioin Engineering，2008，21(2):135－139(in Chinese)

(编 辑:李 晶)

Base excitation and load identification based on Craig reduction model

Zhu Siyan

(School of Aeronautic Science and Engineering，Beijing University of Aeronautics and Astronautics，Beijing 100191，China)

Zhu Liwen

(Beijing Institute of Astronautical Systems Engineering，Beijing 100076，China)

Load identification use measured acceleration of few degree of freedom(DOF)to determine external applied load.Base excitation give movement upon foundation freedom to excite the structure，and then obtain the whole dynamic response.Making use of characteristics of Craig reduction model，and the known motion of interface DOF，the dynamic response of all other DOF could be calculated，thus the base excitation was achieved.The load identification could be performed through the dynamic equation with these motion of interface DOF and the calculated dynamic response of non-interface DOF.The only error of this method for base excitation and load identification is the difference between the Craig reduction model and the original finite-element model.For this reason，frequency response analysis was applied to determine the precision of Craig reduction model through the variation of the frequency response curve.An example was given to demonstrate the validity of the present base excitation and load identification method.

structural dynamics;base excitation;load identification

TB 122

A

1001-5965(2011)06-0665-04

2010-03-10

朱斯岩(1973 －)，女，北京人，讲师，ysz4444@buaa.edu.cn.