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

广义柔度矩阵可利用低阶模态参数近似表示[7],因此受到广泛关注[8,9,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,4,5,6]广泛应用于结构损伤识别的初步预判.模态应变能定义为

MSEji=12φTiKujφi, i=1,2,⋯,k, j=1,2,⋯,ne,   (4)

其中:MSEji表示第j个单元在第i个振型下的模态应变能;Kuj表示结构损伤前第j个单元扩阶后的单元刚度阵.将第j个单元在第i个振型下的模态应变能MSEji除以整个结构在第i个振型下的总模态应变能,可得

nMSEji=MSEji∑j=1neMSEji.   (5)

再令

nmMSEj=∑i=1knMSEjik,   (6)

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

MSEBIj=max{0,(nmMSEj)d−(nmMSEj)u(nmMSEj)u}.   (7)

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

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

MSECj=(nmMSEj)d−(nmMSEj)u.   (8)

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

MMSEBIj=max{0,MSECj−γmax(MSEC)},   (9)

其中γ=(∏j=1ne∣∣MSECj−μ∣∣)1/ne,μ=(∑j=1neMSECj)/ne,γ表示MSECj与其平均值μ之差绝对值的几何平均值.本文设当MMSEBIj>0时,则预判该单元为损伤单元.

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

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

Kuφui=λuiMφui,  i=1,2,⋯,n,   (10)

式中:Ku和M分别为损伤前n阶刚度与质量矩阵;λui=ω2ui和φui分别为第i个固有频率平方与关于质量矩阵M的归一化振型.

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

Kd=Ku−∑j=1neαjKuj,   (11)

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

Kdφdi=λdiMφdi,  i=1,2,⋯,n,   (12)

式中λdi=ω2di和φdi分别为损伤后结构的第i个固有频率平方与关于质量矩阵M的归一化振型.

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

Fgd(α)=Fd(MFd)l=ΦdΛ−1dΦTd(MΦdΛ−1dΦTd)l=ΦdΛ−1−ldΦTd,   (13)

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

Fgmd=∑i=1k1λ2diφdiφTdi,   (14)

式中振型φdi满足φTdiMφdi=1(i=1,2,…,k).使用前k个模态数据可构建损伤结构的广义模态柔度矩阵.

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

Fgmd≈∑i=1k1λ2uiφuiφTui+∑j=1neαj∂∂αj(∑i=1k1λ2diφdiφTdi)|α=0,   (15)

其偏导数部分可展开为

Aj=∑i=1k[−2∂λdi∂αj1λ3uiφuiφTui+1λ2ui∂φdi∂αjφTui+1λ2uiφui(∂φdi∂αj)T]∣∣α=0,   (16)

其中,在α=0处频率与振型对损伤参数αj的导数∂λdi∂αj和∂φdi∂αj可按Nelson方法[11]计算.将方程(12)对损伤参数αj求导,因假设损伤与质量无关,所以∂M∂αj=0,再令α=0可得

(∂Kd∂αj∣∣α=0−∂λdi∂αj∣∣α=0M)φui+(Ku−λuiM)∂φdi∂αj∣∣α=0=0,   (17)

将方程(17)两边左乘φTui,并结合式(10),(11),可得

∂λdi∂αj∣∣α=0=−φTuiKujφui.   (18)

方程变为

(Ku−λuiM)∂φdi∂αj∣∣α=0=(Kuj−φTuiKujφuiM)φui.   (19)

首先确定φui分量中绝对值最大者,设其下标为k0.将(Ku-λuiM)中第k0行和第k0列中的元素全部设为0,但将其第k0行和第k0列交叉位置元素设为1,得到的矩阵记为Bi.将向量(Kuj-φTuiKujφuiM)φui中的第k0行元素设为0,记为bi.解方程组Biatj=bi可得方程(19)的一个特解atj,其通解可写成如下形式:

∂φdi∂αj∣∣α=0=atj+cijφui.   (20)

将方程φTdiMφdi=1两边关于参数aj求导并令α=0,得

φTuiM∂φdi∂αj∣∣α=0=0.   (21)

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

cij=−φTuiMatj.   (22)

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

⎛⎝⎜∂λdi∂αj∂φdi∂αj⎞⎠⎟∣∣∣∣α=0=(−φTuiKujφuiatj−φTuiMatjφui).   (23)

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

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

∑i=1k1λ2uiφuiφTui+∑j=1neαjAj=∑i=1k1λ2eiφeiφTei,   (24)

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

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

Aα=b,   (25)

其中b=vec(∑i=1k1λ2eiφeiφTei)−vec(∑i=1k1λ2uiφuiφTui).这里:将Aj按列拉直成一个n2维列向量,而A是由Aj(j=1,2,…,ne)依次拉直并从左到右排列后形成的一个n2×ne矩阵;α是由各单元的损伤参数组成的一个ne维列向量;vec(∑i=1k1λ2uiφuiφTui)和vec(∑i=1k1λ2eiφeiφTei)分别是矩阵∑i=1k1λ2uiφuiφTui和∑i=1k1λ2eiφeiφTei按列拉直形成的n2维列向量.

3、解法

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

min∥Aα−b∥22.   (26)

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

α=lsqr(A,b,10−10,200)   (27)

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

αj∶=αj1+αj (j=1,2,⋯,ne)   (28)

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

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

步骤1)初始假设r0=ne−−√,r1=∥α∥2.令η=|r0−r1|,当η小于某个值η0时,即认为迭代已达到收敛精度,即可退出程序循环,本文设η0=0.001.

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

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

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

图1损伤识别流程

4、数值实验

4.1模型实验

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

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

φzeij=φeij(1+ε1Rij(max|φui|)),   (29)λzei=λei(1+ε2Ri),   (30)

其中:φzeij和λzei分别为损伤后加入测量噪声的第i阶振型分量与频率平方;而φeij和λei分别为损伤后未加测量噪声的第i阶振型分量与频率平方;ε1和ε2为所加的噪声水平;j表示振型分量;Rij和Ri均是均值为0、方差为1的高斯白噪声产生的随机数;max|φui|是损伤前第i阶振型分量绝对值的最大值.将测量得到的且带有噪声的500组模态数据分别代入广义模态柔度矩阵∑i=1k1λ2eiφeiφTei,然后求其平均值,再代入式(26)求解.

图2框架模型

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

例1考虑如图2所示的框架结构.该框架分为21个单元,基本结构参数如下:每个单元长l=0.5m;弹性模量为E=210GPa;截面惯性矩为I=8.33×10-6m4;密度为ρ=7800kg/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种损伤定位方法定位结果的比较

图4损伤情形1下框架结构3种求解方法损伤识别结果的比较

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

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

图5损伤情形2下框架结构3种损伤定位方法定位结果的比较

图6损伤情形2下框架结构3种求解方法损伤识别结果的比较

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

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

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

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

图7桁架模型

图8损伤情形1下桁架结构3种损伤定位方法定位结果的比较

图9损伤情形1下桁架结构3种求解方法损伤识别结果的比较

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

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

图10损伤情形2下桁架结构3种损伤定位方法定位结果的比较

图11损伤情形2下桁架结构3种求解方法损伤识别结果的比较

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下框架结构两种方法损伤识别结果的比较

图13损伤情形2下框架结构两种方法损伤识别结果的比较

图14损伤情形1下桁架结构两种方法损伤识别结果的比较

图15损伤情形2下桁架结构两种方法损伤识别结果的比较

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

参考文献：

[8]巩文龙,常军,刘大山,等.基于量子粒子群优化算法的结构损伤识别[J].动力学与控制学报,2015,13(5):388-393.

[10]赵博,徐自力,阚选恩.拉杆转子界面局部脱开识别的改进广义柔度矩阵法[J].振动工程学报,2017,30(5):724-729.

[13]李晶.基于广义柔度矩阵的结构损伤识别研究[D].长春:吉林大学,2011.

[15]杨秋伟,王学航,李翠红.基于高次广义柔度灵敏度的结构损伤识别[J].固体力学学报,2019,40(2):157-168.

谢少鹏,吴柏生,钟慧湘.基于广义模态柔度矩阵的结构损伤识别[J].吉林大学学报(理学版),2020,58(03):518-526.

基金:国家自然科学基金(批准号:11672118)