一种针对局部结构的有限元模型修正方法

2020-12-28范新亮王彤夏遵平

航空学报 2020年12期

范新亮，王彤，夏遵平

南京航空航天大学机械结构力学及控制国家重点实验室，南京 210016

有限元模型修正首先在航空航天领域提出，发展至今日已形成了一个庞大的理论体系，并且广泛应用于运载火箭、卫星、航天飞机、飞机、直升机等结构的响应与载荷预示、颤振分析、振动控制[1]。模型修正可基于不同的特征进行，主流的有基于频域模态参数、频响函数(Frequency Response Function,FRF)的方法以及基于时域动响应的方法。其中时域方法虽抗噪性相比前两者稍弱，但获取实验数据简便，且对于包含非线性的结构[2]、运行状态的结构[3]等难以利用模态或频响函数数据进行模型修正的情形仍适用。Modak等[4]通过数值仿真详细对比了基于模态与基于频响函数的方法[5-9]的差异及收敛性。基于模态参数的方法简单有效，在多数情形下能给出较满意的修正结果，但也存在一些缺陷：如模态参数提取过程引入了误差及不确定性[10]，待修正参数数目受到测试模态数量的限制。基于频响函数的方法则避免了提取模态参数带来的误差，且具有大量数据可供修正以解决待修正参数众多的问题。

普遍的模型修正算法均是对整体结构有限元模型进行修正，但有限元模型不确定性常存在于某些关键的局部区域，例如飞机各个部件的连接区域、机械结构的关节点、转子系统的支承等，且这类局部区域所在的子结构(局部结构)无法从整体结构中拆卸，或其组装至整体结构上的动态特性与独立状态有所差异，此时往往需要利用整体结构的测试信息对其有限元模型进行迭代修正，计算效率较低。对此，Weissenburger[11]提出了一种预测局部结构的参数变化对整体结构振动特性影响的方法。Zhu等[12]使用主模态缩聚残余结构，并以局部区域对应的子结构有限元模型中的物理参数和残余结构的模态参数作为优化变量，极小化结构动力学特征方程的残差来修正整体结构的有限元模型。翁顺等[13]提出了一种基于子结构的有限元模型修正方法，当结构仅局部参数变化时只需分析局部区域所对应子结构即可求解整体结构的模态特征量灵敏度矩阵。上述研究都有效提高了大型结构有限元模型修正的效率，但均为基于模态参数的方法。在基于频响函数的方法中，Guo等[14]推导了一种考虑残余结构约束的局部结构应变频响函数的近似计算方法，并根据其灵敏度建立了针对局部区域进行模型修正的方程，该方法与Zhu等[12]所提方法均对残余结构进行了缩聚而保留了完备的局部结构，然而并未对缩聚误差对修正结果的影响进行讨论。综上，若能得到整体结构动响应与独立状态下的局部结构的关系，则可直接利用整体结构测试数据对独立状态下的局部结构模型进行修正，从而显著提高效率。

本文根据频响函数解耦理论得到的整体结构与独立状态的局部结构两者间频响函数的关系式，将包含待修正参数的局部结构动刚度矩阵与残余结构有限元频响函数耦合得到整体结构拟合频响函数，通过极小化其与测量值的残差得到待修正参数的估计值，从而推导了以局部结构动力学矩阵表示的整体结构频响函数残差关于待修正参数的灵敏度方程。因此，当结构建模误差仅发生在局部区域时，利用该方程可对分离出的局部结构单独进行有限元模型修正，将残余结构与更新的局部结构重新装配即得到修正后的整体结构有限元模型。最后，通过仿真和实验算例验证了该方法的有效性和抗噪性。

1 理论背景

1.1 局部结构频响识别方法

为了利用整体结构测试频响函数对独立状态下的局部结构有限元模型进行修正，首先需要建立整体结构与局部结构之间频响函数的关系，即频响函数解耦问题[15-19]。

图1 局部结构与残余结构Fig.1 Local structures and residual structures

整体结构、残余结构及局部结构的输入输出关系分别为

HAFA=XA,HRFR=XR,HLFL=XL

(1)

式中：HA、HR和HL分别为整体结构、残余结构和局部结构的频响函数矩阵。

根据界面力协调条件及位移协调条件有

(2)

设在残余结构VR内部节点中所选测试自由度对应的位移矢量为XI*=PI*TXI，其中PI*为XI*T各列在XIT中的位置矩阵。定义伪测试自由度为VR内部测试自由度与完备界面自由度的并集，则V与VR在伪测试自由度上的位移矢量XA*、XR*为

(3)

(4)

相应载荷矢量为FA*=PA*TFA及FR*=PA*TFR。而VL在完备界面自由度上的位移矢量XL*为

(5)

其载荷矢量为FL*=PLTFL。在上述位移矢量XA*、XR*及XL*上，式(1)成为

(6)

式中：频响函数子矩阵HA*、HR*与HL*分别为

并将式(6)按VR内部测试自由度与完备界面自由度写为分块形式

(7)

(8)

(9)

(10)

其中

(11)

式中：HR*与HC*分别为由残余结构与局部结构有限元模型计算得到的频响函数。

式(10)即为有限元模型整体结构V与残余结构VR、局部结构VL间频响函数的关系式，是后文推导针对局部结构有限元模型修正公式的基础。由式(3)注意到伪测试自由度包含完备的界面自由度XJ，因此计入所有界面自由度时解耦式(10) 是精确的，若忽略某些自由度则界面力协调条件不完全成立，将引起近似误差。

1.2 局部结构有限元模型修正方法

假设整体结构有限元模型中局部区域存在建模误差，其待修正参数为θ，而残余结构不存在或近似不存在误差。由式(10)知局部结构与残余结构频响函数耦合所得整体结构拟合频响函数为

(12)

取由参数θ确定的残差函数为

(13)

式中：ej为激励自由度j在伪测试自由度中的位置向量；HA*ej为整体结构伪测试自由度对应的测试频响函数。

采用极大似然估计来识别参数θ，使残差vj取极小值，即参数θ所决定的整体结构频响函数拟合值与测量值误差最小。极大似然估计需要测量数据的“先验”噪声信息，当满足互不相关的白噪声假设时，“先验”信息可由测量数据的标准方差来代替[20]。对于具有No个输出、Ni个输入的频响函数，共包含NoNi个随机变量。假设所有频率点处的频响函数测量值相互独立，则此NoNi个随机变量的联合概率密度函数为

(14)

式中：H为共轭转置符号；ωk为第k个频率点(为行文简洁，推导中非必要处略去ωk)；v由Ni个残差向量vj组合而得；C为频响函数噪声的协方差矩阵，当各列噪声相互独立时其为Cj组成的对角矩阵，且

(15)

根据极大似然估计原理得到等价极大似然函数为

(16)

式中：Nf为频率点数目。将式(16)简写为

L(θ)=Q(θ)HQ(θ)

(17)

(18)

根据Newton-Gauss方法知参数θ的迭代估计式为

J(θ(r))JH(θ(r))d(r+1)=-J(θ(r))Q(θ(r))

(19)

式中：d(r+1)为第r+1个迭代步的参数增量，即

d(r+1)=θ(r+1)-θ(r)

(20)

(21)

式(12)、式(18)代入式(21)左端雅可比矩阵可化简得

(22)

其中等效频响函数定义为

(23)

式中：I为单位矩阵；DL(θ)为局部结构的动刚度矩阵，其参数化的表达式[12]为

(24)

AL(ω)=

(25)

(26)

式(24)代入式(22)得

(27)

其中

(28)

将式(27)代入式(21)得与ωk及ej对应的迭代方程

(29)

(30)

式(30)即为求解待修正参数的迭代方程，简写为

S(θ(r))d(r+1)=f(θ(r))

(31)

式中:f(θ)为待修正参数θ所确定的整体结构频响函数拟合值与测量值的加权残差；S(θ)为以局部结构的动力学矩阵MLe等表示的残差f(θ)关于参数θ的灵敏度。综上，对整体结构测试得到HA*ej后，再由残余结构计算得到HR*，即可由式(31)对独立状态下的局部结构有限元模型进行修正。

1.3 算法实现相关问题

1.3.1 频率点的筛选准则

实际测试中，反共振区频响函数易受到噪声污染，而共振区不仅噪声小，对参数变化的灵敏度高，并且其幅值能有效反映阻尼参数。因此在进行模型修正的频率范围内选择共振区的频率点，更容易得到稳健的参数识别结果[21]。所提方法以频响函数测试值与拟合值的残差为最小化目标函数，因此选取残差较小的频率点参与参数估计可有效减小测试噪声的扰动影响，故此处以频响函数幅值相关性系数[14]作为筛选频率点的准则

(32)

1.3.2 局部结构模型修正流程

局部结构模型修正流程如图2所示：

1) 将含建模误差的局部结构VL从整体结构V中分离出来。确定完备界面自由度及其对应的位置矩阵PR,J与PL、整体结构伪测试自由度及其对应位置矩阵PA*以及激励自由度及其对应位置向量ej。需要注意的是伪测试自由度并非真实进行测试的自由度，因为完备的界面自由度信息几乎不可能由试验得到。并且伪测试自由度为真实测试自由度与完备界面自由度的并集。

图2 局部结构模型修正流程Fig.2 Process of model updating of local structures

6) 重复步骤3)～5)，直至满足停止准则。迭代结束后将局部结构的修正参数代回整体结构，即得到修正后的整体有限元模型。

注意到当残余结构有限元模型规模较大时，HR*的计算量较大，因此HR*仅计算一次且在迭代中不重复该计算。由此可见，所提方法将无需进行模型修正的残余结构的计算工作置于迭代之外，效率相比传统方法得到了提高。进一步地，可对残余结构有限元模型进行缩聚后计算HR*，但这将引入额外误差。此外，若直接从整体结构测试数据中提取局部结构信息[23]，即

(33)

(34)

而后利用所得HL*对独立状态下的局部结构有限元模型进行修正，将存在以下缺陷：需得到完整的测试频响函数HA*，在实际中工程量巨大甚至难以实现；且提取HL*的过程需求解大量逆矩阵，使噪声的扰动放大而导致抗噪性较差[24]。

2 算例分析

2.1 仿真算例

采用某喷气式飞机模型验证本文方法的有效性及抗噪性。模型整体结构V及局部结构VL、残余结构VR如图3所示，采用Shell63板单元建模。其频带内主要弹性模态为机翼一阶、二阶、三阶弯曲及尾翼的弯曲。设置建模误差存在于机翼所在局部区域，而机身与尾翼所在残余结构VR的模型参数无误差。其中VL又划分为蒙皮、襟翼、内骨架等5个组(图3)，每个组的弹性模量、密度及阻尼系数作为待修正参数θ。

以某组参数η*对应的整体结构有限元模型作为真实结构，在所选激励自由度及伪测试自由度上计算其频响函数后添加10%噪声(白噪声与有色噪声各5%)得到整体结构测试频响函数HA*ej；VR参数与η*中相应分量ζ*(η*=[ζ*θ*])取相同，计算得伪测试自由度上残余结构频响函数HR*；VL的参数θ相对η*中相应分量θ*设置一定的偏差作为待修正参数初始值，计算得初始局部结构的频响函数HL(θ(0))，并与HR*重构得整体结构频响函数初始拟合值，与测量值对比如图4所示。

图3 喷气式飞机局部结构与残余结构Fig.3 Local and residual structures of a jet

图4 修正前整体结构频响函数对比Fig.4 Comparison of global structure FRF before updating

图5 修正后整体结构频响函数对比Fig.5 Comparison of global structure FRF after updating

图6 修正前后频响函数幅值相关性对比Fig.6 FRF magnitude correlation coefficient plot before and after updating

图7为θ(r)相对真实参数比值的迭代历程以及修正前后对比，E1至E5、ρ1至ρ5分别为5个分组的弹性模量及密度。由图知修正后大多参数收敛于真实值附近，仅ρ3及ρ5识别精度较差，推测其原因为目标残差函数对其灵敏度较小，使得该参数易受噪声干扰。也正因此，在所分析频段内，该参数对整体结构的动态特性作用较小，其误差不会影响修正后的整体模型预测动力学响应的准确性。通过该仿真算例表明所提方法能在噪声干扰下有效地对含建模误差的局部区域进行修正，使得整体结构模型的动力学特性与真实情形一致。

图7 迭代历程及修正前后参数值对比Fig.7 Iterative process and comparison of parameters before and after updating

2.2 实验

图8 三角机翼飞机实验系统Fig.8 Experimental system for delta-winged aircraft

为进一步验证所提方法对实际结构局部区域进行模型修正的有效性，对三角机翼飞机进行了实验，如图8所示。厚3 mm的机翼与尾翼通过螺栓与厚6 mm的机身连接，采用Shell63板单元建模并划分为8个组(图9)，且螺栓连接简化为固接。由于该结构在低频范围内主要为机翼的振动，且机翼与机身连接处存在不确定的建模误差，因此将机翼及其连接部分分离出来作为待修正的

图9 三角机翼飞机有限元模型Fig.9 FE model of delta-winged aircraft

局部结构VL，剩余部分则为近似无误差的残余结构VR。实验中测试频率范围设置为100 Hz，取800条谱线。整体结构采用弹性绳悬挂方式模拟自由—自由边界条件，传感器位于两侧机翼上，在均匀分布的41个测试点上进行敲击，获得其频响函数HA*ej。为抑制噪声[23]，须合理选择包含较多结构信息的测试自由度集合。