考虑刚度和质量耦合效应的结构弹性成像方法1)

2023-10-29付君健彭铁川李帅虎周祥曼

力学学报 2023年9期

关键词：特征值弹性函数

付君健彭铁川李然*,2) 李帅虎周祥曼李响

* (水电机械设备设计与维护湖北省重点实验室,湖北宜昌 443002)† (三峡大学机械与动力学院,湖北宜昌 443002)

** (长江三峡通航管理局,湖北宜昌 443002)

††(武汉理工大学交通与物流工程学院,武汉 430063)

引言

弹性成像[1-3]是一种医学影像技术,可以把生物体组织的弹性模量信息转换为可视化图像,显示组织或器官病理变化的位置和大小.对于机械装备结构,损伤会导致结构特定部位弹性模量发生改变,继而使结构力学性能发生改变.如果将弹性模量与机械结构损伤表征关联,可以通过弹性成像方法获取弹性模量分布以实现结构的损伤识别.机械装备结构极易产生缺损、夹杂、结冰和积水等多种损伤,而大多数损伤本质上是结构特定部位的弹性模量发生偏移.因此,研究一种先进的弹性成像方法可实现结构多种损伤类型的全局识别,有助于实现重大机械工程装备的服役安全评价.

目前,结构损伤识别方法主要分为先进成像方法和模型修正方法.先进成像方法是一种无损检测技术,它通过计算机分析处理穿越物体的响应信号,并融合相应算法重构物体内部信息,从而以图像的形式直观显示物体的损伤位置.近年来,许多先进成像方法被用于结构损伤识别,如超声成像[4]、射线成像[5]等.先进成像方法可以对结构进行实时成像,具有损伤直观可视等优点.但大多数成像方法属于局部区域损伤识别技术,需要知道损伤大致区域作为先验信息,对机械结构的损伤识别效率较低.为解决局部损伤识别方法的缺陷,具有全局识别能力的模型修正方法应运而生.模型修正方法以有限元模型为基础,首先定义一组表征损伤位置和程度的修正参数,修正参数如单元刚度折减系数、截面面积和边界条件等,然后通过修正损伤结构的有限元模型将模型响应和实测响应之间差异最小化,最后找到修正参数的最优解得到结构的损伤位置和程度[6-8].模型修正方法对于损伤识别算法、响应、参数的研究较为系统.但大部分模型修正方法仅针对单一损伤类型识别,并且结构的损伤位置和程度仅通过数值显示,直观可视化效果较差[9].而弹性成像方法融合了先进成像方法和模型修正方法的优势,无需任何先验信息,可以从图形角度反映结构全局区域的损伤位置和程度,对结构损伤的几何描述具有较强的直观性.

拓扑优化[10-12]与弹性成像类似,其本质是通过将设计变量与弹性模量耦合实现对整体结构性能的调控.针对拓扑优化,单元密度为1 的单元对结构性能贡献最大,单元密度接近0 的单元对结构性能贡献最小.针对损伤识别,单元密度为1 的单元对应实体结构无损伤区域,单元密度小于1 的单元对应实体结构损伤区域.拓扑优化与损伤识别都是为了找到满足目标函数的弹性模量空间分布.因此,拓扑优化也可以用于结构损伤识别.从损伤表征层面,目前基于拓扑优化的损伤识别方法主要采用水平集法[13-14]、移动变形组件法[15-16]和变密度法[17].水平集法通过水平集函数隐式描述结构,其中零水平集表征结构损伤的边界[18].移动变形组件法通过组件的移动、旋转等方式实现结构的变化,采用拓扑描述函数描述结构,拓扑描述函数值大于0 表示有组件区域,拓扑描述函数值小于0 表示缺损区域[19].变密度法通过单元密度函数描述结构,单元密度为1 表示实体结构无损区域,单元密度接近0 表示实体结构缺损区域[20].水平集法对于损伤几何形状的描述能力强,但一般用于反演结构孔洞损伤的形状[21].对于基于水平集法或移动变形组件法的损伤识别方法,损伤数量需提前预估,初始孔洞或组件数量会影响到最终反演结果[22].从损伤识别目标函数构建层面,其构建可采用静力响应和动力响应数据.静力响应一般包括位移、应变和应力等.然而在实际工程问题中,损伤会导致结构刚度、质量发生改变,而结构的静力响应无法反映结构质量特性的变化,并且静力响应难以测量.但结构的动力响应容易通过传感器测量[23],能同时反映结构刚度、质量特性的变化,并且动力响应在结构发生损伤的早期阶段更加敏感.动力响应的响应类型一般分为频率[24]、振型[25]、频响函数[26]等.在静力响应和动力响应下,基于拓扑优化的损伤识别方法都表现出了较好的损伤识别效果[27].

上述研究对于损伤的表征比较单一,均只针对实体结构中造成结构刚度和质量减小的损伤类型(缺损、夹杂).然而在航空航天领域中,由于复杂自然环境[28-29]及各种意外情况[30]的影响,机械装备结构除会产生造成刚度和质量减小的损伤类型(缺损、夹杂)外,具有多孔特征的结构[31-32]还会产生结冰或积水损伤,这两类损伤类型分别会造成结构刚度和质量增大、结构仅质量增大.大多数损伤会造成结构弹性模量发生改变.因此,研究一种以拓扑优化驱动弹性成像的方法,实现结构缺损、夹杂、结冰、积水多种损伤类型的全局识别,对于丰富弹性成像理论内涵和实际工程应用具有重要意义.

本文提出了一种考虑刚度和质量耦合效应的结构弹性成像方法.首先,构建缺损、夹杂、结冰和积水多种损伤类型的统一表征.其次,构建损伤表征、力学模型和特征值响应之间的映射关系.再建立基于特征值响应的弹性成像模型,并建立弹性成像反演的数值实施方式.最后,分别考虑缺损、结冰、积水3 类损伤,探讨对二维结构进行弹性成像实现损伤识别的有效性,并探讨三维结构弹性成像问题验证本方法的通用性.

1 多种损伤类型的统一表征

为了解决结构缺损、夹杂、结冰和积水多种损伤类型的全局识别问题,关键在于建立不同损伤类型的统一表征.在损伤表征过程中,损伤会造成结构刚度和质量发生改变,并且不同损伤类型对结构刚度和质量的影响不同.例如,缺损和夹杂会同时造成结构刚度和质量减小,结冰会同时造成结构刚度和质量增大,积水不会影响结构刚度但会造成结构质量增大.弹性模量仅能反映结构刚度变化,单一的弹性模量不足以表征多种损伤类型,因而还需引入能反映结构质量变化的参数,构建同时考虑结构刚度和质量效应的损伤表征.受动力学拓扑优化理论启发,本文以结构离散单元的弹性模量系数和材料密度系数作为成像参数,将成像参数与弹性模量、材料密度关联,构建考虑结构刚度和质量耦合效应的损伤表征.结构的损伤类型可分为以下3 类:

(1)第1 类损伤类型(缺损、夹杂): 损伤同时造成结构的刚度和质量减小,继而导致成像参数减小;

(2)第2 类损伤类型(结冰): 损伤同时造成结构的刚度和质量增大,继而导致成像参数增大;

(3)第3 类损伤类型(积水): 损伤对结构刚度无影响,仅造成结构质量增大,继而导致成像参数增大.

基于以上分类,定义弹性模量系数x和材料密度系数y两种成像参数,建立成像参数与损伤表征的数学模型.如图1 所示为多种损伤类型表征,在同一套有限元网格定义了两套成像参数,第一套成像参数xj针对弹性模量,第二套成像参数yj针对材料密度,两套成像参数可反映不同损伤类型对结构刚度和质量的影响.

图1 多种损伤类型表征Fig.1 Characterization of damages

式中,x和y为成像参数,x和y最小值取0.001 为避免刚度和质量矩阵奇异.当无损结构为实体结构时,b取值为1;当无损结构为具有多孔特征的结构时,b取值范围为 0.001 ≤b＜1.

基于拓扑优化理论的固体各向同性材料惩罚模型(solid isotropic material with penalization,SIMP)[33],分别构建成像参数x(弹性模量系数)与弹性模量、成像参数y(材料密度系数)与材料密度的非线性插值模型如下

式中,E0为实体材料的弹性模量,E为单元弹性模量,ρ0为实体材料的密度,ρ 为单元材料密度.p和q为惩罚指数,其可加快优化模型收敛,然而在动力学问题中,p＞q这个条件会导致在拓扑优化过程中产生局部模态,有研究提出采用更大的质量惩罚指数 (p＜q) 可以控制局部模态[24],因此本文设置p＜q.

2 损伤对结构性能的影响

为了考虑损伤对结构性能的影响,将成像参数嵌入到无阻尼自由振动系统的有限元模型.

将结构离散为n个单元,则

式中,x,y为结构成像参数矩阵.

单元刚度、质量矩阵可写为

式中,Kj为第j个单元的刚度矩阵,B为应变位移矩阵,D为结构材料弹性矩阵,为实体材料单元刚度矩阵,Mj为第j个单元的质量矩阵,N为形函数矩阵,为实体材料单元质量矩阵,Ω 表示结构域,Ωj表示第j个单元的结构域.

全局刚度、质量矩阵可表示为

式中,K为关于成像参数x的结构全局刚度矩阵,M为关于成像参数y的结构全局质量矩阵.

进一步构建损伤表征、力学模型和特征值响应之间的映射关系,首先对结构无阻尼自由振动系统进行有限元分析,得到无阻尼自由振动系统的有限元平衡方程

式中,x,y为成像参数矩阵,λi为第i阶特征值,其中λi=ωi2,ωi为第i阶固有圆频率,ωi=2πFi,Fi为第i阶特征频率,φi为第i阶特征向量.

3 基于特征值响应的弹性成像建模与求解

3.1 弹性成像模型

特征值响应是最容易得到的模态参数之一,在建立弹性成像模型时,本文仅考虑结构特征值响应.当使用特征值作为响应参数进行损伤识别时有以下两个特点[35]: (1)当损伤出现在对称结构的对称位置时,可能会引起特征值响应的变化量相同;(2)结构各阶特征值响应对各处损伤的敏感性不同.因此,为避免损伤识别结果出现多解的问题,提高损伤识别精度,应适当选用多阶特征值构造弹性成像目标函数.本文以数字模型和真实结构模型的特征值响应变化率的平方和为目标函数,以有限元平衡方程和成像参数上下限为约束条件,构建如下结构弹性成像模型

式中,x,y分别为元素取值在之间的成像参数矩阵,xmin和xmax分别为成像参数xj可取的最小值和最大值,ymin和ymax分别为成像参数yj可取的最小值和最大值,其中xmin,xmax,ymin和ymax的取值由式(1) 可知,n为成像参数的数量,f(x,y) 表示弹性成像目标函数,m为模态阶数,wi为第i阶特征值的权重,权重总和为1,λi为数字模型的特征值响应,λi*为真实结构模型的特征值响应.特征值权重是可调参数,通过改变特征值权重可灵活调整各阶特征值对损伤识别的敏感性.

3.2 灵敏度分析

为了求解弹性成像模型,本文推导了目标函数关于成像参数的导数.

对目标函数关于单元成像参数xj,yj分别求导可得

对结构特征问题方程(11)关于单元成像参数xj和yj分别求导可得

式中,φi是正则化后的第i阶特征向量,T 代表矩阵转置.

(1)当考虑对刚度和质量有耦合效应的损伤类型(第1 类或第2 类损伤类型)时,由式(1)可知y=x.将y=x代入式(13)～式(19)得到目标函数关于单元成像参数的导数为

其中,特征值关于单元成像参数的导数为

再将式(21)代入式(20)得到最终目标函数关于单元成像参数的导数为

(2)当考虑仅对质量有影响的损伤类型(第3 类损伤类型)时,由式(1)可知x=b,因此无需对单元成像参数xj求导.将x=b代入式(13)～式(19)得到目标函数关于单元成像参数的导数为

其中,特征值关于单元成像参数的导数为

由于K为常数,因此可得

此时特征值关于单元成像参数的导数为

再将式(26)代入式(23)得到最终目标函数关于单元成像参数的导数为

4 数值实施

图2 所示为考虑刚度和质量耦合效应的结构弹性成像流程图.本方法需获得损伤结构的特征值响应,并作为唯一已知的输入数据,结构内部损伤的位置、大小和数量均未知.为此,本文采用数值实验替代实际实验获得损伤结构的特征值响应,通过在结构上开孔(预设缺损)或设置其他损伤,模拟可能存在的损伤结构,并通过有限元法计算损伤结构的特征值响应.然后求解数字模型(无损结构模型)的结构特征问题方程得到特征值和特征向量,并正则化特征向量用于之后灵敏度分析.之后,计算目标函数和进行灵敏度分析.再采用基于梯度的优化算法更新成像参数,通过式(42)的收敛准则判断目标函数是否收敛,优化模型不断迭代直至满足收敛准则,最后输出成像结果.

图2 考虑刚度和质量耦合效应的结构弹性成像流程Fig.2 Structural elastography process considering the coupling effect of stiffness and mass

本文采用基于梯度的优化算法,考虑成像参数的约束条件,成像参数x和y的迭代更新策略分别如下

式中,r为迭代步数,α 为迭代步长,min 与 max 分别代表取括号中最小值与最大值.

然后将第r步敏度和第r+1 步敏度相加取平均值,以保持迭代曲线的稳定性

最终成像参数x和y的迭代更新策略分别如下

(1)当结构存在对刚度和质量有耦合影响的损伤类型(第1 类或第2 类损伤类型)时,由式(1)可知y=x,显然此时h1=h2,因此成像参数的迭代更新策略可统一为

(2)当结构存在仅对质量有影响的损伤类型(第3 类损伤类型)时,由式(1)可知x为常数,因此x无需迭代更新.此时成像参数的迭代更新策略为

定义结构弹性成像模型的收敛准则如下

式中,ε 代表一个极小正数,rmax为最大迭代步数.

5 数值算例

本节通过4 个数值算例来验证考虑刚度和质量耦合效应的结构弹性成像方法的有效性和通用性.在本节所有算例中,二维/三维实体结构的弹性模量E0=200 GPa,材料密度 ρ0=7800 kg/m3,泊松比ν=0.3(如无特殊说明,本文仅考虑各向同性材料),弹性成像模型的成像参数初始值为无损结构成像参数值.为了解决难以区分对称结构对称位置损伤的问题,需要在结构特定位置施加合适的非结构集中质量[36].

5.1 存在第1 类损伤类型的结构弹性成像

为验证本方法识别结构中第1 类损伤类型(缺损)的有效性,定义图3 所示的二维半MBB (Messerschmitt-Bölkow-Blohm)梁结构,初始无损结构为实体结构,结构域长宽比为2:1,长度为4 m,宽度为2 m,结构左端支撑约束,右下角支撑约束,结构被离散为80×40 个双线性四边形单元,在结构长度方向5/8 位置及宽度方向1/2 位置设置一个尺寸为0.6 m×0.6 m的缺损区域,缺损区域材料弹性模量和材料密度分别设置为0.001E0,0.001ρ0,在结构右下角位置节点上施加数值大小为结构域质量20%的非结构集中质量.本算例仅考虑结构前6 阶特征值,特征值权重w=[0.25 0.15 0.15 0.15 0.15 0.15].

图3 内置缺损的二维半MBB 梁结构Fig.3 Two-dimensional half-MBB beam structure with built-in defect

目标函数迭代曲线如图4 所示,随着优化的进行,目标函数值下降速度由快到慢,当弹性成像模型迭代287 步后目标函数值趋于平稳,结构成像图中低成像参数值区域逐渐清晰.从图4 中可以看出迭代曲线光滑,表明算法稳定性良好.在弹性成像模型经过616 次迭代后,得到如表1 所示的无损模型的特征频率,并得到如图5 所示的A 缺损结构的成像结果,成像结果下方颜色条的刻度值对应成像参数值.由表1 可知最终无损模型的前6 阶特征频率与损伤模型的前6 阶特征频率基本一致.为了方便将损伤识别结果与真实损伤区域进行对比,图5 中黑色正方形边线表示真实损伤结构的损伤区域边界.由图5 可知,最终成像结果显示的损伤区域位置和大小与设定损伤区域基本一致.由于本文是基于变密度拓扑优化法构建了成像参数与弹性模量、材料密度的非线性插值模型,因而成像参数变量存在中间值,导致成像结果的边界会存在模糊的过渡,故本方法对损伤形状的成像是一种近似,无法严格识别出损伤的精准形状.由于噪声数据影响导致成像结果左下方出现些许高亮区域,但从成像结果中可以直观看出高亮区域的成像参数值与损伤区域的成像参数值明显相差较大,因此可以辨别出高亮区域不是损伤区域.

表1 算例1 中损伤模型和无损模型的特征频率(Hz)Table 1 Eigen-frequencies of damaged model and undamaged model in the example 1 (Hz)

图4 二维半MBB 梁结构弹性成像迭代曲线Fig.4 Iterative curve of elastography of two-dimensional half-MBB beam structure

图5 二维半MBB 梁结构成像结果Fig.5 Two-dimensional half-MBB beam structure imaging results

5.2 存在第2 类损伤类型的结构弹性成像

为验证本方法识别结构中第2 类损伤类型(结冰)的有效性,定义如图6 所示的二维悬臂梁结构,初始无损结构为具有多孔特征的结构,该结构由80×40 个多孔单胞构成,结构的弹性模量和材料密度分别设置为0.3E0,0.3ρ0,结构域长度为4 m,宽度为2 m,结构左端固定约束,结构被离散为80×40 个双线性四边形单元,在结构长度方向1/4 位置及宽度方向1/2 位置设置一个尺寸为0.6 m×0.6 m 的结冰区域,结冰区域的弹性模量和材料密度分别设置为E0,ρ0,在结构右上角位置节点上施加数值大小为结构域质量20%的非结构集中质量.本算例考虑与算例5.1 相同的特征值阶次和特征值权重.

图6 内置结冰的二维悬臂梁结构Fig.6 Two-dimensional cantilever beam structure with built-in ice

如图7 所示,随着优化的进行,目标函数值下降速度很快,弹性成像模型迭代157 步后目标函数开始趋于平稳,迭代900 步后满足收敛准则,得到如表2 所示的无损模型的特征频率,并得到如图8 所示的B 结冰结构的成像结果.由表2 可知最终无损模型的前6 阶特征频率与损伤模型的前6 阶特征频率基本一致.成像结果中出现一个高成像参数值的区域,该区域为识别的结冰区域,与设定的结冰区域位置和大小近乎完全一致,表现出了良好的弹性成像效果.本算例采用与算例5.1 不同的边界条件,结果表明在不同边界条件下本方法仍具有良好的损伤识别效果.

表2 算例2 中损伤模型和无损模型的特征频率(Hz)Table 2 Eigen-frequencies of damaged model and undamaged model in the example 2 (Hz)

图7 具有多孔特征的二维悬臂梁结构弹性成像迭代曲线Fig.7 Iterative curve of elastography of two-dimensional cantilever beam structure with porous features

图8 具有多孔特征的二维悬臂梁结构成像结果Fig.8 Imaging results of two-dimensional cantilever beam structure with porous features

5.3 存在第3 类损伤类型的结构弹性成像

为进一步验证本方法识别结构第3 类损伤类型(积水)的有效性,如图9 所示初始无损结构为具有多孔特征的结构,该结构由60×20 个多孔单胞构成,结构弹性模量和材料密度分别设置为0.3E0,0.3ρ0,结构长度和宽度分别为3 m,1 m,采用与算例5.1 相同的边界条件,结构被离散为60×20 个双线性四边形单元,在结构长度方向3/10 位置和宽度方向3/4 位置、长度方向3/4 位置和宽度方向1/4 位置分别设置尺寸为0.2 m×0.2 m 的积水区域,由于积水仅影响结构质量,因此积水区域的材料弹性模量不变,材料密度设置为 ρ0.本算例施加非结构集中质量的位置和数值大小与算例5.1 保持一致,并且考虑与算例5.1 相同的特征值阶次和特征值权重.

图9 内置积水的二维半MBB 梁结构Fig.9 Two-dimensional half-MBB beam structure with built-in water accumulation

由图10 可知,在弹性成像模型最初的几步迭代中目标函数值迅速下降,然后目标函数值在弹性成像模型接下来的迭代步中下降变缓,最后目标函数值在弹性成像模型迭代165 步后趋于平稳.弹性成像模型迭代389 步后得到如表3 所示的无损模型的特征频率,由表3 可知最终无损模型的前6 阶特征频率与损伤模型的前6 阶特征频率基本一致.图11为弹性成像模型迭代389 步后结构的成像结果,从图中可以看出,最终成像结果清晰,可以几乎无偏差地识别出两个积水区域各自的位置和大小,说明本方法用于多损伤条件下也具有良好的损伤识别效果.

表3 算例3 中损伤模型和无损模型的特征频率(Hz)Table 3 Eigen-frequencies of damaged model and undamaged model in the example 3 (Hz)

图10 具有多孔特征的二维半MBB 梁结构弹性成像迭代曲线Fig.10 Iterative curve of elastography of two-dimensional half-MBB beam structure with porous features

图11 具有多孔特征的二维半MBB 梁结构成像结果Fig.11 Imaging results of two-dimensional half-MBB beam structure with porous features

5.4 三维结构弹性成像

为进一步验证本方法的通用性,开展三维结构弹性成像的研究.如图12 所示的三维悬臂梁结构,结构长为2 m、宽为0.6 m、高为1.2 m,结构左端面施加固定约束,结构被离散为20×6×12 个六面体单元,在结构内部中心位置设置一个尺寸为0.4 m×0.4 m×0.4 m 的损伤区域.分别考虑损伤区域为缺损、结冰、积水三类损伤类型,定义如下3 种存在损伤的三维悬臂梁结构: (1)设置初始无损结构为实体结构,损伤区域为缺损,缺损区域材料弹性模量设置为0.001E0,材料密度设置为0.001ρ0;(2)设置初始无损结构为具有多孔特征的结构,该结构由20×6×12 个多孔单胞构成,结构弹性模量和材料密度分别设置为0.5E0和0.5ρ0,损伤区域为结冰,结冰区域的弹性模量和材料密度分别设置为E0和ρ0;(3)同样设置初始无损结构为具有多孔特征的结构,该结构也由20×6×12 个多孔单胞构成,结构弹性模量和材料密度也分别设置为0.5E0和0.5ρ0,但损伤区域为积水,由于积水仅影响结构质量,因此积水区域的材料弹性模量不变,材料密度设置为ρ0.本算例在结构上端面与后端面交线长度方向3/4 位置的节点上施加数值大小为结构域质量20%的非结构集中质量.三维结构相较于二维结构的模型更加复杂,逆向反演三维结构的损伤也更为复杂,为了提高损伤识别精度,本三维算例选取的特征值阶数比二维算例略多.本算例考虑结构前10 阶特征值,每阶特征值权重均设置为0.1.

图12 内置损伤的三维悬臂梁结构Fig.12 Three-dimensional cantilever beam structure with built-in damaged

如图13 所示,弹性成像迭代曲线光滑且目标函数最终收敛,弹性成像模型经过680 次迭代后得到如图14 所示的F 缺损结构的成像结果.如图15 所示,随着优化的进行,目标函数值下降速度由快到慢并且很快收敛,经过531 次迭代后得到如图16 所示的F 结冰结构的成像结果.如图17 所示,弹性成像模型迭代197 步后目标函数开始趋于平稳,经过256 次迭代后得到如图18 所示的F 积水结构的成像结果.由于损伤位于结构内部导致不方便观察,因此通过显示半边结构空间的成像结果来解决这一问题.从3 个三维结构成像结果图中可以看出,最终成像结果所呈现的损伤区域位置及大小与设定损伤区域基本一致,表明本方法用于三维结构同样具有良好的损伤识别效果.虽然三维结构弹性成像问题相较于二维结构弹性成像问题计算成本略高,但是三维成像只需定义新的三维结构无阻尼自由振动系统有限元模型,损伤表征、弹性成像模型和求解算法均与二维成像保持一致.