基于子结构解耦的连接特性识别新方法1)

范新亮 王 彤 夏遵平

(南京航空航天大学振动工程研究所,机械结构力学及控制国家重点实验室,南京 210016)

引言

机械结构中广泛存在着各种连接形式,如螺栓连接、铆接、焊接以及胶接等.连接结构的动力学特性决定着装配体的动响应预测精度,因此准确地对其进行辨识显得尤为重要.而由于连接结构的多样性和复杂性[1-3],难以建立其一般的解析模型[4],因此常基于子结构解耦技术并利用模态或频响函数等测试数据对连接特性进行辨识.该技术从装配体和残余结构的动力学行为中提取连接结构的动态特性[5],其中基于模态试验数据的方法中振型的测试误差极易导致识别结果不理想[6],因此更适于有限元环境下应用;而基于频响函数试验数据的方法则由于测试简单、数据信息充足且精度较高,得到广泛关注与研究[7-12].

Okubo 和Miyazaki[13]最早将基于频响函数的子结构解耦技术应用于连接特性识别,经过不断地发展与改善,许多方法被提出,并在一些实际结构中得到了应用[14-17].虽然形式不尽相同,但几乎所有方法的辨识公式都与Ren 和Beards[18]提出的基本公式相似或等价[19],其思路均为利用实测的装配体频响函数与实测或有限元模型计算的残余结构频响函数来获取连接特性.连接特性识别所面临的主要难题之一是装配体结构界面自由度频响函数的测量[19].对此,Čelič等[4-5]在Ren 等所提基本公式上进行了改进,详细讨论了忽略部分界面自由度及内部自由度对辨识结果的影响,并对比了逐频求逆的直接法与对各个频率点进行平均的最小二乘法,然而其在忽略部分界面自由度时的辨识结果不佳.此外,采用旋转自由度的测试技术[16]及未测频响函数估计技术[19-20]是解决该问题的可行方案,但却增加了测试或计算的成本,且额外地引入了测试误差或估计误差.连接特性识别的另一难题是辨识方程对测试噪声的高度敏感性,即测试数据微小的误差也可能导致结果与真实特性完全无关,且不同形式的求解方程在辨识精度上有着较大差异[18],这是由过多的矩阵求逆以及系统方程本身的严重病态造成的.对含噪数据进行降噪是减小噪声干扰最直接的方法,Peeters 等[7]采用模态模型对测试频响函数进行滤波并应用于子结构解耦,尽管测试噪声得以滤除,但滤波后的频响函数对真实值的拟合程度仍难以控制.替代降噪技术的有Zhang 等[21-22]提出的一种噪声自校准参数辨识方法,以谐波基对噪声校准量参数化并与系统参数通过误差函数共同进行估计,有效提高了参数辨识的准确性.在优化辨识公式的形式方面,Wang 等[23]减少了对噪声有放大作用的矩阵求逆运算,使得方法对实测数据的适用性增强,但由于公式中的简化使得在高频时辨识效果不佳.Wang等[19,24]根据基本公式推导了统一形式的辨识方程,并利用该方程将实测与估算的未测频响函数一起引入辨识过程抑制噪声,提高了识别结果的准确性,然而由于采用了逐频求逆的直接法,在由动刚度矩阵序列估计连接结构物理参数的过程中将产生偏差.

为此,本文在子结构解耦基本方程基础上,提取其测试自由度分量并推导了显含连接动刚度矩阵的辨识方程,该方程直接利用可测的装配体频响函数进行辨识,而无需针对其界面自由度信息进行未测频响函数的估计或旋转自由度及部分不可测界面自由度的测量.而后将该方程写为增量格式,得到了连接特性识别的迭代形式,其残差相对全量格式显著减小,因而数值稳定性更好[25],且通过迭代可使得由连接动刚度与残余结构频响函数重构的装配体频响函数拟合值不断逼近于测量值.最后,将动刚度曲线以多项式进行拟合,得到了拟合系数的估计方程,并筛选合适的频率点联立方程组求解,减少了病态矩阵的求逆运算且改善了抗噪性.弹簧质量系统数值算例和双梁连接结构试验表明该方法辨识效果较好,能有效地对实际连接特性进行识别.

1 理论背景

1.1 连接结构动力学特性识别方程

连接结构、残余结构及其组成的装配体结构如图1 所示,令I 代表残余结构或装配体结构部分内部节点的自由度,J 代表残余结构或连接结构全部界面节点的自由度,A 和R 分别代表装配体结构和残余结构,则有输入-输出关系

图1 连接结构,残余结构与装配体结构Fig.1 Joint,residual and assembly structure

其中XI,XJ与FI,FJ分别为装配体结构在内部自由度上、界面自由度上的位移矢量与的力矢量,而X与F为由其组合而成的位移矢量与力矢量,顶标-与^则分别代表残余子结构与连接结构的相应物理量.HA,HR和HJ分别为装配体结构、残余结构和连接结构的频响函数矩阵.

由界面自由度上位移与力的协调关系有

可得基于子结构解耦识别连接特性的基本公式[19]

其中PR,J为筛选HR中界面自由度所在列的布尔矩阵,即PR,JTHR=[HRJIHRJJ]

连接结构动力学特性识别问题一般忽略连接结构的内部自由度(需满足连接件内部自由度上无外力的假设),且由XJ所指为完备界面自由度知

设实际对装配体结构进行振动测试的自由度上位移矢量为X*,且X*=P*TX,其中P*亦为布尔矩阵.取P*T左乘式(5)得

式(8)经移项后左乘广义逆(P*THRPR,J)+得

注意此处求取广义逆必须满足行大于列的条件,即要求测试自由度的数目应大于界面自由度数目.上式两端左乘HC并将式(6)和式(7)代入得

移项并在等式两端左乘DJ得

而由AX=B的极小二乘解为X=(A)+B知上式等价于

将式(9)记为

式(10)与式(8)形式一致,其中左、右等效频响为

式(10)即为显含连接动刚度DJ的辨识方程,是后文推导增量格式以及引入多项式拟合的基础.其中HR通过校准的残余子结构有限元模型计算得到,因此对式(10)两端右乘装配体结构振动测试的激励自由度布尔矩阵PI后,辨识过程所需测试数据仅为装配体在所选测试自由度上的频响函数P*THAPI,而无需测量或估计其界面自由度的量.且式(10)由式(5)提取测试自由度上的平衡方程后经矩阵变换得到,界面协调条件仍成立,因此仍是精确的.此外,所选测试自由度能准确辨识DJ的必要条件为拥有足够多有效信息.

1.2 识别方程增量形式

连接结构界面自由度数目较多时,由于待求未知量的数目增加而使得辨识精度明显下降,且噪声对辨识结果的扰动显著增强.为提高方程数值稳定性并通过迭代逼近最佳辨识值,将式(10)变换为增量形式.若给定连接结构动刚度初值DJ0,欲识别其真实值DJ,由式(10)分别有

又由式(11)知右等效频响函数的差量为

将其代入式(14)可得

其中I为单位矩阵,将上式记为

因此,由式(5),式(16)和式(17)可得识别连接特性的迭代方程

其中各迭代量为

式(18)和式(19)右乘PI后,给定有限元模型计算数据HR、装配体测试数据P*THAPI与初值后,即可由迭代更新连接结构动刚度.

全量形式的辨识公式(10)以残余结构与装配体频响函数的残差为极小化目标函数,增量式(18)则以范数远小于前者的装配体频响函数测试值与初始拟合值的残差为目标函数,参数识别的数值稳定性得到提高[25],且减小了噪声对辨识结果的干扰.

1.3 动刚度曲线多项式拟合

若采用直接法[5]逐个频率点依次求解式(18),对含噪声矩阵(且在界面自由度较多时常因测试自由度信息不足而呈现病态)求逆后将使误差放大.且连接结构对装配体结构动力学特性的影响在部分频率点几乎可以忽略,使相应频率点的辨识结果对噪声极为敏感[16].为此,以多项式函数拟合动刚度曲线,代入式(18)得到关于拟合系数的方程,并选择合适的频率点联立求解,则原问题转化为频域内极小化残差函数的参数估计问题.

连接结构动刚度矩阵第m行第n列(m,n=1,2,···,NJ,NJ为界面自由度数目)对应元素DJmn(ωk)可用多项式函数表示为(1+iωkβ)α0mn-ωk2α2mn,其中β为刚度比例阻尼系数,ωk为第k个频率点,而拟合系数α0mn,α2mn实际上为连接结构刚度矩阵、质量矩阵的对应分量,则DJ的增量可写为

取ej为振动测试时筛选位移矢量X中第j个激励自由度所在行的位置向量,在式(18)两端右乘ej得

将式(20)代入(21)得

对上式两端进行矩阵的列拉直得

式中V(·) 表示矩阵的列拉直即为第g个迭代步的残差向量,而为残差关于拟合系数的灵敏度矩阵,γ(g)则为拟合系数增量.

式(22)对应于频率点ωk以及第j个激励自由度,对k(k=1,2,···,Nf,Nf为频率点数目)及j(j=1,2,···,NI,NI为激励自由度数目)组合可得频域上多输入多输出的参数估计方程组f(g)=S(g)γ(g).当所联立的频率点及激励自由度包含足够多系统有效信息时,可有效减小待估计参数较多(即连接结构界面自由度较多)时灵敏度矩阵的病态程度,且通过联合全频域上所有测点与激励点的信息可提高辨识的准确性.

方法迭代收敛性的证明[21]过程为,若所选测试自由度集、激励自由度集、频率点集使得形成的S(g)满秩,记且

αT为真实拟合系数时,相应残差f(αT)为0,即αT是迭代方程(23)的不动点.而

且由S(g)满秩即S(g)TS(g)可逆及有界性知存在常数M1,M2使得

故||α(g+1)-αT||≤M1M2||α(g)-αT||2=d(g)||α(g)-αT||,当任意α(g)满足||α(g)-αT||＜1/M1M2时,有d(g)＜1,迭代方程(23)为压缩映射,α(g)将收敛至不动点αT.由此表明给定了合适的连接特性初值后,方法将收敛于真实值.

1.4 算法实现

1.4.1 方程归一化及频率点的选择

式(22)中系数矩阵Sj的各列之间存在显著的幅值差异,导致其条件数较大,因此按式(24)(25)对其各列进行归一化

其中[],r与[]r,分别代表矩阵第r列与第r行,norm(.)代表对列取范数.

实际测试试验中,共振区受噪声污染小,且对参数变化灵敏度更高[26],故信号信噪比较低时可仅选择共振区的频率点参与辨识过程,当信噪比较高时则可选择给定带宽的共振区邻域.此外,由于式(22)以装配体频响函数测量值与拟合值的残差范数为极小化目标,因此可将其幅值相关性系数[27]作为筛选频率点的准则: 相关性较好的频率点,其残差较小且可信度较高;反之,相关性较低则说明该频率点可能受噪声污染较为严重,且残差较大不利于准确辨识.幅值相关性系数αa计算式为

其中,ℜ 代表取实部,σmin为所允许的相关性最小值,实践中发现取σmin＞0 时识别效果较好.

1.4.2 连接结构动力学特性识别流程

连接特性辨识流程如图2 所示,具体过程如下.

图2 连接结构动力学特性辨识流程Fig.2 Workflow of the joint dynamic properties identification

步骤1: 确定完备界面自由度及其对应的布尔矩阵PR,J,装配体结构的测试自由度及其对应的布尔矩阵P*,激励自由度及其对应的位置向量ej(j=1,2,···,NI).

步骤2: 根据激励自由度ej及测试自由度P*对装配体结构进行振动测试,得到测试频响函数P*THAej.选择共振区(或其邻域)内的频率点序列L1并利用校准后的残余子结构有限元模型计算相应频响函数HRP*,HRPR,J.置g=0.

根据式(26)在 L1中选择参与辨识过程的频率点 L2.

步骤4: 由式(11)、式(19)计算所选频率点序列 L2对应的左、右等效频响函数HEl*,HEr*ej及转换矩阵TE*(g) .

步骤5: 计算频率点序列 L2中ωk处的灵敏度矩阵及残差向量代入式(22),对k和j组合形成超定方程组f(g)=S(g)γ(g),并进行归一化.求解拟合系数增量γ(g)后得到ΔDJ(g),更新连接动刚度矩阵DJ(g +1)=DJ(g)+ΔDJ(g),置g=g+1.

步骤6: 重复步骤3～5,直至满足停止准则,即||f(g)||≤ε,输出连接结构动刚度矩阵的拟合系数.

当残余结构有限元模型规模较大时,其频响函数HRP*,HRPR,J的计算量较大且是决定算法效率的关键一环,因此需在迭代前计算频率点序列 L1对应的HR.且该模型须经过精细的校准,否则HR将包含较大的误差,导致错误的识别结果.

2 算例分析

2.1 仿真算例

为验证所提方法的正确性及抗噪性,以具有10 个平动自由度的简单线性多自由度系统(图3)进行数值仿真.将其视为由连接结构与两个残余子结构组成的装配体(残余子结构1、2 分别在自由度处与连接结构耦合),且该连接结构同时包含质量、刚度及阻尼效应.系统各参数为1.0×107N/m,残余结构与连接结构的刚度比例阻尼系数β均取1.0×10-7.

图3 连接结构与残余结构耦合Fig.3 Coupling of joint and residual structure

装配体结构计算的频响函数添加噪声后作为测试值HAej,残余结构计算的频响函数作为所需HR.连接动刚度矩阵所取初值DJ0为

在无噪声情形下,由本文方法利用测试频响函数P*THAej(P*不包含界面自由度按1.4 中所述流程进行辨识,并与利用完备频响函数HAej由逐频求逆的直接法[5]进行对比.由图4(a)知两种方法辨识的DJ与HR重构的装配体频响函数均与真实值完全一致,验证了所提方法仅利用所选测试自由度(无需包含完备的界面自由度)上的装配体频响函数识别连接特性的正确性.

为检验所提方法的抗噪性,在前述算例的频响函数真实值中加入10%的白色噪声与10%的有色噪声后再次进行识别,由图4(b)知两种方法在测试噪声较小的共振区识别效果均较好,而在其余噪声较大的频率点,所提方法的准确度高于直接方法.对比两者所识别DJ的对角元素与非对角元素可知(图5),本文方法辨识结果受噪声的影响远小于直接法.

图4 装配体结构频响函数重构值与真实值对比Fig.4 Comparison of reconstructed and accurate assembly FRF

图5 所辨识连接结构动刚度实部的对比Fig.5 Comparison of real part of the identified joint dynamic stiffness

2.2 实验验证

2.2.1 胶接

为验证所提方法识别实际连接结构的有效性,首先对一个T 形结构进行了实验.该装配体结构由两根相同的长为500 mm、宽为50 mm 的梁(残余结构)在与一铁块在界面处胶接而成,且忽略连接界面的非线性行为.实验中采用了N-Modal 模态测试系统,PCB 力锤以及PCB 三轴传感器.测试频率范围设置为1000 Hz,取1600 条谱线.如图6 所示,装配体结构采用弹性绳悬挂的方式模拟自由-自由边界条件,传感器位于梁I 中心位置,在均匀分布的126 个测试点上进行锤击,获得其频响函数P*THAPI.辨识公式中(P*THRPR,J)+要求测试自由度数目大于连接结构自由度数目,然而为抑制噪声[24]并避免P*THRPR,J出现零列或呈现病态,应布置更多包含结构有效信息的测试自由度.

图6 含一处连接的耦合结构及其测试装置Fig.6 The coupled structure with one joint and its test setup

采用Shell63 板单元对残余结构建模,并按材料参数划分为两个组.利用自由状态的测试数据对单根梁进行模型修正[27-31]后,有限元模型计算的频响函数与其测试值匹配程度较高(图7),由此得到校准的残余结构有限元模型.其材料参数如表1 所示,其中D与E代表密度与弹性模量,且残余结构与连接结构的刚度比例阻尼系数均取为1.0 × 10-7.由该校准模型计算得残余结构频响HRP*和HRPR,J.

表1 梁I 和II 的模型材料参数Table 1 The material parameters of Beam I and Beam II

图7 校准的残余子结构有限元模型计算频响与测试值对比Fig.7 Comparison of the residual structure's FRF: calibrated FEM and experimental model

残余结构与连接结构的耦合如图8 所示.连接结构两端的虚拟节点分别通过MPC 约束依赖于梁I、梁II 被胶接区域内的节点,使得虚拟节点的平移及转动自由度上的位移代表了胶接区域内节点的位移均值,且虚拟节点上的力被均匀分配至胶接区域内的节点上[1].初值DJ0的取法为: 以密度为90 000 kg/m3、弹性模量为3.0×107Pa 的梁单元连接两虚拟节点,其形成的动刚度矩阵作为与HR所重构的装配体频响函数初始拟合值与测试值对比如图9 所示.

图8 T 形梁的连接结构与残余结构耦合Fig.8 Coupling of Joint and residual structure of T-Beam

图9 迭代前装配体频响函数重构值与测试值对比Fig.9 Comparison of reconstructed and experimental assembly FRF before iteration

按1.4 中所述流程识别得到的DJ拟合系数常数项、二次项作为连接结构的刚度、质量矩阵,与残余结构有限元模型重组可得到修正的装配体有限元模型,其频响函数结果与测试值匹配较好(图10),表明本文方法能有效识别实际连接结构动态特性.值得注意的是,由于实际连接结构刚度、质量矩阵的对称性、稀疏性及其他特性将对DJ的拟合系数产生一系列约束,需在辨识方程中予以考虑.迭代过程中共振区序列 L1中各频率点的装配体重构频响与测试频响的幅值相关性如图11 所示,可见其随着迭代步的增加而不断提高,即残差范数||f(g)||逐渐减小,验证了算法的收敛性.

图10 迭代后装配体频响函数重构值与测试值对比Fig.10 Comparison of reconstructed and experimental assembly FRF after iteration

图11 频响函数相关性的迭代历程Fig.11 Iterative process of FRF correlation

2.2.2 螺栓连接

为进一步验证所提方法对较复杂连接动力学特性的识别能力,对一个包含两处螺栓连接的双梁结构进行了实验.该装配体结构由两根相同的长度为500 mm 的L 形梁(残余结构)通过螺栓连接而成,测试中使用的数据采集系统如图12 所示.测试频率范围设置为2500 Hz,取1600 条谱线.装配体结构以海绵垫支撑以模拟自由-自由边界条件,三轴传感器位于螺栓附近,在均布的138 个测试点上进行锤击,获得其频响函数P*THAPI.采用Shell63 板单元对单根L 形梁建模,并按材料参数划分为3 个组.经模型修正后得到校准的残余结构有限元模型,其计算的频响函数与测试值吻合较好(图13).材料参数如表2 所示,且残余结构与连接结构的刚度比例阻尼系数均取为1.0 × 10-7.利用该模型可较准确地计算出包含完备界面自由度信息的残余结构频响函数.

表2 梁I 和II 的模型材料参数Table 2 The material parameters of Beam I and Beam II

图13 残余结构校准后模型频响与测试值对比Fig.13 Comparison of the residual structure's FRF: calibrated FEM and experimental model

残余结构有限元模型与连接结构的耦合如图14所示,与实验一类似,通过MPC 约束将梁I、梁II 的螺栓压紧区域内的节点与连接结构耦合,其中压紧区域取螺栓直径的两倍范围.该结构包含两处相同的连接,取约束方程使得两者连接参数相同,即由于两处连接无共同自由度,因此总体连接动刚度矩阵为初值的取法为:以密度为7900 kg/m3、弹性模量为1.8×1011Pa的梁单元连接虚拟节点,其动刚度矩阵作为与HR重构的装配体频响函数初始拟合值与测量值对比如图15 所示.

图14 L 形梁的连接结构与残余结构耦合Fig.14 Coupling of Joint and residual structure of L-Beam

图15 迭代前装配体频响函数重构值与测试值对比Fig.15 Comparison of reconstructed and experimental assembly FRF before iteration

同样地,所辨识连接结构与残余结构重组得到装配体结构有限元模型,其频响函数结果与测量值对比如图16 所示,尽管在550～ 650 Hz 的高频段有所差异,但150～ 550 Hz 的频率范围内均匹配较好.由此证明了所提方法能有效地对较复杂的多处连接特性进行识别.图17 为动刚度曲线相对初值归一化的拟合系数随迭代的变化,可以发现部分参数振荡较为明显,装配体结构动响应对这些参数通常有低敏感性的特点因而易受噪声干扰且难以准确辨识,也正因该特点,其辨识误差对对连接结构预测装配体动力学特性的准确性影响较小[16].

图16 迭代后装配体频响函数重构值与测试值对比Fig.16 Comparison of reconstructed and experimental assembly FRF after iteration

图17 拟合系数的迭代历程Fig.17 Iterative process of fitting coefficient

3 结论

本文提出一种新的连接特性识别方法,直接利用可测的装配体频响函数进行辨识,避免了旋转自由度及部分不可测界面自由度的测量或估计.通过具有收敛性质的增量方程增加了界面自由度数目较多时辨识的数值稳定性;并以多项式拟合动刚度曲线,在所选频率点上联立方程得到参数估计公式,减少了对病态矩阵的求逆,提高了辨识结果的准确性.通过数值仿真算例了验证所提方法的正确性及抗噪性,通过含胶接、螺栓连接等结构的试验验证了方法识别较复杂的实际连接结构特性的有效性.