基于广义模态柔度矩阵的结构损伤识别

谢少鹏, 吴柏生, 钟慧湘

(广东工业大学 机电工程学院, 广州 510006)

结构损伤检测近年来受到广泛关注.Pandey等[1]利用结构柔度矩阵的变化检测结构损伤; Jaishi等[2]基于模态柔度残量的有限元模型修正检测损伤.由于实际工程中不可能测得所有阶模态参数, 因此用模态截断的柔度矩阵识别方法有较大误差.损伤定位是结构损伤识别的一个重要组成部分, 文献[3]基于模态应变能的定位方法利用损伤前后单元的模态应变能变化量作为损伤识别的指标, 文献[4-6]对该方法进行了深入研究.但上述方法使用的模态阶数较多, 定位的损伤单元数远超实际损伤单元数或者遗漏损伤单元.

广义柔度矩阵可利用低阶模态参数近似表示[7], 因此受到广泛关注[8-10].本文基于广义柔度矩阵引入广义模态柔度矩阵, 利用Nelson方法[11]计算其对损伤参数的敏感性.本文提出的损伤定位方法, 相比于现有方法, 在振型数据不完整及含有噪声的情况下, 该方法具有使用模态参数少, 定位的损伤单元数比实际损伤单元数略多且不遗漏损伤单元等特点.从损伤定位确定的损伤单元出发, 本文使用最小平方正交三角分解法(LSQR)[12]求解最小二乘问题, 逐步确定损伤位置与程度.该方法仅涉及要求导的模态, 可降低测量噪声的影响.

1 振型扩充及损伤单元预判

1.1 振型扩充

在结构有限元模型的每个节点上通常仅可放置适当数量的传感器, 此时可用某些振型扩充技术近似计算出那些未测量自由度的振型分量.本文采用如下公式完成振型扩充[13]:

(1)

其中：K和M分别表示n阶结构有限元模型的刚度和质量矩阵；λi为第i阶固有频率的平方;m和s分别表示已测量与未测量的自由度.由式(1)的第二式可得

Ksmφm+Kssφs=λi(Msmφm+Mssφs),

(2)

则未测量的振型部分可表示为

φs=-(Kss-λiMss)-1(Ksm-λiMsm)φm.

(3)

1.2 损伤单元预判

基于模态应变能(MSE)方法[3-6]广泛应用于结构损伤识别的初步预判.模态应变能定义为

(4)

(5)

再令

(6)

其中k表示预先确定所取的模态阶数.利用上述公式可分别计算出各单元损伤前后的参数(nmMSEj)u和(nmMSEj)d, 其中单元刚度矩阵均用损伤前的相应矩阵计算.Seyedpoor[5]利用完整模态数据给出了模态应变能指数MSEBIj:

(7)

文献[5]假设当MSEBIj大于某个值时, 则该单元刚度发生变化, 即该单元发生损伤.该方法虽然对损伤处敏感, 但在振型测量不完整及含有噪声的模态数据情形下, 存在较多误判损伤单元.

为了准确定位损伤位置并尽量减少误判单元数量, 同时不遗漏损伤单元, 本文提出一种新的判别损伤位置的方法.定义参数：

MSECj=(nmMSEj)d-(nmMSEj)u.

(8)

由于噪声等因素的影响, 使得某些未损伤单元的MSEBIj>0, 导致误判的单元数量较多, 故改进如下: 引入修正的模态应变能指数

(9)

2 确定损伤程度的控制方程

对损伤前的结构, 在不考虑阻尼情形下, 其频率与振型满足如下特征方程：

Kuφui=λuiMφui,i=1,2,…,n,

(10)

假设损伤前后质量矩阵M保持不变, 则损伤后刚度矩阵可表示为

(11)

其中:ne为单元个数;Kuj和αj(0≤αj≤1)分别为第j个单元损伤前的扩阶单元刚度矩阵与损伤程度, 损伤程度0表示未损伤, 而1表示完全损伤.损伤后结构的频率与振型满足如下方程：

Kdφdi=λdiMφdi,i=1,2,…,n,

(12)

令Fd和Φd分别表示损伤结构的柔度矩阵与模态矩阵, 则广义柔度矩阵[7]为

(13)

其中l=1,2,….当l=0时为通常的柔度矩阵, 当l=1时为一阶广义柔度矩阵[7].随着l增大, 高阶模态参数对广义柔度矩阵的影响迅速衰减.假设k个最低频率及其对应的振型可用.由方程(13)并设l=1, 引入广义模态柔度矩阵：

(14)

将式(14)中的广义模态柔度矩阵在结构损伤前, 即在α=0处做一阶Taylor展开并忽略高阶项, 有

(15)

其偏导数部分可展开为

(16)

(17)

(18)

方程变为

(19)

(20)

(21)

将式(20)代入式(21), 可得

(22)

于是, 特征值与特征向量导数为

(23)

利用式(16),(23)即可求得Aj(j=1,2,…,ne).

令有限元模型预测损伤结构的广义模态柔度矩阵与利用测量模态参数计算得到的广义模态柔度矩阵相等, 可得

(24)

式中φei和λei分别为损伤后测量所得的第i阶关于M归一化振型和频率的平方.

矩阵方程(24)可写成如下形式：

Aα=b,

(25)

3 解 法

不失一般性, 方程(25)可转化为如下最小二乘问题：

(26)

传统最小二乘法对右端项b存在噪声, 且当矩阵A条件数较大时会失效.Yan等[14]基于模态柔度矩阵, 使用非负最小二乘法识别结构损伤; 杨秋伟等[15]提出了反馈奇异值截断法识别结构损伤； Paige等[12]首次提出了一种基于Lanczos迭代法求解病态最小二乘问题的方法----最小平方QR分解(LSQR)法, 可提高计算精度.本文使用MATLAB软件包中LSQR函数：

α=lsqr(A,b,10-10,200)

(27)

求解方程(26).为提高解的精度, 可适当调高收敛精度和迭代次数.在得到损伤参数α后, 用加速公式[15]：

(28)

对数据进一步处理, 以得到更准确的结果.

当方程组维数较少时, 上述方法可获得满意的结果.而对于较复杂的结构, 在多个单元损伤及测试噪声影响下, 上述方法无法准确确定损伤程度.因此, 本文提出如下迭代策略：

步骤2) 令r0=r1, 选出α中大于0.03(该数值在识别程度较小的损伤单元过程中可适当下调)且小于1的分量, 并将其下标归入集合Θ.保留系数矩阵A中对应集合Θ的列向量, 去掉其余的列向量.

步骤3) 调用MATLAB中的LSQR函数, 用新系数矩阵重新计算得到新损伤程度值, 通过加速式(28)再次得到一组新损伤值.计算加速后损伤参数的2-范数r1, 返回步骤1).

通过上述步骤可过滤对求解结果无用甚至产生干扰的信息.在迭代过程中, 系数矩阵规模(列向量数)越来越小, 对于大规模最小二乘问题, 可提高计算效率.算法流程如图1所示.

图1 损伤识别流程Fig.1 Flow chart of damage identification

4 数值实验

4.1 模型实验

根据测量的自由度信息, 抽取广义柔度灵敏度矩阵Aj中与测量自由度对应的行和列元素, 并按列拉直成一个列向量, 依次组合相应灵敏度矩阵对应的列向量后, 形成结构灵敏度矩阵, 仍记为A.右端项损伤前后的广义柔度矩阵也分别进行类似处理, 得到的向量仍记为b, 再用本文方法进行求解.

为模拟实际工程中的测量误差, 分别给频率和振型施加1.5%和5%的高斯白噪声, 施加噪声公式[15]为

(29)

(30)

图2 框架模型Fig.2 Model of frame

本文使用框架和桁架两个算例, 比较损伤定位指数MSEBI[5]、损伤定位指数nMSEBI(normalized modal strain energy based index)[6]及本文方程(9)的损伤定位指数MMSEBI.先基于本文提出的定位方法确定损伤单元, 再分别利用传统最小二乘法、LSQR方法及迭代LSQR方法(本文方法)计算单元的损伤程度.每个算例分别设定两种损伤情形, 以验证本文方法的有效性.

例1考虑如图2所示的框架结构.该框架分为21个单元, 基本结构参数如下： 每个单元长l=0.5 m； 弹性模量为E=210 GPa； 截面惯性矩为I=8.33×10-6m4； 密度为ρ=7 800 kg/m3； 横截面面积为A=10-2m2； Poisson比为ν=0.3.在实际工程中可将损伤程度小于5%的单元视为未损伤单元.

例1仅利用结构前两阶模态.考虑各节点的轴向位移, 整个结构共有60个自由度.测量前两阶频率及对应振型的第10、第12和第14节点处的垂直分量与第2、第4、第6、第8、第15、第17、第19和第21节点处的轴向位移分量.

损伤情形1： 第8个单元损失15%的刚度.

此时, 3种损伤定位方法的定位结果如图3所示.由图3可见, 3种方法均可判断出损伤单元的位置, 基于MSEBI方法判断的损伤单元较多, 由MMSEBI方法判断的损伤单元数次之, 而利用nMSEBI方法定位的损伤单元较少.用3种方法计算的损伤程度如图4所示.由图4可见, 迭代LSQR算法的计算结果最接近实际损伤值.

图3 损伤情形1下框架结构3种损伤定位方法定位结果的比较Fig.3 Comparison of location results of three damage location methods for frame structure in damage case 1

图4 损伤情形1下框架结构3种求解方法损伤识别结果的比较Fig.4 Comparison of damage identification results of three solution methods for frame structure in damage case 1

损伤情形2： 第5、第11和第18个单元分别损失15%,20%和20%的刚度.

此时, 3种损伤定位方法定位结果如图5所示.基于nMSEBI方法在多损伤单元情形下, 出现了第18个损伤单元漏判.而基于MMSEBI方法仍然能识别出全部的损伤单元, 且与基于MSEBI方法相比定位损伤单元较少.基于上述损伤定位, 损伤程度计算结果如图6所示.由图6可见, 迭代LSQR方法仍然可较准确地确定损伤程度.

图5 损伤情形2下框架结构3种损伤定位方法定位结果的比较Fig.5 Comparison of location results of three damage location methods for frame structure in damage case 2

图6 损伤情形2下框架结构3种求解方法损伤识别结果的比较Fig.6 Comparison of damage identification results of three solution methods for frame structure in damage case 2

例2考虑如图7所示的桁架结构.该模型是一个由31个杆单元组成的桁架结构.假设直杆和斜杆单元长分别为l1=1 m和l2=1.414 m, 密度ρ=2 770 kg/m3, 横截面面积A=10-4m2, 弹性模量为E=70 GPa.

对该桁架结构, 测量前三阶固有频率及对应振型的第3、第5、第7、第8、第9和第11节点的水平与垂直振型分量.

损伤情形1： 第5个单元损失15%的刚度.

此时, 3种损伤定位方法所得结果如图8所示.由于所布置的传感器远离第5个单元, 故基于nMSEBI的方法漏判第5个损伤单元, 而本文方法虽然误判单元稍多, 但未出现漏判现象.而基于MSEBI的判断结果仍略多于前两种方法.3种方法损伤识别结果如图9所示.由图9可见, 本文方法求解结果准确, 而其他两种方法均出现了一些误判为损伤的单元.

图7 桁架模型Fig.7 Model of truss

图8 损伤情形1下桁架结构3种损伤 定位方法定位结果的比较Fig.8 Comparison of location results of three damage location methods for truss structure in damage case 1

图9 损伤情形1下桁架结构3种求解 方法损伤识别结果的比较Fig.9 Comparison of damage identification results of three solution methods for truss structure in damage case 1

损伤情形2： 第7、第15和第28个单元分别损失20%,25%和30%的刚度.

此时, 3种损伤定位方法所得结果如图10所示.由图10可见, nMSEBI方法定位丢失了第15个损伤单元, 而MMSEBI方法与MSEBI方法一样, 没有丢失损伤单元信息, 且误判单元数量比MSEBI方法少.图11比较了3种求解方法所得结果.由图11可见, 迭代LSQR法的计算结果最接近实际情形, 而另外两种方法则均出现了损伤单元误判.

图10 损伤情形2下桁架结构3种损伤 定位方法定位结果的比较Fig.10 Comparison of location results of three damage location methods for truss structure in damage case 2

图11 损伤情形2下桁架结构3种求解 方法损伤识别结果的比较Fig.11 Comparison of damage identification results of three solution methods for truss structure in damage case 2

4.2 与模态柔度矩阵法的比较

为进一步说明本文方法的有效性, 下面将本文方法与模态柔度矩阵法所得结果进行比较, 两种方法均为同一噪声水平下先完成损伤定位, 再用迭代LSQR法求解.

对于模态柔度矩阵法, 只需将式(16)改成如下形式：

(31)

将式(25)中的b改成如下形式：

(32)

首先, 考虑例1.

损伤情形1： 第2个单元损伤15%的刚度.

损伤情形2： 第7、第15和第20个单元分别损伤15%,20%和20%的刚度.

其次， 考虑算例2.

损伤情形1： 第4个单元损伤15%的刚度.

损伤情形2： 第5、第13和第24个单元分别损伤15%,20%和20%的刚度.

损伤识别结果如图12～图15所示.从识别精度上看, 模态柔度矩阵法在仅有低阶模态数据的情形下, 有时易产生损伤单元误判(图15), 而广义模态柔度矩阵法所得结果更接近实际损伤.

图12 损伤情形1下框架结构两种 方法损伤识别结果的比较Fig.12 Comparison of damage identification results of two methods for frame structure in damage case 1

图13 损伤情形2下框架结构两种 方法损伤识别结果的比较Fig.13 Comparison of damage identification results of two methods for frame structure in damage case 2

图14 损伤情形1下桁架结构两种 方法损伤识别结果的比较Fig.14 Comparison of damage identification results of two methods for truss structure in damage case 1

图15 损伤情形2下桁架结构两种 方法损伤识别结果的比较Fig.15 Comparison of damage identification results of two methods for truss structure in damage case 2

综上, 本文提出了一种新的损伤定位方法(MMSEBI), 减少了未知量的个数并避免损伤单元的漏判.采用部分低阶模态参数构建广义模态柔度矩阵, 利用Nelson方法推导了其敏感性, 进而将损伤识别问题转化为最小二乘问题.因求解该问题常用的最小二乘法失效, 故本文应用LSQR法逐步确定损伤位置和程度以提高识别精度.数值实验验证了本文方法适用于多种不同的结构损伤, 且与基于模态柔度矩阵的损伤识别法相比, 其计算结果精度更高.