引 言
发动机安装节是连接飞机机体和发动机的重要结构部件．发动机产生的推力通过安装节、吊挂传递至机体；同时由于发动机低压和高压转子的不平衡，产生的振动也通过安装节、吊挂传递至机体结构，共同引发机体结构振动，从而在机舱内产生结构噪声．因此，飞机舱内噪声控制研究最重要的工作是识别并获取发动机安装节上发动机安装点的动态力，由此精确计算机体结构的振动以及由于发动机振动产生的舱内噪声，进而采用特定的方式抑制这类振动和噪声[1]．由此可知，精确识别发动机安装节上发动机安装点的动态力是振动与噪声控制工作中的一项重要工作．但是，由于力传感器尺寸和安装条件的限制，嵌入式力传感器必将改变原有系统的动力特性，测量结果常出现偏差[2]，而在飞行条件下更无法采用直接方法获取安装节上的动态力．
鉴于此，基于测量得到的结构振动响应，采用求逆运算的方式对安装节激励点的动态力进行计算，不仅在噪声与振动控制工程中具有重大意义，对结构动力学分析也具有重要应用价值[3]．目前，利用测量得到结构激励与响应点之间的FRF函数(frequency response function, 即频率响应函数)矩阵，采用最小二乘法对激励点动态力进行估计的方法在动态力识别研究中已被广泛应用，且该方法已被应用于解决特定的工程问题[4-5]．具体而言，该方法利用在结构上一个或多个点测量得到的振动响应，同时配合激励点至响应点之间的FRF函数，通过将测量得到的振动响应向量与FRF函数伪逆矩阵相乘来计算激励力．1979年，Bartlett等[6]利用最小二乘法，在直升机模型中确定了主桨毂上的垂直和侧向动态力，同时讨论了该方法在工程应用中的正确性和可行性．此外，Okubo等[7]将该方法应用于不同结构，成功识别到了这些结构上的动态力．
本文首先介绍采用最小二乘法识别结构动态力的基本理论，然后应用该方法对某型飞机飞行状态下安装节上发动机安装点的动态力进行识别．利用发动机安装点与吊挂之间的FRF函数矩阵，以及发动机安装点与安装节上部点的FRF函数矩阵，结合飞行时测量得到的吊挂和安装节上部点的振动加速度响应，对安装节上发动机安装点的载荷进行了两次计算．通过对比两次计算的矩阵条件数，确定由发动机安装点与吊挂之间的FRF函数矩阵计算得到的动态力相对更为准确，并对两次计算得到的动态力相对误差进行了对比．
1 最小二乘法动态力识别理论
1.1 最小二乘法动态力识别
最小二乘法是动态力识别研究中应用最广泛的方法．针对单自由度振动系统，振动响应y(x,t)与x0位置的激励力f(x0,t)可用时域振动方程表示为[5]
y(x,t)=
A(x|x0,t-τ)f(x0,τ)dτ，
(1)
式中，A表示时间域上动态力激励位置与振动响应测量位置之间的传递函数．对式(1)左右两侧进行Fourier变换，该方程在频域中可表达为
(2)
式中，
和
分别代表频率ω处的复变函数．
针对多自由度振动系统，式(2)可用矩阵形式表示为
(3)
其中m表示振动响应测量的数量，n表示需要被识别的动态力数量，
表示动态力激励位置xj与振动响应测量位置xi之间的FRF函数．对于动态力识别问题，式(3)中的
和是已知量．若m=n,可直接通过利用FRF函数逆矩阵
计算得到激励点的动态力．但是，在实际工程应用中，m与n通常不完全相等，同时在实际测量过程中还可能产生振动响应测量误差，故需采用最小二乘法，通过逆运算识别激励点动态力[8]．一般而言，利用最小二乘法，可通过增加振动响应测量点的方式来降低振动响应测量误差对激励点动态力识别精度的影响，从而有效提高激励点动态力识别的精度[8]．利用式(3)，可将识别的激励点动态力表示为
(4)
式中，
表示激励点动态力
表示动态力激励位置xj与振动响应测量位置xi之间的FRF函数
表示方阵
的伪逆矩阵,可表示为
(5)
其中，角标H表示矩阵的转置共轭运算．
结构激励点与响应点之间的FRF函数是结构的固有属性，与外界激励无关[9]．通常，可通过力锤激励法或激振器法对FRF函数进行测量[10]．测量得到FRF函数以及实际工作时的振动响应后，即可采用式(4)计算激励点的动态力．计算过程中，对FRF函数进行矩阵逆运算将影响激励点动态力识别的精度．实际应用中，由于结构在某些频率存在共振特性，导致FRF函数的伪逆矩阵可能出现病态，从而导致结构振动响应测量
中的微小误差被逐步放大，最终严重影响动态力的识别精度．鉴于此，需要利用逆运算过程中的矩阵条件数对计算过程进行评估：矩阵条件数越低，则识别的动态力精度越高．
1.2 矩阵条件数准则
以向量
表示振动响应测量误差．根据式(3)，包含误差的动态力识别可表示为
(6)
式中，
表示识别的动态力精确值，
表示包含振动响应测量误差
的振动响应．根据奇异值分解[11]，可将FRF函数矩阵分解为
(7)
式中，U和V为酉矩阵，VH为矩阵V的伴随矩阵，Σ为对角矩阵，其中对角矩阵中的主对角线元素为矩阵的奇异值．因此，对角矩阵Σ可表示为
(8)
式中，s1,s2,…,sn为按照降幂排列的非负奇异值，diag(s1,s2,…,sn)表示以这些非负奇异值为主对角元素的对角矩阵．结合式(4)，包含振动响应测量误差的动态力识别结果可表示为
(9)
式中，ui为矩阵U的列向量，vi为矩阵V中的第i阶，ηi为向量的第i阶，si为的奇异值矩阵中的第i阶．如果由于结构共振导致矩阵在某些频率处出现病态，则最大与最小奇异值的比值将非常大，从而导致微小的振动响应测量误差被放大[5,8]．工程应用中，一般采用FRF函数矩阵的条件数对动态力识别的准确度进行评估[12]．矩阵条件数可通过最大奇异值与最小奇异值的比值进行计算：
(10)
式中，
表示矩阵条件数．如果矩阵在某些频率处出现病态，则的值在该频率处会非常高[8]．工程应用中，对κth通常设定100作为其阈值[13]，如果大于100，可利用截断奇异值分解或Tikhonov正则化方法处理特定频率范围内的逆运算，使得该逆运算能够继续进行，实现矩阵求逆数值收敛．本文通过合理选择发动机安装节与吊挂之间FRF函数的测量位置以及降低振动响应的测量误差，使得计算频率范围内动态力识别的矩阵条件数的最大值未超过100，从而满足了实际工程中的应用要求．
2 发动机安装节与吊挂间的FRF函数
采用力锤激励法或激振器法对结构的FRF函数进行测量．对多自由度振动系统，FRF函数矩阵可表示为
(11)
式中,Gff(ω)为激励力的自谱矩阵，Gfx(ω)表示激励点x与振动响应点之间的互谱矩阵．目前，各类形式的数据采集系统对FRF函数测试均有标准的测量模板，根据测量得到的激励力信号和振动响应加速度信号，能够自动根据式(11)计算得到响应的FRF函数．本文采用力锤激励法，对发动机刚性安装节与吊挂之间的FRF函数矩阵进行测量．图1所示为发动机安装节和吊挂FRF函数测量的力锤激励点和振动加速度响应点．
图1 发动机安装节和吊挂FRF函数测量示意图
Fig. 1 The FRF measurement diagram of the engine mount and pylon
针对发动机安装节和吊挂结构的特征，本文测量了两组FRF函数矩阵．第一组为发动机安装节上/下部安装激励点与吊挂振动加速度响应点之间的FRF函数矩阵；第二组为发动机安装节上/下部安装激励点与安装节上振动加速度响应点之间的FRF函数矩阵．飞机处于静止状态下，利用力锤激励法完成上述两组FRF函数矩阵测试，图2～4给出了测量得到的第一组FRF函数，图5～7给出了第二组FRF函数．
图2 安装节下部安装激励点X方向的FRF函数 图3 安装节下部安装激励点Y方向的FRF函数
Fig. 2 FRFs at the lower point in the X-direction Fig. 3 FRFs at the lower point in the Y-direction
图4 安装节下部安装激励点Z方向的FRF函数 图5 安装节上部安装激励点X方向的FRF函数
Fig. 4 FRFs at the lower point in the Z-direction Fig. 5 FRFs at the upper point in the X-direction
图6 安装节上部安装激励点Y方向的FRF函数 图7 安装节上部安装激励点Z方向的FRF函数
Fig. 6 FRFs at the upper point in the Y-direction Fig. 7 FRFs at the upper point in the Z-direction
3 飞行状态下振动加速度响应测量
为了对飞行状态下发动机低压和高压转子引发的动态力进行计算，需要测量巡航状态下对应FRF函数振动响应点的振动加速度．
图8和图9为发动机安装节和吊挂上三向振动加速度计的安装图，与图1中FRF函数测量的振动响应点位置对应．
图8 安装节上三向振动加速度计安装 图9 吊挂梁上三向振动加速度计安装
Fig. 8 The tri-axial accelerometer on the mount Fig. 9 The tri-axial accelerometer on the pylon
测量时，飞机飞行高度为35 000 ft(1 ft=0.304 8 m)，速度为0.78Ma．安装节和吊挂梁上测量得到的振动加速度频谱如图10和图11所示．
图10 安装节振动加速度频谱
Fig. 10 Acceleration spectra of the engine mount
图11 吊挂梁振动加速度频谱
Fig. 11 Acceleration spectra of the pylon
4 发动机安装节动态力识别
根据本文中描述的最小二乘法，利用飞行状态下测量得到的振动加速度响应和对应的两组FRF函数矩阵，对发动机安装节上部安装激励点和下部安装激励点的动态力进行识别．图12、13为利用安装节上部安装点、下部安装点与吊挂振动加速度响应点之间的FRF函数矩阵以及吊挂振动加速度响应，识别到的安装节上的动态力．
图14、15为利用安装节上部安装点、下部安装点与安装节上振动加速度响应点之间的FRF函数矩阵以及安装节上振动加速度响应，识别到的安装节上的动态力．
通过对比图12～15，在识别的动态力频谱中可以发现，利用两组FRF函数矩阵和对应的振动加速度响应，都识别到了由发动机低压转子引起的90 Hz附近的动态力和由发动机高压转子引起的260 Hz附近的动态力，识别到的动态力数值相近．两组FRF函数矩阵识别动态力过程中的矩阵条件数如图16所示．
图12 发动机安装节下部安装点动态力频谱
Fig. 12 Identified force spectra at the lower point
图13 发动机安装节上部安装点动态力频谱
Fig. 13 Identified force spectra at the upper point
图14 发动机安装节下部安装点动态力频谱
Fig. 14 Identified force spectra at the lower point
图15 发动机安装节上部安装点动态力频谱
Fig. 15 Identified force spectra at the upper point
从图16可知，对于两组FRF函数矩阵动态力识别过程，在0～2 000 Hz范围内矩阵条件数均低于100．通过对比，利用安装节上部安装点、下部安装点与吊挂振动加速度响应点之间的FRF函数矩阵(第一组)以及吊挂振动加速度响应识别到的动态力对应的矩阵条件数在整个频率范围内都低于利用安装节上部安装点、下部安装点与安装节上振动加速度响应点之间的FRF函数矩阵(第二组)以及安装节上振动加速度响应识别到的动态力对应的矩阵条件数，即识别到的第一组动态力比第二组动态力更为精确．本文通过计算两组识别到的动态力RMS(root mean square,均方根)值，计算识别的动态力相对误差，如表1所示．
图16 矩阵条件数对比
Fig. 16 Comparison of norm condition numbers
表1 两组FRF函数矩阵发动机安装节动态力识别相对误差
Table 1 Relative errors of the identified dynamic forces for the 1st and 2nd FRF groups
parameterpositionlower excitation point on the mountfX/NfY/NfZ/Nupper excitation point on the mountfX/NfY/NfZ/Nidentified force RMS(the 1st FRF matrix)50.7652.9624.0126.80100.0222.08identified force RMS(the 2nd FRF matrix)38.8450.6128.1016.8480.3821.76relative error δ/%23.484.4417.0337.1613.3214.49
通过对比表1中两组FRF函数矩阵计算得到的发动机安装节动态力相对误差，可以得出，对于安装节上部X方向动态力，相对误差达到了37.16%，误差相对较大，但误差绝对值不大；其他位置和方向的动态力相对误差都相对较小．两组FRF函数矩阵计算得到的飞行状态下安装节动态力均符合工程应用要求．
在动态力的频谱域上，因为发动机最关注的低压转子动态力频率为80～90 Hz，高压转子动态力频率为260～280 Hz，因此主要对300 Hz以下频率范围的动态力误差进行研究．其中，安装节下部动态力误差如图17所示，安装节上部动态力误差如图18所示．
图17和图18表明，安装节上部Y方向动态力在260～280 Hz频率上误差波动较大，其余位置和方向上的误差都在动态力识别工程可接受的范围内．
图17 安装节下部动态力误差图 图18 安装节上部动态力误差图
Fig. 17 Relative errors of dynamic forces Fig. 18 Relative errors of dynamic forces at the lower point at the upper point
5 结 论
飞机发动机安装节是用于连接发动机与机体的重要结构部件．为了对发动机引起的机体结构振动和发动机引发的机体结构振动噪声进行预测和抑制，需要获取发动机安装节上发动机安装点的动态力．但是，在飞行状态下，该动态力无法通过直接测量的方式获取，所以本文详细描述了利用最小二乘法对动态力识别的方法．本文结合实际的工程应用，选取发动机刚性安装点、吊挂与机身结构等具有响应线性关系位置进行测试，利用两组发动机安装节安装点与结构振动响应点之间的FRF函数以及飞行状态下的振动加速度响应，成功获取了飞行状态下发动机安装节上部安装点和下部安装点的动态力．同时，利用矩阵条件数准则，对识别到的两组动态力进行了对比分析，并得到了以下结论：
1) 基于最小二乘法的动态力识别方法能成功用于飞机发动机安装节飞行状态动态力识别，获取了飞行状态下发动机工作在安装节产生的动态力．
2) 动态力识别属于逆向问题，本文利用矩阵条件数对动态力识别过程进行评估，利用矩阵条件数来判定识别到的动态力的精确程度．
3) 在实际工程应用中，由于飞行状态下绝对精确的动态力无法获得，因此本文通过使用两组FRF函数矩阵及对应的振动加速度响应，对同一飞行工况下的发动机安装节动态力进行估计，并对两组识别到的发动机安装节动态力进行相对误差计算．计算结果表明，对于安装节上部X方向动态力，相对误差达到了37.16%，但动态力误差绝对值不大，对于其他位置和方向，动态力相对误差都较小，满足工程应用要求．
4) 本文详细叙述了基于最小二乘法的动态力识别过程，其他类似的工程问题也可参考本文进行应用，并利用矩阵条件数对识别到的动态力进行评估．
[1] 杜建镔. 结构优化及其在振动和声学设计中的应用[M]. 北京: 清华大学出版社, 2015.(DU Jianbin. Structural Optimization and Its Application in Vibration and Acoustic Design[M]. Beijing: Tsinghua University Press, 2015.(in Chinese))
[2] 张磊, 曹跃云, 杨自春, 等. 总体最小二乘正则化算法的载荷识别[J]. 振动与冲击, 2014, 33(9): 159-164.(ZHANG Lei, CAO Yueyun, YANG Zichun, et al. Load identification using CG-TLS regularization algorithm[J]. Journal of Vibration and Shock, 2014, 33(9): 159-164.(in Chinese))
[3] 姚学锋, 杨桂, 姚振汉, 等. 先进复合材料自行车架的动力学特性分析[C]//第六届全国结构工程学术会议论文集. 北京, 1997.(YAO Xuefeng, YANG Gui, YAO Zhenhan, et al. Dynamic characteristics of advanced composite bicycle frame[C]//The 6th National Academic Conference on Structural Engineering. Beijing, 1997.(in Chinese))
[4] 陆秋海, 李连友, 向律楷, 等. 非平稳环境激励下结构工作模态参数识别法[J]. 清华大学学报(自然科学版), 2013, 53(3): 389-393.(LU Qiuhai, LI Lianyou, XIANG Lükai, et al. Operational modal parameter identification of structures for non-stationary ambient excitation[J]. Journal of Tsinghua University (Science and Technology), 2013, 53(3): 389-393.(in Chinese))
[5] XU X, OU J P. Force identification of dynamic systems using virtual work principle[J]. Journal of Sound and Vibration, 2014, 337: 71-94.
[6] BARTLETT F D, FLANNELLY W G. Model verification of force determination for measuring vibratory loads[J]. Journal of the American Helicopter Society, 1979, 24(2): 10-18.
[7] OKUBO N S, TATSUNO T. Identification of forces generated by a machine under operating condition[C]//Proceedings of International Modal Analysis Conference. Chicago, 1985.
[8] LIU Y, JR SHEPARD S W. Dynamic force identification based on enhanced least squares and total least-squares schemes in the frequency domain[J]. Journal of Sound and Vibration, 2005, 282(1/2): 37-60.
[9] 缑百勇, 陆秋海, 王波, 等. 利用固有频率异常值分析法检测螺栓拧紧力[J]. 振动与冲击, 2015, 34(23): 77-82.(GOU Baiyong, LU Qiuhai, WANG Bo, et al. Bolt tightening force detection using outlier analysis of structural natural frequencies[J]. Journal of Vibration and Shock, 2015, 34(23): 77-82.(in Chinese))
[10] 白会彦, 杜建镔. 某减速机构刚度分析及测试[J]. 计算机辅助工程, 2016, 25(5): 7-11.(BAI Huiyan, DU Jianbin. Stiffness analysis and test on decelerating mechanism[J]. Computer Aided Engineering, 2016, 25(5): 7-11.(in Chinese))
[11] GOLUB G H, VAN LOAN C F. Matrix Computation[M]. 4th ed. Baltimore: Johns Hopkins University Press, 2013.
[12] KARLSSON S E S. Identification of external structure loads from measured harmonic responses[J]. Journal of Sound and Vibration, 1996, 196(1): 59-74.
[13] LMS Egineering Innovation. Inverse force identification: tutorial[Z]. 2012.