基于模型缩聚-频响函数型模型修正的子结构损伤识别方法

方有亮, 娄佳琪, 张 颖*, 李宗娆, 侯童非

(1.河北大学建筑工程学院， 保定 071002； 2.河北省土木工程监测与评估技术创新中心，保定 071002)

随着社会的发展进步，实际工程中结构的规模越来越大，结构形式越来越复杂，对大型结构的损伤识别方法的研究是必不可少的[1]。大型复杂结构在实际试验中面临着测试信息少、修正参数过多、计算量过大等问题。运用传统损伤识别的方法对整个结构进行求解，会造成计算难度增加，迭代计算消耗大量时间，而且还可能造成难以收敛等问题。为了解决上述困难，可以将整体结构分成多个子结构进行损伤识别。由于子结构是相互独立的，因此运用子结构方法可单独识别各子结构。侯吉林等[2]对约束子结构模型修正方法的提出与研究，使得利用少量信息就可准确识别子结构损伤情况；翁顺等[3]利用子结构在主模态的特征信息进行模型修正；Li等[4]提出了一种基于频域动态响应重构的子结构损伤识别方法；周林仁等[5]提出了针对大跨斜拉桥结构的子结构模型修正方法；雷鹰等[6]对大型结构在部分观测下的子结构损伤诊断法进行了研究。基于此，采用子结构的Guyan静力缩聚方法，该方法无须考虑子结构模态坐标的概念，不进行子结构主模态分析，方法简单，操作性更强。

对于频响函数型模型修正方法而言，采用直接实测得到的频响函数而避免了模态法进行模态参数识别的误差，应用更加广泛。Hu等[7]提出了不依赖于灵敏度分析的正交正则模态修正方法的思想；徐张明等[8]对一种灵敏度分析的频响函数型模型修正方法进行研究；刘金梅等[9]利用基于灵敏度分析的模型修正方法对井架结构进行分析识别；Wang等[10]将Kriging模型引入算法，提出了一种新的基于加速度频响函数的模型修正的方法。对于频响函数型的模型修正方法，要采用实测频响函数。而实际实验中，许多自由度信息是无法得到的，如梁的转角自由度，这就会出现有限元分析模型中的自由度数目远大于实测的自由度数目。因此，将理论模型进行自由度缩聚，从而满足理论自由度与实际的相匹配。张以帅[11]提出将缩聚模型代入到模型修正的方法中进行迭代修正。该方法能够准确识别参数，同时避免了灵敏度分析带来的误差和计算效率。但经过多次计算发现，对于测点数目少的识别，使用上述方法，在迭代计算过程中会出现奇异。

为此，将子结构方法与模型缩聚-频响函数型模型修正迭代方法结合，运用到损伤识别中。首先，将结构划分出多个子结构，采用子结构Guyan静力缩聚方法，缩聚掉不识别子结构的内部自由度，与待识别的子结构组装形成组装子结构;然后,再次进行自由度缩聚，仅保留测点自由度信息;最后，运用模型缩聚-频响函数性模型修正方法识别子结构。

1 子结构的Guyan静力缩聚

将整体结构划分成N个子结构，按内部和边界自由度对子结构的刚度矩阵和质量矩阵进行分区，则第i(1,2,…,N)个子结构的刚度矩阵、质量矩阵如式(1)所示：

(1)

式(1)中：上标I表示第i个子结构的内部自由度；上标B表示第i个子结构对接界面自由度。

在Guyan缩聚中选择对接界面自由度为主自由度(保留)，内部自由度为副坐标(缩减)，则子结构动力方程为

(2)

忽略副坐标上质量对整体的惯性力，由式(2)可得:

(3)

(4)

式(4)中:Ti为Guyan缩聚的转换矩阵,其公式为

(5)

(6)

这时，子结构在经过静力缩聚变换后的动力方程为

(7)

将子结构的缩聚方程组集在一起，利用对接界面位移协调和力的平衡条件就可以形成系统缩聚后的动力方程：

(8)

式(8)中：MG、KG分别为子结构缩聚组集后的系统质量矩阵和刚度矩阵；uB为对接界面位移向量。

2 模型缩聚-频响函数型模型修正的子结构损伤识别方法

模型缩聚-频响函数型模型修正迭代方法为

(9)

由于在测点较少的情况下，使用模型缩聚-频响函数型模型修正迭代方法，在迭代计算中较容易出现奇异，所以将子结构与该方法结合，进行损伤识别。结构的节点编号由阿拉伯数字“1,2，…”表示，单元号由“(1)，(2)，…”表示，测点编号为“①，②，…”。如图1所示，为一个三层平面刚架，该结构共有15个节点，18个单元，节点1、11、12处为固定端。以该三层平面刚架为例，推导模型缩聚-频响函数型模型修正的子结构损伤识别方法。

图1 整体结构节点和单元图Fig.1 The nodes and elements of global structure

假设在实验中测点为节点5、6、7、8、9、10的水平方向，则只能测得保留这6点水平自由度信息的实验频率ωt和实验频响函数Ht(ωt)，已知待识别子结构为刚度损伤，进行损伤定位和损伤定量。

(1)首先将整体结构划分成2个子结构，分别对两个子结构的节点和单元重新编号，如图2所示，子结构1中共有11个节点，12个单元，子结构2中有7个节点，6个单元。其次运用子结构Guyan静力缩聚，对子结构2进行简化，将内部自由度缩聚到边界自由度上。子结构2的刚度矩阵、质量矩阵分别用K2、M2表示，Guyan缩聚的转换矩阵为Ta2，则缩聚后子结构的刚度矩阵、质量矩阵为

(10)

(2)将子结构2的内部自由度缩聚到附加边界自由度上，即缩聚到图2(b)中节点5、6、7上。再将缩聚后的子结构2的刚度矩阵和质量矩阵附加到子结构1刚度矩阵和质量矩阵的相应位置，形成组装子结构，如图3所示。组装子结构的质量矩阵、刚度矩阵分别用KK、MM表示。

图3 组装子结构Fig.3 Assembley substructure

(11)

设损伤的子结构1刚度矩阵为Kt1，则对于实验模型质量矩阵不变，刚度矩阵为

(12)

(3)组装子结构的有限元模型再次简化，进行Guyan缩聚，仅保留图3中节点5、6、7、8、9、10水平向自由度信息，得到最终的有限元简化模型并进行了测点编号，如图4所示。

图4 简化模型及测点编号Fig.4 Simplified model and the number of measuring points

缩聚的转换矩阵为Ta1，缩聚后的质量矩阵、刚度矩阵为

(13)

实验模型缩聚的转换矩阵为Tt1，实验模型缩聚后的刚度矩阵为

(14)

(4)由式(10)～式(14),可得：

(15)

(16)

式中：m为子结构1单元数(m=1,2,…,12)；Δαj为子结构1各单元刚度矩阵修正量；Kj为子结构1各单元刚度矩阵。

(5)将式(15)、式(16)代入式(9)，忽略阻尼得：

(17)

(6)判断修正量是否小于一个小值ε。如果不小于，则将各单元的修正量代入式(18)、式(19)，对组装子结构的刚度矩阵及子结构1中的各单元刚度矩阵进行修正。

(18)

(19)

子结构损伤识别流程图，如图5所示。

图5 子结构损伤识别流程图Fig.5 The flow chart of substructure damage identification

3 数值模拟

图6为一平面刚架。刚架总高0.84 m，总宽 0.56 m，每根梁的横截面积为0.031 m×0.008 m，梁的弹性模量E=2.1×1011Pa，质量密度ρ=7 850 kg/m3。

图6 三层平面刚架Fig.6 Three-layer plane frame

该刚架结构共有18个单元、15个节点、45个自由度。将该刚架结构拆分成两个子结构，如图2所示。子结构1共有12个单元、11个节点。子结构2共有6个单元、7个节点、9个内部自由度、9个界面自由度。子结构2进行子结构Guyan静力缩聚，得到组装子结构，如图3所示。最终进行自由度缩聚得到简化模型，如图4所示，保留节点5、6、7、8、9、10水平自由度信息。对比这三个模型的前6阶固有频率，如表1所示。

表1 各阶段固有频率对比Table 1 Comparison of natural frequencies in different stage

由表1可知，经过缩聚后，前三阶频率基本一致，误差不到0.2%，越高阶频率则误差越大，因此，取前三阶频响函数进行计算。频响函数中加入1%噪声干扰。

3.1 单一单元损伤

假设子结构1中2单元刚度损伤31%(弹性模量减少模拟刚度损伤)。

如图7所示，在测点①处激励，测点③得到的单损伤前后的频响函数曲线。曲线中的峰值的频率即为表2中对应的频率值。由图7、表2可以看出，单一单元损伤前后频响函数发生了变化。将每次迭代求出的子结构1的单元刚度矩阵修正量求和，得到子结构1中每个单元的损伤程度为

图7 单一单元损伤前后频响函数对比Fig.7 Comparison of FRFs before and after single element damage

表2 单一单元损伤前后频率对比Table 2 Frequency comparison before and after single element damage

Δαj=[-0.012 2 0.296 7 -0.018 1 -0.028 7

0.064 7 -0.056 2 0.019 1 -0.012 3

0.038 4 -0.023 8 -0.004 0 -0.025 0]T

(20)

单一单元损伤结果如图8所示，能够准确识别出损伤的位置及损伤程度，最终识别结果为子结构1中2单元刚度损伤29.67%，识别结果较准确。其他未损伤单元的识别结果最大误差在5%左右。

图8 单一单元损伤识别结果Fig.8 Single element damage recognition result

3.2 多单元损伤

3单元损伤19%，5单元损伤25%，11单元损伤10%。

如图9所示，在测点①处激励，测点③得到的多损伤前后的频响函数曲线。曲线中的峰值的频率即为表3中对应的频率。可以看出，多单元损伤前后频响函数发生了变化，前三阶频率分别降低了0.28、1.35、1.73。将每次迭代求出的子结构1的单元刚度矩阵修正量求和，得到子结构1每个单元的损伤程度：

图9 多单元损伤前后频响函数对比Fig.9 Comparison of FRFs before and after multi-element damage

表3 多单元损伤前后频率对比Table 3 Frequency comparison before and after multi-element damage

Δαj=[-0.032 8 0.040 2 0.152 9 -0.026 3

0.254 8 -0.005 6 -0.017 6 0.045 6

-0.043 5 -0.030 3 0.134 6 0.01]T

(21)

多单元损伤结果如图10所示，能够准确识别出损伤的位置及损伤程度，最终识别结果为子结构1中3单元损伤15.29%，5单元损伤25.48%，11单元损伤13.46%。能够较准确识别出单元损伤程度。其他未损伤单元的识别结果最大误差在4%左右。

图10 多单元损伤识别结果Fig.10 Multi-element damage recognition result

算例中设定的收敛精度为ε=1×10-5，单一单元损伤识别经过7次迭代后收敛至该精度。如图11所示，损伤单元从第4次迭代后参数识别进入收敛，直至7次迭代全部单元收敛到设定精度。损伤单元最终识别的刚度损伤为每次迭代结果之和(图11)。多单元损伤识别经过9次迭代，识别参数收敛至设定结果，如图12所示。从表4中的数据可以明确每次迭代的参数修正量，最终得出单元刚度的损伤程度。

图11 单一单元损伤迭代收敛情况Fig.11 Single element damage iteration convergence

图12 多单元损伤迭代收敛情况Fig.12 Multi-elements damage iterative convergence

表4 多单元损伤损伤单元收敛过程Table 4 The convergence of damaged elements in multi-element damage

4 结论

提出了基于模型缩聚-频响函数型模型修正的子结构损伤识别方法。并通过数值模拟算例，对该方法进行了分析和验证，得出以下结论。

(1)子结构的Guyan静力缩聚方法无须考虑子结构模态坐标的概念，不进行子结构主模态分析，方法简单，操作方便。

(2)通过对三层刚架数值模拟，在含噪声干扰下，对本文方法进行了验证，结果表明该方法简便，易于计算，结果准确。

(3)由于划分了子结构，因此可以识别任意子结构，相比于识别整体结构，更加高效，方便。