基于拓扑优化的结构弹性成像方法1)

付君健 *, 李帅虎 李 好 ** 高 亮 ** 周祥曼 *, 田启华 *,

* (水电机械设备设计与维护湖北省重点实验室,湖北宜昌 443002)

† (三峡大学机械与动力学院,湖北宜昌 443002)

** (华中科技大学数字制造装备与技术国家重点实验室,武汉 430074)

引言

弹性成像源自医学影像领域,是一种通过弹性模量表征生物体组织物理特性的成像方法[1-2].弹性成像能够获得常规成像技术难以获取的弹性模量信息,辅助临床诊断生物体器官病理变化的位置和形状.对于工程中的机械装备结构,结构损伤会造成特定部位弹性模量发生改变,继而影响结构的整体物理特性.将弹性模量与损伤表征关联,发展结构缺陷识别的弹性成像技术,对于机械工程装备的服役安全评价具有重要意义.

与超声成像[3]、电阻抗成像[4-5]、热成像[6-7]、太赫兹成像[8-9]等先进的成像方法类似,弹性成像是电子技术、计算机技术和反演数学相结合的产物,它利用计算机对穿越物体的物理信号进行处理,通过算法重建物体内部信息.弹性成像的基本原理是对生物体组织施加一个外部静态或动态激励.根据胡克定律,生物体组织在外力作用下将产生位移、应变等物理响应,结合超声波测距法和数字信号处理技术,可以预测生物体内部组织的弹性模量信息.尽管超声成像、X 射线等无损检测技术可实现结构的缺陷和损伤识别,但多数检测技术属于局部识别技术,对环境要求高,需要先验信息,难以检测整体结构性能缺陷.对于大型机械装备结构,损伤一般发生在结构内部或难以观测位置,发展全局性的损伤识别技术具有重要的工程意义.作为一种全局性的损伤识别方法,弹性成像是成像参数、外界激励和反演模型三大要素的综合表达.弹性成像的参数是材料的弹性模量,一般情况下弹性阶段的力学响应即可实现成像,如结构静态位移[10]和动力响应[11],不必拓展至塑性分析阶段.通过对比仿真响应和实测响应,可以建立弹性成像的反演模型.弹性成像属于约束优化问题,一般作为反问题求解,其关键技术是弹性成像的反演算法.

拓扑优化是一种在给定设计空间和约束下,寻求满足性能最优的材料分布形式的结构优化方法[12-14].拓扑优化与弹性成像均属于约束优化问题,其本质都是为了寻找一种满足目标函数的弹性模量分布.因此,拓扑优化也可用于弹性成像,逆向反演结构的内部损伤所对应的弹性模量.在进行结构损伤识别时,由于拓扑优化几何描述能力强,不需要预设损伤形状,可自动反演复杂的损伤形状.例如对于带有孔洞缺陷的结构,采用水平集法[15-17]或移动变形组件法[18]及类似方法[19],可有效描述缺陷几何形状.虽然在静力响应和动力响应下,拓扑优化方法都表现出了较好的损伤识别效果[10,20],但仍存在两种问题需要解决.其一,水平集法具有清晰几何边界,一般用于带有孔洞缺陷的结构损伤识别.当结构存在多种不同弹性模量的材料时,或结构存在缺损、夹杂、裂纹等多种损伤形式,水平集法则难以对结构缺陷进行有效识别,需要在拓扑优化框架中同时引入几何形状和弹性模量作为反演变量[21-22].其二,目前采用水平集法或移动变形组件法进行结构损伤识别时,虽然无需预设缺陷形状,但需预估缺陷数量,初始孔洞或初始组件数量对反演结果有重要影响[23].因此,研究适用于多损伤并存、初始缺陷不敏感的弹性成像方法,提高弹性成像的局部表征和全局识别能力,具有更广泛的理论和工程意义.

由于弹性成像的数学本质是优化驱动的反问题,与结构拓扑优化方法数学原理相同,本文提出了一种基于拓扑优化的弹性成像方法.通过构建损伤表征、结构模型与物理响应的映射关系,建立并求解结构弹性成像的优化模型,并探讨了在不同结构维度、损伤类型、缺陷数量下结构弹性成像的效果.

1 弹性属性与缺陷表征

拓扑优化本质是寻求满足设计需求的最优材料分布,基于人工材料插值模型将单元密度与弹性模量耦合,调控单元弹性模量对结构性能的贡献[24-26],这种对弹性模量的全局调控机制与结构弹性成像技术不谋而合.弹性模量可以衡量物体抵抗弹性变形的能力,结构缺陷会造成特定部位材料弹性模量发生改变,继而影响结构的整体物理特性.若将弹性模量与损伤表征关联,则有望通过弹性模量的空间分布识别结构多种损伤,反映损伤对结构物理特性的综合影响.

为了解决结构变形、裂纹、缺损等多损伤的识别问题,不同损伤类型的统一表征是关键.受拓扑优化理论启发,采用结构离散单元的相对密度(弹性模量系数)作为弹性成像参数来表征损伤程度,建立结构缺陷、弹性属性和力学模型之间的耦合关系,成像参数与损伤表征的数学定义见式(1).如图1 所示,块状分布的中间弹性模量值可表征夹杂,带状分布的中/低弹性模量值可表征裂纹,块状分布的低弹性模量值可表征孔洞.通过缺陷的形态和成像参数的数值,即可实现对结构缺陷的评估.

图1 损伤表征Fig.1 Damage characterization

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

式中,x为成像参数,取值范围为 0 .001 ≤x≤1 ,x最小值取0.001 为避免刚度矩阵奇异.p为惩罚系数,其作用是为了加速优化模型收敛.E0为实体材料的弹性模量,E为被惩罚后的单元弹性模量.

进一步构建弹性属性和力学模型之间的关系,将融合了结构缺陷信息的弹性模量引入结构有限元平衡方程,以结构静力学问题为例,有限元平衡方程如下

式中,K为关于成像参数x的结构全局刚度矩阵,U为结构的全局位移响应,F为激励载荷向量,xe为单元e的成像参数,ke为单元的刚度矩阵,B为应变位移矩阵,D为实体材料弹性矩阵,Ω 代表结构域,Ωe代表单元e的结构域.通过求解式(3),可得到受到缺陷影响的结构位移响应,继而构建了损伤表征、结构模型与物理响应的映射关系.

2 弹性成像优化模型

为体现缺陷对结构物理响应的影响,以损伤结构响应和无损结构响应的最小二乘为目标函数.再定义与结构自由度维度相同的测量位置向量 ψa,通过在最小二乘函数前乘以测量位置向量,可实现结构特定位置响应误差的最小化

式中,J(x) 为关于成像参数的目标函数,ψa为结构响应的测量位置向量(在自由度a处元素为1,其余元素为0),a代表结构响应测量位置的自由度编号,上标 T 代表矩阵转置,U0为无损结构响应,U为含缺陷信息的损伤结构响应.

考虑到结构自由度上的响应值的数量级可能不同,以结构响应变化率的平方为目标函数,定义如下

在研究静力学结构弹性成像理论时,仅考虑结构位移响应,以成像模型和真实模型位移响应变化率的平方为目标函数,以静力学平衡方程为约束条件,以成像参数上下限为约束条件,构建式(8)所示的结构弹性成像模型

式中,x为元素取值在[0.001,1]之间的成像参数向量,N代表成像参数的数量,J(x) 表示目标函数,U0为初始结构(无损结构)位移响应,U为含缺陷信息的结构(损伤结构)位移响应,K为结构全局刚度矩阵,F为全局载荷向量.

从优化模型(8)可知,传统拓扑优化方法属于正向设计,即给定设计目标计算新型结构形式[28-31];而弹性成像属于逆向设计,通过已知结构响应反演结构的拓扑构型(弹性模量分布),从而实现结构缺陷形态和位置、损伤程度的识别.

3 优化模型求解

3.1 敏度分析

为了求解弹性成像优化模型,本文基于伴随法推导目标函数关于成像参数的导数[32],构造拉格朗日函数如下

式中,λ 为静力学平衡方程的拉格朗日乘子列向量.

求解拉格朗日乘子列向量得

则拉格朗日函数关于成像参数的导数为

3.2 数值实施

考虑成像参数的有界性约束,成像参数的迭代更新策略如下

为了增加收敛曲线的数值稳定性,将第n步敏度和n+1步敏度加权平均,以修正灵敏度信息

经过敏度修正后的成像参数更新策略如下

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

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

4 数值算例

4.1 二维均质材料结构弹性成像

为验证弹性成像的有效性,首先对二维均质材料结构进行弹性成像研究.如图2 所示的二维悬臂梁结构,结构域长宽为 2 m×1 m,悬臂梁左端固定,右端中心处施加有F=-1 N 的集中载荷,将结构域划分为 4 0×20 个双线性四边形单元,材料弹性模量E=100 GPa ,泊松比v=0.3 (如无特殊说明,本文仅考虑各向同性材料) .本算例的位移响应测量点设置在悬臂梁上、下边缘中点处.

图2 二维均质悬臂梁结构Fig.2 Two-dimensional homogeneous cantilever beam structure

对于均质各向同性材料结构,因为只存在唯一成像参数,需对敏度信息进行均匀化处理

为了得到响应点的位移响应,图3 所示,本文仅通过有限元法进行结构变形的数值模拟,得到两个测量点X和Y方向的位移响应作为模型输入.优化模型的初始弹性模量设置为 2 00 GPa,成像参数初始值为1,迭代步长为 ξ =0.1 .由图4 可知,经过280 步迭代后,目标函数收敛至6.594 × 10-5,成像参数收敛至0.5.将成像参数与初始弹性模量相乘,得到此时的结构弹性模量约为 1 00 GPa,反演得到的结构弹性模量与设置的弹性模量一致.对于均质材料结构,结构响应与成像参数之间呈线型关系,即使以少量测量点的位移响应为参照,也能准确反演整体结构的弹性模量,此方法可在结构材料未知的情况下对弹性模量进行成像.

图3 二维均质悬臂梁弹性成像结构Fig.3 Two-dimensional homogeneous cantilever beam elastography structure

图4 二维均质悬臂梁弹性成像迭代图Fig.4 Iterative graph of two-dimensional homogeneous cantilever beam elastography

4.2 二维非均质材料结构弹性成像

为验证本方法对非均质材料结构成像的有效性,如图5 所示,采用算例4.1 中的结构参数和边界条件,材料弹性模量设置为E=200 GPa,泊松比均v=0.3 ,将结构域划分为4 0×20 个双线性四边形单元.弹性成像的最小分辨率为1 个单元尺寸,在悬臂梁长度方向1/4 和3/4 的位置及宽度方向1/2 位置设置尺寸为 0 .2 m×0.2 m 的孔洞缺陷,孔洞缺陷区域材料的弹性模量设置为 0 .001E,迭代步长为 ξ =0.1,分别设置单孔洞缺陷和双孔洞缺陷形式,进行非均质结构的弹性成像.本算例的位移响应测量点设置在悬臂梁上、下边缘处.

图5 内置缺陷的二维悬臂梁结构Fig.5 Two-dimensional cantilever beam structure with built-in damaged

图6(a)所示为存在A缺陷结构的成像结果,优化模型迭代259 步后目标函数收敛,结构内部出现较低的弹性模量区域,结构弹性成像结果较清晰,最终弹性成像呈现的缺陷位置和大小与设定缺陷基本一致.图6(b)所示为存在B缺陷结构的成像结果,优化模型迭代93 步后目标函数收敛,结构内部缺陷位置的弹性模量与设定缺陷有一定误差.虽然B缺陷弹性成像结果较模糊,但已初步表明缺陷存在的大致位置.图7 所示为同时考虑A和B两处缺陷的结构成像过程.在优化模型迭代至100 步时,左边A缺陷的弹性模量先识别出来,右边B缺陷的弹性模量识别稍慢.随着迭代的进行,A缺陷的弹性模量成像逐渐清晰,可正确识别A缺陷的位置和大小,B缺陷弹性模量成像的清晰度一般,但仍可以表征此处存在一定缺陷,成像结果与A、B缺陷单独成像的结果基本一致.在前述研究中[10,20-21,23],若结构存在多个缺陷,在算法识别时需要提前设定缺陷数量.而本文在单元上定义了成像参数,成像参数的数量与单元数量相同,每一处的单元都可用以表征缺陷.因此,本文方法无需提前人为指定缺陷数量,算法可自动识别缺陷处的弹性模量,通过弹性模量的变化判断是否存在缺陷.

图6 二维悬臂梁结构单缺陷弹性成像结果Fig.6 Single damage elastography results of the two-dimensional cantilever beam structure

图7 二维悬臂梁结构多缺陷弹性成像过程Fig.7 Multi-damage elastography process of the two-dimensional cantilever beam structure

由图6～ 图8 可知,不同位置缺陷的弹性成像质量有一定差异,主要存在两方面的原因.其一,成像质量与缺陷对结构性能的影响程度有关.以静力学为例,不同位置的材料对测量点处结构变形的贡献度不同,若缺陷对测量点处的结构变形影响较小(如缺陷B),则弹性成像效果会偏弱,若缺陷对结构测量点处的变形影响较大,则成像质量效果会偏强(如缺陷A).在本例中若将约束置于右边界,载荷置于左边界中点处,则B缺陷的识别效果会更好.其二,弹性成像质量与测量点数量有关.不同于均质材料结构的弹性成像,非均匀材料结构的结构响应与成像参数之间呈现高度非线性,采用较少的测量点无法准确对含有缺陷的结构进行弹性成像.本算例仅采用悬臂梁上、下边界来测量结构响应,没有用到结构内部的测量点.在实际工程中,测量点的数量有限、位置特殊,有限的测量点会一定程度影响弹性成像的结果.

图8 二维悬臂梁结构弹性成像迭代图Fig.8 Iterative diagram of elastography of two-dimensional cantilever beam structure

4.3 边界条件对成像结果的影响

为验证弹性成像方法的通用性,研究边界条件对弹性成像的影响.定义图9 所示的二维米歇尔(Michell)梁结构,结构域长宽为 2 m×1 m,米歇尔梁左下角固定约束,右下角支撑约束.在下端中心处施加F=-1 N 的集中载荷,将结构域划分为 4 0×20 个双线性四边形单元.材料弹性模量E=200 GPa,泊松比v=0.3 .在米歇尔梁长度方向1/4 和3/4 的位置及宽度方向1/2 位置设置尺寸为 0 .2 m×0.2 m 的孔洞缺陷,孔洞缺陷区域材料的弹性模量设置为 0 .001E,分别设置单孔洞缺陷和双孔洞缺陷形式,进行不同边界条件的非均质结构的弹性成像.

图9 内置缺陷的二维Michell 梁结构Fig.9 Two-dimensional Michell-like beam structure with built-in damaged

对比图6 和图10 可知,同样的二维长方形结构,在米歇尔梁的约束条件下,依然可以得到缺陷A和缺陷B的弹性模量图,说明本方法在不同的边界条件下本方法仍具有良好的成像特性.在不同的边界条件下,弹性成像结果会存在微小的差异,因为同样的缺陷在不同边界条件下,对结构性能的影响不同,进而会影响弹性成像优化模型的反演.

图10 二维Michell 梁结构弹性成像结果Fig.10 Two-dimensional Michell-like beam structure elastography results

4.4 三维结构弹性成像

为进一步验证本方法对非均质材料结构成像的通用性,对三维非均质材料结构进行弹性成像研究.如图11 所示的三维悬臂梁结构,结构域长、宽、高为L×H×W=2 m×1 m×1 m,悬臂梁左端固定,右端面中心处沿Y轴方向施加F=-1 N 的线载荷,将结构域划分为 4 0×20×20 个六面体单元,材料弹性模量E=200 GPa ,泊松比v=0.3 .在悬臂梁中心设置尺

图11 内置缺陷的三维悬臂梁结构Fig.11 Three-dimensional cantilever beam structure with built-in damaged

寸为 0 .2 m×0.2 m×0.2 m 的缺陷,考虑孔洞和夹杂两种缺陷形式,孔洞缺陷弹性模量设置为 0 .001E,夹杂缺陷弹性模量设置为 0 .5E.在悬臂梁上、下边缘处设置位移响应测量点,进行三维非均质结构的弹性成像.

以悬臂梁上、下端面节点的X和Y自由度位移响应为参照.图12(a)所示为夹杂缺陷的三维弹性成像结果,缺陷处的成像参数接近0.5.图12(b)所示为孔洞缺陷的三维弹性成像结果,缺陷处的成像参数接近0.由于缺陷位于结构内部,图12 仅显示半边结构空间的成像结果以方便观察.图13 所示优化迭代曲线光滑,表明算法具有良好的稳定性.由图12 可知,本方法可对不同损伤程度的缺陷形式进行弹性成像,不仅能识别缺陷的位置和相对大小,还能识别缺陷的近似弹性模量值.由二维弹性成像扩展至三维弹性成像时,仅需要定义新的三维结构有限元模型.三维弹性成像计算成本高于二维问题,弹性成像的变量定义、优化模型和求解算法均与二维弹性成像一致.

图12 三维悬臂梁结构弹性成像结果Fig.12 Three-dimensional cantilever beam structure elasticity elastography

图13 三维悬臂梁结构弹性成像迭代图Fig.13 Three-dimensional cantilever beam structure elastography iteration diagram

5 结论

本文提出一种基于拓扑优化的结构弹性成像方法,构建了损伤表征、结构模型与物理响应的映射关系,建立并求解了结构弹性成像的优化模型.研究结果表明,弹性成像方法在二维、三维问题上具有通用性,可有效实现均质/非均质、单/多缺陷结构的弹性成像,且弹性成像结果不依赖于特定边界条件.本文虽然仅提出了静态响应下的弹性成像理论及求解方法,但该方法可进一步拓展至复杂动态响应的结构弹性成像.