基于不完备静态数据的结构损伤两阶段识别方法

杨万锋,杨 骁

(上海大学土木工程系,上海200444)

各类结构在服役期因受各种因素的影响而不可避免地产生性能的退化,因此对结构进行定期健康检测和损伤识别尤为重要.学者们针对不同的工程结构提出了许多新的损伤识别方法[1-6].这些方法的共性之处在于需要得到结构模型的完整测量数据.但在实际工程中,由于传感器数量或者测量条件等的限制,实际测量的结构自由度小于结构分析模型的自由度.因此,如何在不完备的实测数据以及考虑随机测量噪声引起的实测数据误差条件下,得到满足一定精度要求的损伤识别结果是一个值得关注和研究的课题.

不完备测试数据下的结构损伤识别问题已有部分研究成果.Sheena等[7]不考虑噪声影响,基于不完备测量数据,利用样条插值函数得到了未测量数据,据此提出了一种利用静态测量数据修正结构模型的方法.Snayaei等[8]基于不完备测量数据,以外部荷载与内力之差最小为目标函数,提出了一种确定结构损伤的优化迭代策略.崔飞等[9]研究了基于静态应变及位移测量的结构损伤识别算法,并结合Gauss-Newton法和梯度法提高了算法的鲁棒性.Kim等[10]在精确有限元模型的基础上,提出了一种不完备振型下基于结构固有频率变化的损伤识别方法.李蕊等[11]建立了有限测量数据条件下,将有限元模型缩聚理论与自适应序贯非线性最小二乘法相结合的结构损伤识别方法.赵建华等[12]在振型扩阶的基础上,利用Moore-Penrose广义逆法,提出了一种不完备测量的结构损伤识别方法.Seyed等[13]分别利用不完备的模态和静态位移数据建立了相应的目标函数,利用模拟退火算法优化目标函数,从而确定了损伤位置和损伤程度.聂振华等[14]将结构动态响应以相空间的形式展开,根据损伤前后相空间拓扑结构的变化,提出了新的损伤因子.李世龙等[15]基于模态扩充思想,利用损伤结构实测模态扩充数据与完整结构计算模态数据,建立了结构损伤控制方程,并提出运用信赖域优化算法进行结构各单元的刚度损伤识别.Chen等[16]研究了传感器的最佳位置确定,并基于多损伤定位准则,给出了基于模态频率变化的结构损伤程度识别方法.Sun等[17]利用不完备模态数据提出了基于自适应随机游走马尔可夫链蒙特卡罗方法的模型参数算法.

与结构动力学方程相比,结构静力平衡方程只与结构的刚度有关,与阻尼、质量等动力特性无关.因此,基于结构静力响应的结构损伤识别可有效避免因质量,特别是阻尼等的不确定性对损伤识别精度的影响.由于目前各种静态测量设备和技术已先进成熟,因而用较低的成本即可得到结构相当准确的变形或应变值.考虑到实际工程中实测数据有限而待识别参数较多,本工作结合基于静态应变能的损伤定位指标和基于有限元缩聚法损伤程度计算,提出了一种两阶段损伤识别方法.本方法利用静态应变能变化的损伤定位指标,对结构可能的损伤进行定位识别,在确定可能损伤位置的基础上,利用缩聚法建立损伤状态方程,并采用模拟退火算法求解损伤状态方程,从而得到结构的准确损伤位置和损伤程度.为验证本方法的有效性和可靠性,对一个5单元超静定,一个13单元静定和一个10单元超静定平面桁架结构的损伤识别进行了数值模拟.结果表明,在测量位移不完备的情况下,本工作所提出的两阶段结构损伤识别方法可对结构损伤位置和损伤程度进行准确识别,且对测量噪声具有较强的鲁棒性.

1 损伤定位

从理论上讲,静力识别中待识别的结构损伤参数可以是结构单元刚度矩阵中的任何变量.但一般而言,结构材料弹性模量E的变化最能直接反映损伤对结构及其构件受力所产生的影响.Banks等[18]和Luo等[19]给出了各向同性弹性材料结构的损伤:

式中,Eu和Ed分别为结构损伤前后材料的弹性模量,d(x,y)为结构材料局部损伤状态的损伤分布函数,其取值范围在0∼1之间,其中d(x,y)=0表示结构材料未损伤,d(x,y)=1表示结构材料完全损伤.

通常,结构静态平衡有限元方程可表示为

式中,[K]为结构的总刚度阵,{u}为结构总节点位移矢量,{F}为结构等效节点载荷.

式中,i和k分别表示荷载工况和单元,M和N分别表示荷载工况数和单元数,[Kku]为单元k的未损伤单元刚度矩阵,为未损伤时单元k在荷载i作用下的节点位移矢量,上标T表示转置.

显然,若单元k未损伤,则损伤前后该单元的刚度矩阵不变;若单元k出现损伤,则其损伤后的刚度矩阵为其中dk为单元k的平均损伤程度,是一个需确定的未知量.为简化计算,用未损伤时的单元刚度矩阵来代替损伤后的单元刚度矩阵这样,在静态荷载i作用下,损伤结构每个单元的静态应变能可近似为

荷载i作用下单元k的应变能增量∆Eik为

由此可定义荷载i作用下单元k相对于未损伤结构总应变能Eui的相对应变能增量∆Erik为

为减小随机噪声对节点位移的影响,定义M种荷载工况下的平均应变能增量为

与文献[20]类似,定义单元k的损伤定位指标为

由式(8)定义的损伤定位指标Dk具有如下性质:①若待识别结构是静定的,根据相对应变能增量∆的定义,则对于损伤单元k,Dk>0,对于未损伤单元k,Dk=0;②若待识别结构是超静定的,由于超静定结构内部单元损伤时会引起内力重分布,因此Dk=0表明该单元未损伤,而Dk>0包含真实的损伤单元以及若干未损伤单元.也就是说,Dk>0只表明可能的损伤单元,但该单元是否损伤,需利用单元k的局部损伤状态参量dk进行判断,即dk>0表示单元损伤,dk=0表示单元未损伤.

同时需要指出的是,在上述计算方法中,未损伤结构的节点位移由未损伤结构的有限元模型获得,而损伤结构未测量自由度上的位移由静态凝缩法获得[21-22].

2 损伤程度

假定由损伤定位指标Dk>0确定L个可能的损伤单元为k=k1,k2,···,kL,相应的损伤状态参量为dkr,r=1,2,···,L,此时,损伤结构的静力平衡方程为

式中,[Kd]=[Kd(dk1,dk2,···,dkL)]为损伤结构的总刚度矩阵,{ud}为损伤结构的节点位移矢量.

考虑到实际测量时,由于传感器数量有限或受测量条件的限制,不能测量结构所有的自由度位移.此时,将结构的节点位移矢量{ud}拆分为测量自由度上的位移{}和未测量自由度上的位移{},而相应的等效节点载荷矢量{F}拆分成相应的{Fm}和{Fs},其中下标m和s分别表示测量自由度和未测量自由度.

基于此,式(9)可改写为

由此可得

式(11)为确定状态参量dkr(r=1,2,···,L)的非线性损伤状态方程.考虑测量噪声的影响,将式(11)的求解转化为如下的最小值问题:在约束0 6 dk16 1,0 6 dk26 1,···,0 6 dkL6 1条件下,寻求dkr(r=1,2,···,L),使得

式中,d=(dk1,dk2,···,dkL),k·k表示矢量的模.

数值求解最小值问题式(13)的方法很多,但为了避免最小值解的搜索过程陷入局部最小解,可采用模拟退火法[23-24]进行求解,从而得到损伤状态参量dkr(r=1,2,···,L).对于单元kr0,若dkr0=0,则单元kr0未损伤;若dkr0>0,则单元kr0损伤,且损伤程度为dkr0.

3 数值算例

为验证本工作所建立的结构损伤两阶段识别方法的有效性和可靠性,首先对文献[20]中的5单元超静定平面桁架进行损伤定位,并与文献[20]的结果进行对比.5单元超静定平面桁架的有限元模型单元(编号①∼⑤)和节点(编号1∼4)如图1所示.设结构材料的弹性模量Eu=210 GPa,横截面积A=5×10−5m2,且假定结构的单元②、单元③、单元④分别存在20%,25%和10%不同程度的损伤,损伤前后该桁架的位移由有限元模拟给出.

图1 5单元超静定平面桁架结构Fig.1 5 element statically indeterminate planar truss structure

图2给出了文献[20]和本工作提出的损伤定位方法对该桁架在上述损伤工况下的损伤定位结果.可见,文献[20]的损伤定位方法出现了误判,即把单元④判定为未损伤单元,与实际的损伤情况不符;而本方法能初步判断单元②、单元③、单元④、单元⑤为可能损伤单元.该判断结果包含了所有的损伤单元.可见,本方法能更全面的给出可能损伤单元,为后续更准确的损伤判断奠定了基础.

图2 5单元超静定平面桁架损伤定位Fig.2 Damage localization of the 5 element statically indeterminate planar truss

为考察测量噪声对损伤识别结果的影响,假定测量位移由损伤结构有限元模拟位移叠加测量噪声构成,即若{}为损伤结构的有限元模拟位移,则测量位移为

式中,η为噪声程度,Nm为测量位移的自由度数,rand()为在区间[0,1]的均布随机数.

选用2个桁架结构进行损伤位置识别和损伤程度确定:一为13单元静定平面桁架结构,二为10单元超静定平面桁架结构,且2个算例均考虑了2种损伤工况,以期能较全面了解损伤单元数目及位置对识别结果的影响.

3.1 13单元静定平面桁架结构

13单元静定平面桁架结构的有限元模型单元(编号①∼13)和节点(编号1∼8)如图3所示.设结构材料的弹性模量Eu=210 GPa,横截面积A=5×10−5m2,且假定结构存在如表1所示的2种损伤工况Ⅰ和Ⅱ.为了减少随机噪声的影响,考虑了如表2所示的2种荷载情况,且2种损伤工况Ⅰ和Ⅱ中节点7的竖向位移和节点8的竖向及水平位移均未测量.

图3 13单元静定平面桁架结构Fig.3 13 element statically determinate planar truss structure

表1 13单元静定平面桁架损伤工况Table 1 Damage cases of the 13 element statically determinate planar truss

表2 13单元静定平面桁架荷载情况Table 2 Load cases of the 13 element statically determinate planar truss

3.1.1 损伤定位

图4和5给出了2种损伤工况Ⅰ和Ⅱ下13单元静定平面桁架结构单元的损伤定位结果.由图4可见,在无噪声的情况下,只有单元⑥的损伤定位指标D6>0,其他单元的Dk=0(k 6=6),与假定损伤工况Ⅰ的损伤单元一致.同样,从图5可以看出,在无噪声的情况下,只有单元③、单元⑤和单元⑥的损伤定位指标Dk>0(k=3,5,6),其他单元的Dk=0(k 6=3,5,6),亦与假定损伤工况Ⅱ的损伤单元一致.对于分别具有3%和5%均布随机噪声情况的损伤定位,从图4和5中也可以分别看出,噪声对损伤定位指标Dk的数值有一定的影响,但影响程度有限,不会造成误判.当噪声水平在3%时,损伤工况Ⅰ和Ⅱ下单元的损伤定位指标Dk值与无噪声时的最大相对误差分别为2.41%和3.46%;当噪声水平在5%时,相应的Dk值的最大相对误差为4.60%和5.48%.同时,未损伤单元的Dk值较损伤单元小一个数量级,很容易被剔除.因此,即便存在测量噪声,由缩聚法[21]得到的未测量位移仍能对静定平面桁架的损伤位置进行较准确判断.

图4 13单元静定平面桁架损伤工况Ⅰ的损伤定位Fig.4 Damage localization for the damage caseⅠof the 13 element statically determinate planar truss

图5 13单元静定平面桁架损伤工况Ⅱ的损伤定位Fig.5 Damage localization for the damage caseⅡof the 13 element statically determinate planar truss

3.1.2 损伤程度

在确定损伤位置的基础上,图6和7给出了2种损伤工况Ⅰ和Ⅱ下13单元静定平面桁架结构单元的损伤程度识别结果.可见,在没有噪声干扰时可准确判断2种工况的损伤程度,而当存在3%和5%噪声水平干扰时,2种工况下的损伤识别结果存在一定的误差:损伤工况Ⅰ下的损伤程度识别结果与真实损伤程度之间的最大相对误差分别为3.2%和4.8%,损伤工况Ⅱ下的最大相对误差分别为4.0%和4.4%.数值结果同时表明:噪声水平越高,损伤程度识别误差越大,但对噪声具有较强的鲁棒性.

图6 13单元静定平面桁架损伤工况Ⅰ的损伤程度Fig.6 Damage degree for the damage caseⅠof the 13 element statically determinate planar truss

图7 13单元静定平面桁架损伤工况Ⅱ的损伤程度Fig.7 Damage degree for the damage caseⅡof the 13 element statically determinate planar truss

3.2 10单元超静定平面桁架结构

10单元超静定平面桁架结构的有限元模型单元(编号①∼⑩)和节点(编号1∼6)如图8所示.设结构材料的弹性模量Eu=206.8 GPa,横截面积A=3.226×10−4m2,且假定结构存在如表3所示的2种损伤工况Ⅰ和Ⅱ,并考虑了如表4所示的2种荷载情况,且2种损伤工况中节点4的竖向及水平位移均未测量.

表3 10单元超静定平面桁架损伤工况Table 3 Damage cases of the 10 element statically indeterminate planar truss

表4 10单元超静定平面桁架荷载情况Table 4 Load cases of the 10 element statically indeterminate planar truss

图8 10单元超静定平面桁架结构(mm)Fig.8 10 element statically indeterminate planar truss structure(mm)

3.2.1 损伤定位

图9和10给出了2种工况下10单元超静平面桁架结构单元的损伤定位结果.对于损伤工况Ⅰ,由图9可见,结构的可能损伤单元为单元④、单元⑤和单元⑩,而损伤工况Ⅰ的实际损伤单元只是单元⑤和单元⑩,因此单元④为误判单元.造成误判的原因是当超静定结构中某个或某几个单元损伤时,结构会出现内力重分布,导致某些未损伤单元的静态应变能增加,使其损伤定位指标Dk>0,从而造成误判.对于损伤工况Ⅱ,由图10可见,单元④、单元⑤、单元⑥为损伤单元,与假设损伤工况Ⅱ的损伤单元一致.此时,虽然损伤结构也存在内力重分布,但未引起单元损伤的误判.另外,从图9和10中还可以发现,噪声不会引起损伤单元的误判.

3.2.2 损伤程度

对于10单元超静定平面桁架结构的损伤工况Ⅰ,由损伤定位可知,单元④、单元⑤和单元⑩为可能的损伤单元.在此基础上,采用模拟退火法计算单元④、单元⑤和单元⑩的损伤状态参量dk(k=4,5,10),结果如图11所示.可见:单元④的损伤状态参量d4几乎为0,因此该单元未损伤;单元⑤和单元⑩为损伤单元,且当噪声水平不超过5%时,识别的单元损伤程度与真实损伤程度的相对误差未超过5%.对于损伤工况Ⅱ,图12给出了损伤程度识别结果.可以看出,在无噪声干扰的情况下,损伤程度能被准确识别,但当有3%和5%测量噪声干扰时,损伤程度判断有一定的误差,然而识别的单元损伤程度与真实损伤程度的相对误差小于5%.

图9 10单元超静定平面桁架损伤工况Ⅰ的损伤定位Fig.9 Damage localization for the damage caseⅠof the 10 element statically indeterminate planar truss

图10 10单元超静定平面桁架损伤工况Ⅱ的损伤定位Fig.10 Damage localization for the damage caseⅡof the 10 element statically indeterminate planar truss

图11 10单元超静定平面桁架损伤工况Ⅰ的损伤程度Fig.11 Damage degree for the damage caseⅠof the 10 element statically indeterminate planar truss

若不先进行损伤定位的初步判断,而直接利用文献[13]的方法对10单元超静定平面桁架结构所有单元的损伤状态参量进行求解,则会出现计算收敛速度慢、误差大的情况,甚至出现误判.对于工况Ⅰ在3%噪声干扰的情况下,单元⑤的损伤程度的相对误差为9.5%,且得出单元①有5.5%的损伤程度;在5%噪声干扰的情况下,单元⑤和单元⑩的损伤程度判断误差分别为11.5%和13.5%,远大于本方法的相对误差.对于工况Ⅱ在3%噪声干扰的情况下,单元④、单元⑤和单元⑥的损伤程度判断误差分别为5.8%,6.5%和16.3%,并且会出现单元⑧和单元⑩的误判,显示单元⑧的损伤程度为13.3%,单元⑩的损伤程度为7.9%;在5%噪声干扰的情况下,实际损伤单元的损伤程度判断误差分别为10.8%,7.5%和22.8%,同时还误判单元⑧和单元⑩的损伤程度分别为20.2%和10.6%.可见,本工作提出的结构损伤两阶段识别方法具有较明显的优势.

图12 10单元超静定平面桁架损伤工况Ⅱ的损伤程度Fig.12 Damage degree for the damage caseⅡof the 10 element statically indeterminate planar truss

4 结束语

本工作提出了一种不完备静态变形数据条件下的结构损伤两阶段识别方法,即利用损伤定位指标确定结构的可能损伤位置,建立损伤状态方程,利用模拟退火法得到准确的损伤位置和损伤程度.本方法通过第一步可能损伤位置的确定,减少了损伤状态参量的计算,达到了减少计算量、提高识别精度的目的.具体数值算例结果表明:本方法能较准确地识别出结构的损伤位置和损伤程度,同时在噪声水平不大于5%时具有较强的鲁棒性.与文献[20]的损伤定位指标相比,本工作所提出的损伤定位指标能更全面地识别出结构可能的损伤位置,同时剔除了部分未损伤的单元.

本方法在大型复杂结构损伤识别中的具体应用仍有待进一步研究,包括如何确定最优的测量自由度数目及测量位置等.这些问题将在后续工作中进一步研究.