EI、Scopus 收录
中文核心期刊

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

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

付君健, 彭铁川, 李然, 李帅虎, 周祥曼, 李响. 考虑刚度和质量耦合效应的结构弹性成像方法. 力学学报, 2023, 55(9): 1971-1982. DOI: 10.6052/0459-1879-23-220
引用本文: 付君健, 彭铁川, 李然, 李帅虎, 周祥曼, 李响. 考虑刚度和质量耦合效应的结构弹性成像方法. 力学学报, 2023, 55(9): 1971-1982. DOI: 10.6052/0459-1879-23-220
Fu Junjian, Peng Tiechuan, Li Ran, Li Shuaihu, Zhou Xiangman, Li Xiang. Structural elastography method considering the coupling effect of stiffness and mass. Chinese Journal of Theoretical and Applied Mechanics, 2023, 55(9): 1971-1982. DOI: 10.6052/0459-1879-23-220
Citation: Fu Junjian, Peng Tiechuan, Li Ran, Li Shuaihu, Zhou Xiangman, Li Xiang. Structural elastography method considering the coupling effect of stiffness and mass. Chinese Journal of Theoretical and Applied Mechanics, 2023, 55(9): 1971-1982. DOI: 10.6052/0459-1879-23-220
付君健, 彭铁川, 李然, 李帅虎, 周祥曼, 李响. 考虑刚度和质量耦合效应的结构弹性成像方法. 力学学报, 2023, 55(9): 1971-1982. CSTR: 32045.14.0459-1879-23-220
引用本文: 付君健, 彭铁川, 李然, 李帅虎, 周祥曼, 李响. 考虑刚度和质量耦合效应的结构弹性成像方法. 力学学报, 2023, 55(9): 1971-1982. CSTR: 32045.14.0459-1879-23-220
Fu Junjian, Peng Tiechuan, Li Ran, Li Shuaihu, Zhou Xiangman, Li Xiang. Structural elastography method considering the coupling effect of stiffness and mass. Chinese Journal of Theoretical and Applied Mechanics, 2023, 55(9): 1971-1982. CSTR: 32045.14.0459-1879-23-220
Citation: Fu Junjian, Peng Tiechuan, Li Ran, Li Shuaihu, Zhou Xiangman, Li Xiang. Structural elastography method considering the coupling effect of stiffness and mass. Chinese Journal of Theoretical and Applied Mechanics, 2023, 55(9): 1971-1982. CSTR: 32045.14.0459-1879-23-220

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

基金项目: 国家自然科学基金(51775308)和湖北省自然科学基金(2021CFB236)资助项目
详细信息
    通讯作者:

    李然, 正高级工程师, 主要研究方向为大型金属结构损伤识别与寿命评判、内河通航枢纽工程技术. E-mail: rangly@126.com

  • 中图分类号: O342

STRUCTURAL ELASTOGRAPHY METHOD CONSIDERING THE COUPLING EFFECT OF STIFFNESS AND MASS

  • 摘要: 弹性成像是一种医学影像技术, 可以把生物体组织的弹性模量信息转换为可视化图像. 为了将弹性成像用于机械结构的损伤识别, 实现机械结构缺损、夹杂、结冰和积水多种损伤类型的全局识别, 提出一种考虑刚度和质量耦合效应的结构弹性成像方法. 弹性模量仅能反映结构刚度变化, 单一的弹性模量不足以表征多种损伤类型, 因而还需引入能反映结构质量变化的参数. 受结构拓扑优化理论启发, 以结构离散单元的弹性模量系数和材料密度系数作为成像参数, 将成像参数与弹性模量和材料密度关联, 构建考虑结构刚度和质量耦合效应的损伤表征. 基于无阻尼自由振动系统的有限元模型, 构建损伤表征、力学模型和特征值响应之间的映射关系. 以数字模型和真实结构模型的特征值响应变化率的平方和为目标函数, 以有限元平衡方程和成像参数上下限为约束条件, 建立基于特征值响应的结构弹性成像模型, 推导弹性成像目标函数关于成像参数的导数, 采用基于梯度的优化算法求解弹性成像模型. 数值算例表明, 针对结构缺损、夹杂、结冰、积水等损伤, 该方法无需先验信息, 可有效实现结构损伤位置、数量和形状的精确量化, 并通过三维结构弹性成像算例验证了该方法的通用性.
    Abstract: Elastography is a medical imaging technology that can convert the elastic modulus information of biological tissues into visual images. In order to use elastography for damage identification of mechanical structures and achieve global identification of various damage types such as mechanical structural defects, inclusions, ice, and water accumulation. A structural elastography method considering the coupling effect of both stiffness and mass is proposed. The elastic modulus can only reflect the change in structural stiffness and is not enough to characterize various damage types. It is also necessary to introduce parameters that reflect the change in the mass of the structure. Inspired by the structural topology optimization theory, the proposed method takes the elastic modulus coefficients and material density coefficients of the discrete elements of the structure as the imaging parameters. The imaging parameters are correlated with the elastic modulus and material density to construct a damage characterization considering the coupling effect of structural stiffness and mass. The mapping relationship between damage characterization, mechanical model, and eigenvalue response is constructed based on the finite element model of the undamped free vibration system. The sum of the squares of the eigenvalue response change rates of the digital model and the real structural model is used as the objective function. The finite element equilibrium equation as well as the upper and lower limits of the imaging parameters are used as constraints. Eigenvalue response-based structural elastography model is established based on the objective function and constraints. The derivatives of the elastography objective function with respect to the imaging parameters are derived. The elastography model is solved using a gradient-based optimization algorithm. Numerical examples show that this method can effectively achieve accurate quantification of structural damage location, quantity, and shape without any prior information for damages such as structural defects, inclusions, ice, and water accumulation. The generality of this method is further verified by three-dimensional structural elastography examples.
  • 弹性成像[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类损伤, 探讨对二维结构进行弹性成像实现损伤识别的有效性, 并探讨三维结构弹性成像问题验证本方法的通用性.

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

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

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

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

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

    图  1  多种损伤类型表征
    Figure  1.  Characterization of damages
    $$ \left.\begin{split} &x = y = b\text{ }无损伤\\ &x = y = 0.001\text{ }缺损\\ &0.001 < x = y < b\text{ }夹杂\\ &b < x = y\leqslant 1\text{ }结冰\\ &x = b,b < y\leqslant 1\text{ }积水\end{split}\right\} $$ (1)

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

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

    $$ E(x) = {x^p}{E^0} $$ (2)
    $$ \rho (y) = {y^q}{\rho ^0} $$ (3)

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

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

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

    $$ {\boldsymbol{x}} = [{x_1}{\text{ }}{x_2}{\text{ }}\cdots{\text{ }}{x_n}] $$ (4)
    $$ {\boldsymbol{y}} = [{y_1}{\text{ }}{y_2}{\text{ }}\cdots{\text{ }}{y_n}] $$ (5)

    式中, ${\boldsymbol{x}}$, ${\boldsymbol{y}}$为结构成像参数矩阵.

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

    $$ {{\boldsymbol{K}}_j}({x_j}) = \int_{{\varOmega ^j}} {{{\boldsymbol{B}}^{\text{T}}}{\boldsymbol{D}}({x_j}){\boldsymbol{B}}{\text{d}}\varOmega } = {x_j^p}{\boldsymbol{K}}_j^0 $$ (6)
    $$ {{\boldsymbol{M}}_j}({y_j}) = \int_{{\varOmega ^j}} {\rho ({y_j}){{\boldsymbol{N}}^{\text{T}}}{\boldsymbol{N}}{\text{d}}\varOmega } = {y_j^q}{\boldsymbol{M}}_j^0 $$ (7)

    式中, $ {{\boldsymbol{K}}_j} $为第$j$个单元的刚度矩阵, $ {\boldsymbol{B}} $为应变位移矩阵, $ {\boldsymbol{D}} $为结构材料弹性矩阵, $ {\boldsymbol{K}}_j^0 $为实体材料单元刚度矩阵, $ {{\boldsymbol{M}}_j} $为第$j$个单元的质量矩阵, $ {\boldsymbol{N}} $为形函数矩阵, ${\boldsymbol{M}}_j^0$为实体材料单元质量矩阵, $ \varOmega $表示结构域, $ {\varOmega ^j} $表示第$j$个单元的结构域.

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

    $$ {\boldsymbol{K}}({\boldsymbol{x}}) = \sum\limits_{j = 1}^n {{x_j^p}{\boldsymbol{K}}_j^0} $$ (8)
    $$ {\boldsymbol{M}}({\boldsymbol{y}}) = \sum\limits_{j = 1}^n {{y_j^q}{\boldsymbol{M}}_j^0} $$ (9)

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

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

    $$ {\boldsymbol{M}}\left( {\boldsymbol{y}} \right){\boldsymbol{\ddot U}}{\text{ + }}{\boldsymbol{K}}\left( {\boldsymbol{x}} \right){\boldsymbol{U}}{{ = {\boldsymbol{0}}}} $$ (10)

    式中, $ {\boldsymbol{U}} $和$ {\boldsymbol{\ddot U}} $为全局位移向量及其二阶导数. 假设式(10)的解为${\boldsymbol{\varphi }}{{\text{e}}^{{\rm{a}}\omega t}}$(${\rm{a}} = \sqrt { - 1}$), 然后通过转换得到结构特征问题方程[34]并进行求解, 从而获得结构的特征值响应. 结构特征问题方程如下

    $$ \left[ {{\boldsymbol{K}}({\boldsymbol{x}}) - {\lambda _i}{\boldsymbol{M}}({\boldsymbol{y}})} \right]{{{{\boldsymbol{\varphi}} }}_i} = {\boldsymbol{0}} $$ (11)

    式中, ${\boldsymbol{x}}$, ${\boldsymbol{y}}$为成像参数矩阵, ${\lambda _i}$为第$i$阶特征值, 其中${\lambda _i} = {\omega _i}^2$, ${\omega _i}$为第$i$阶固有圆频率, ${\omega _i} = 2{\text{π }}{F_i}$, ${F_i}$为第$i$阶特征频率, ${{\boldsymbol{\varphi }}_i}$为第$i$阶特征向量.

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

    $$ \left.\begin{split} & {\text{find [}}{\boldsymbol{x}},{\boldsymbol{y}}{\text{] = [}}{x_1}{\text{ }}{x_2}{\text{ }}\cdots{\text{ }}{x_n},{y_1}{\text{ }}{y_2}{\text{ }}. ..{\text{ }}{y_n}{\text{]}} \\ & {\text{min }}f({\boldsymbol{x}},{\boldsymbol{y}}) = \sum\limits_{i = 1}^m {{w_i}} {\left( {\frac{{{\lambda _i}({\boldsymbol{x}},{\boldsymbol{y}}) - {\lambda _i}^*}}{{{\lambda _i}^*}}} \right)^2} \\ & {\text{s}}{\text{.t}}{\text{. }}{\boldsymbol{M}}\left( {\boldsymbol{y}} \right){\boldsymbol{\ddot U}}{\text{ + }}{\boldsymbol{K}}\left( {\boldsymbol{x}} \right){\boldsymbol{U}} = {{{\boldsymbol{0}} }} \\ & {{ 0 < }}{x_{\min }} \leqslant {x_j} \leqslant {x_{\max }} \\ & {{ 0 < }}{y_{{\text{min}}}} \leqslant {y_j} \leqslant {y_{\max }}{\text{ }} \\ & {\text{for }}j = {\text{1, 2, }}\cdots{\text{, }}n \end{split} \right\} $$ (12)

    式中, ${\boldsymbol{x}}$, $ {\boldsymbol{y}} $分别为元素取值在$\left[ {{x_{\min }},{x_{\max }}} \right]$, $\left[ {{y_{{\text{min}}}},{y_{{\text{max}}}}} \right]$之间的成像参数矩阵, ${x_{\min }}$和$ {x_{\max }} $分别为成像参数$ {x_j} $可取的最小值和最大值, ${y_{{\text{min}}}}$和$ {y_{{\text{max}}}} $分别为成像参数$ {y_j} $可取的最小值和最大值, 其中${x_{\min }}$, $ {x_{\max }} $, ${y_{{\text{min}}}}$和$ {y_{{\text{max}}}} $的取值由式(1)可知, $ n $为成像参数的数量, $ f({\boldsymbol{x}},{\boldsymbol{y}}) $表示弹性成像目标函数, $m$为模态阶数, $ {w_i} $为第$i$阶特征值的权重, 权重总和为1, $ {\lambda _i} $为数字模型的特征值响应, $ {\lambda _i}^* $为真实结构模型的特征值响应. 特征值权重是可调参数, 通过改变特征值权重可灵活调整各阶特征值对损伤识别的敏感性.

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

    对目标函数关于单元成像参数$ {x_j} $, ${y_j}$分别求导可得

    $$ \frac{{\partial f({\boldsymbol{x}},{\boldsymbol{y}})}}{{\partial {x_j}}} = 2\sum\limits_{i = 1}^m {{w_i} {\frac{{\partial {\lambda _i}({\boldsymbol{x}},{\boldsymbol{y}})}}{{\partial {x_j}}}} } {\frac{{{\lambda _i}({\boldsymbol{x}},{\boldsymbol{y}}) - {\lambda _i}^*}}{{{\lambda _i}^{{*^2}}}}} $$ (13)
    $$ \frac{{\partial f({\boldsymbol{x}},{\boldsymbol{y}})}}{{\partial {y_j}}} = 2\sum\limits_{i = 1}^m {{w_i} {\frac{{\partial {\lambda _i}({\boldsymbol{x}},{\boldsymbol{y}})}}{{\partial {y_j}}}} } {\frac{{{\lambda _i}({\boldsymbol{x}},{\boldsymbol{y}}) - {\lambda _i}^*}}{{{\lambda _i}^{{*^2}}}}} $$ (14)

    对结构特征问题方程(11)关于单元成像参数$ {x_j} $和${y_j}$分别求导可得

    $$ \begin{split} & \left( {\frac{{\partial {\boldsymbol{K}}({\boldsymbol{x}})}}{{\partial {x_j}}} - \frac{{\partial {\lambda _i}({\boldsymbol{x}},{\boldsymbol{y}})}}{{\partial {x_j}}}{\boldsymbol{M}}({\boldsymbol{y}}) - {\lambda _i}\frac{{\partial {\boldsymbol{M}}({\boldsymbol{y}})}}{{\partial {x_j}}}} \right){{\boldsymbol{\varphi }}_i} + \\ &\qquad \left( {{\boldsymbol{K}}({\boldsymbol{x}}) - {\lambda _i}{\boldsymbol{M}}({\boldsymbol{y}})} \right)\frac{{\partial {{\boldsymbol{\varphi }}_i}}}{{\partial {x_j}}} = {\boldsymbol{0}} \end{split} $$ (15)
    $$ \begin{split} & \left( {\frac{{\partial {\boldsymbol{K}}({\boldsymbol{x}})}}{{\partial {y_j}}} - \frac{{\partial {\lambda _i}({\boldsymbol{x}},{\boldsymbol{y}})}}{{\partial {y_j}}}{\boldsymbol{M}}({\boldsymbol{y}}) - {\lambda _i}\frac{{\partial {\boldsymbol{M}}({\boldsymbol{y}})}}{{\partial {y_j}}}} \right){{\boldsymbol{\varphi }}_i} + \\ &\qquad \left( {{\boldsymbol{K}}({\boldsymbol{x}}) - {\lambda _i}{\boldsymbol{M}}({\boldsymbol{y}})} \right)\frac{{\partial {{\boldsymbol{\varphi }}_i}}}{{\partial {y_j}}} = {\boldsymbol{0 }}\end{split} $$ (16)

    其中, 从特征问题方程可知$\left( {{\boldsymbol{K}}({\boldsymbol{x}}) - {\lambda _i}{\boldsymbol{M}}({\boldsymbol{y}})} \right)\dfrac{{\partial {{\boldsymbol{\varphi }}_i}}}{{\partial {x_j}}}$项和$\left( {{\boldsymbol{K}}({\boldsymbol{x}}) - {\lambda _i}{\boldsymbol{M}}({\boldsymbol{y}})} \right)\dfrac{{\partial {{\boldsymbol{\varphi }}_i}}}{{\partial {y_j}}}$项均为0. 使用正则化后的特征向量进行质量矩阵归一化可得

    $$ {{\boldsymbol{\varphi }}_i}^{\text{T}}{\boldsymbol{M}}({\boldsymbol{y}}){{\boldsymbol{\varphi }}_i} = 1 $$ (17)

    式中, $ {{\boldsymbol{\varphi }}_i} $是正则化后的第$i$阶特征向量, $ {\text{T}} $代表矩阵转置.

    将式(15)和式(16)两边同乘以${{{{\boldsymbol{\varphi}} _i^{\text{T}}}}}$可得特征值关于单元成像参数$ {x_j} $和${y_j}$的导数分别为

    $$ \frac{{\partial {\lambda _i}({\boldsymbol{x}},{\boldsymbol{y}})}}{{\partial {x_j}}} = {{\boldsymbol{\varphi }}_i^{\text{T}}}\left( {\frac{{\partial {\boldsymbol{K}}({\boldsymbol{x}})}}{{\partial {x_j}}} - {\lambda _i}\frac{{\partial {\boldsymbol{M}}({\boldsymbol{y}})}}{{\partial {x_j}}}} \right){{\boldsymbol{\varphi }}_i} $$ (18)
    $$ \frac{{\partial {\lambda _i}({\boldsymbol{x}},{\boldsymbol{y}})}}{{\partial {y_j}}} = {{\boldsymbol{\varphi }}_i^{\text{T}}}\left( {\frac{{\partial {\boldsymbol{K}}({\boldsymbol{x}})}}{{\partial {y_j}}} - {\lambda _i}\frac{{\partial {\boldsymbol{M}}({\boldsymbol{y}})}}{{\partial {y_j}}}} \right){{\boldsymbol{\varphi }}_i} $$ (19)

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

    $$ \frac{{\partial f({\boldsymbol{x}})}}{{\partial {x_j}}} = 2\sum\limits_{i = 1}^m {{w_i} {\frac{{\partial {\lambda _i}({\boldsymbol{x}})}}{{\partial {x_j}}}} } {\frac{{{\lambda _i}({\boldsymbol{x}}) - {\lambda _i}^*}}{{{\lambda _i}^{{*^2}}}}} $$ (20)

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

    $$ \frac{{\partial {\lambda _i}({\boldsymbol{x}})}}{{\partial {x_j}}} = {{\boldsymbol{\varphi }}_i^{\text{T}}}\left( {\frac{{\partial {\boldsymbol{K}}({\boldsymbol{x}})}}{{\partial {x_j}}} - {\lambda _i}\frac{{\partial {\boldsymbol{M}}({\boldsymbol{x}})}}{{\partial {x_j}}}} \right){{\boldsymbol{\varphi }}_i} $$ (21)

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

    $$ \begin{split} & \frac{{\partial f({\boldsymbol{x}})}}{{\partial {x_j}}} = 2\sum\limits_{i = 1}^m {\left[ {{w_i}{{\boldsymbol{\varphi }}_i^{\text{T}}}\left( {\frac{{\partial {\boldsymbol{K}}({\boldsymbol{x}})}}{{\partial {x_j}}} - {\lambda _i}\frac{{\partial {\boldsymbol{M}}({\boldsymbol{x}})}}{{\partial {x_j}}}} \right){{\boldsymbol{\varphi }}_i} \cdot } \right.} \\ &\qquad \left. { {\frac{{{\lambda _i}({\boldsymbol{x}}) - {\lambda _i}^*}}{{{\lambda _i}^{{*^2}}}}} } \right] \end{split} $$ (22)

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

    $$ \frac{{\partial f({\boldsymbol{y}})}}{{\partial {y_j}}} = 2\sum\limits_{i = 1}^m {{w_i}{\frac{{\partial {\lambda _i}({\boldsymbol{y}})}}{{\partial {y_j}}}} } {\frac{{{\lambda _i}({\boldsymbol{y}}) - {\lambda _i}^*}}{{{\lambda _i}^{{*^2}}}}} $$ (23)

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

    $$ \frac{{\partial {\lambda _i}({\boldsymbol{y}})}}{{\partial {y_j}}} = {{\boldsymbol{\varphi }}_i^{\text{T}}}\left( {\frac{{\partial {\boldsymbol{K}}}}{{\partial {y_j}}} - {\lambda _i}\frac{{\partial {\boldsymbol{M}}({\boldsymbol{y}})}}{{\partial {y_j}}}} \right){{\boldsymbol{\varphi }}_i} $$ (24)

    由于${\boldsymbol{K}}$为常数, 因此可得

    $$ \frac{{\partial {\boldsymbol{K}}}}{{\partial {y_j}}} = {\boldsymbol{0}} $$ (25)

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

    $$ \frac{{\partial {\lambda _i}({\boldsymbol{y}})}}{{\partial {y_j}}} = {{\boldsymbol{\varphi }}_i^{\text{T}}}\left( { - {\lambda _i}\frac{{\partial {\boldsymbol{M}}({\boldsymbol{y}})}}{{\partial {y_j}}}} \right){{\boldsymbol{\varphi }}_i} $$ (26)

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

    $$ \frac{{\partial f({\boldsymbol{y}})}}{{\partial {y_j}}} = 2\sum\limits_{i = 1}^m {{w_i}{{\boldsymbol{\varphi }}_i^{\text{T}}}\left( { - {\lambda _i}\frac{{\partial {\boldsymbol{M}}({\boldsymbol{y}})}}{{\partial {y_j}}}} \right){{\boldsymbol{\varphi }}_i}} {\frac{{{\lambda _i}({\boldsymbol{y}}) - {\lambda _i}^*}}{{{\lambda _i}^{{*^2}}}}} $$ (27)

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

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

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

    $$ {{\boldsymbol{x}}^{r + 1}} = \max \left({x_{\min }},\min \left({x_{\max }},{{\boldsymbol{x}}^r} - \alpha \frac{{\partial f({\boldsymbol{x}},{\boldsymbol{y}})}}{{\partial {\boldsymbol{x}}}}\right)\right) $$ (28)
    $$ {{\boldsymbol{y}}^{r + 1}} = \max \left({y_{\min }},\min \left({y_{\max }},{{\boldsymbol{y}}^r} - \alpha \frac{{\partial f({\boldsymbol{x}},{\boldsymbol{y}})}}{{\partial {\boldsymbol{y}}}}\right)\right) $$ (29)

    式中, $r$为迭代步数, $ \alpha $为迭代步长, ${\text{min}}$与${\text{max}}$分别代表取括号中最小值与最大值.

    为了使优化模型的收敛速度保持恒定, 分别将$\dfrac{{\partial f({\boldsymbol{x}},{\boldsymbol{y}})}}{{\partial {\boldsymbol{x}}}}$和$\dfrac{{\partial f({\boldsymbol{x}},{\boldsymbol{y}})}}{{\partial {\boldsymbol{y}}}}$进行归一化处理

    $$ {{\boldsymbol{S}}_1} = \dfrac{{\dfrac{{\partial f({\boldsymbol{x}},{\boldsymbol{y}})}}{{\partial {\boldsymbol{x}}}}}}{{\max \left( {\dfrac{{\partial f({\boldsymbol{x}},{\boldsymbol{y}})}}{{\partial {\boldsymbol{x}}}}} \right) - \min \left( {\dfrac{{\partial f({\boldsymbol{x}},{\boldsymbol{y}})}}{{\partial {\boldsymbol{x}}}}} \right)}} $$ (30)
    $$ {{\boldsymbol{S}}_2} = \dfrac{{\dfrac{{\partial f({\boldsymbol{x}},{\boldsymbol{y}})}}{{\partial {\boldsymbol{y}}}}}}{{\max \left( {\dfrac{{\partial f({\boldsymbol{x}},{\boldsymbol{y}})}}{{\partial {\boldsymbol{y}}}}} \right) - \min \left( {\dfrac{{\partial f({\boldsymbol{x}},{\boldsymbol{y}})}}{{\partial {\boldsymbol{y}}}}} \right)}} $$ (31)

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

    $$ {{\boldsymbol{h}}_1}^{r + 1} = \frac{{{{\boldsymbol{S}}_1}^{r + 1} + {{\boldsymbol{S}}_1}^r}}{2} $$ (32)
    $$ {{\boldsymbol{h}}_2}^{r + 1} = \frac{{{{\boldsymbol{S}}_2}^{r + 1} + {{\boldsymbol{S}}_2}^r}}{2} $$ (33)

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

    $$ {{\boldsymbol{x}}^{r + 1}} = \max ({x_{\min }},\min ({x_{\max }},{{\boldsymbol{x}}^r} - \alpha {{\boldsymbol{h}}_1}^{r + 1})) $$ (34)
    $$ {{\boldsymbol{y}}^{r + 1}} = \max ({y_{\min }},\min ({y_{\max }},{{\boldsymbol{y}}^r} - \alpha {{\boldsymbol{h}}_2}^{r + 1})) $$ (35)

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

    $$ {{\boldsymbol{x}}^{r + 1}} = \max ({x_{\min }},\min ({x_{\max }},{{\boldsymbol{x}}^r} - \alpha {{\boldsymbol{h}}_1}^{r + 1})) $$ (36)

    其中

    $$\qquad\quad {{\boldsymbol{h}}_1}^{r + 1} = \frac{{{{\boldsymbol{S}}_1}^{r + 1} + {{\boldsymbol{S}}_1}^r}}{2} $$ (37)
    $$\qquad\quad {{\boldsymbol{S}}_1} = \dfrac{{\dfrac{{\partial f({\boldsymbol{x}})}}{{\partial {\boldsymbol{x}}}}}}{{\max \left( {\dfrac{{\partial f({\boldsymbol{x}})}}{{\partial {\boldsymbol{x}}}}} \right) - \min \left( {\dfrac{{\partial f({\boldsymbol{x}})}}{{\partial {\boldsymbol{x}}}}} \right)}} $$ (38)

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

    $$ {{\boldsymbol{y}}^{r + 1}} = \max ({y_{\min }},\min ({y_{\max }},{{\boldsymbol{y}}^r} - \alpha {{\boldsymbol{h}}_2}^{r + 1})) $$ (39)

    其中

    $$ {{\boldsymbol{h}}_2}^{r + 1} = \frac{{{{\boldsymbol{S}}_2}^{r + 1} + {{\boldsymbol{S}}_2}^r}}{2} $$ (40)
    $$ {{\boldsymbol{S}}_2} = \dfrac{{\dfrac{{\partial f({\boldsymbol{y}})}}{{\partial {\boldsymbol{y}}}}}}{{\max \left( {\dfrac{{\partial f({\boldsymbol{y}})}}{{\partial {\boldsymbol{y}}}}} \right) - \min \left( {\dfrac{{\partial f({\boldsymbol{y}})}}{{\partial {\boldsymbol{y}}}}} \right)}} $$ (41)

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

    $$ \left| {{f^{r + 1}} - {f^r}} \right| \leqslant \varepsilon {\text{ or }}r \geqslant {r_{\max }} $$ (42)

    式中, $\varepsilon $代表一个极小正数, ${r_{\max }}$为最大迭代步数.

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

    为验证本方法识别结构中第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梁结构
    Figure  3.  Two-dimensional half-MBB beam structure with built-in defect

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

    图  4  二维半MBB梁结构弹性成像迭代曲线
    Figure  4.  Iterative curve of elastography of two-dimensional half-MBB beam structure
    表  1  算例1中损伤模型和无损模型的特征频率(Hz)
    Table  1.  Eigen-frequencies of damaged model and undamaged model in the example 1 (Hz)
    Eigen-frequencyVibration order
    123456
    damaged model54.34197.25235.42371.65537.65616.35
    undamaged model54.81194.30235.36368.93548.82621.28
    下载: 导出CSV 
    | 显示表格
    图  5  二维半MBB梁结构成像结果
    Figure  5.  Two-dimensional half-MBB beam structure imaging results

    为验证本方法识别结构中第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的结冰区域, 结冰区域的弹性模量和材料密度分别设置为${E^0}$, ${\rho ^0}$, 在结构右上角位置节点上施加数值大小为结构域质量20%的非结构集中质量. 本算例考虑与算例5.1相同的特征值阶次和特征值权重.

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

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

    图  7  具有多孔特征的二维悬臂梁结构弹性成像迭代曲线
    Figure  7.  Iterative curve of elastography of two-dimensional cantilever beam structure with porous features
    表  2  算例2中损伤模型和无损模型的特征频率(Hz)
    Table  2.  Eigen-frequencies of damaged model and undamaged model in the example 2 (Hz)
    Eigen-frequencyVibration order
    123456
    damaged model35.84236.80411.77478.93883.40988.36
    undamaged model35.88236.06410.98479.23883.72987.73
    下载: 导出CSV 
    | 显示表格
    图  8  具有多孔特征的二维悬臂梁结构成像结果
    Figure  8.  Imaging results of two-dimensional cantilever beam structure with porous features

    为进一步验证本方法识别结构第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的积水区域, 由于积水仅影响结构质量, 因此积水区域的材料弹性模量不变, 材料密度设置为${\rho ^0}$. 本算例施加非结构集中质量的位置和数值大小与算例5.1保持一致, 并且考虑与算例5.1相同的特征值阶次和特征值权重.

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

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

    图  10  具有多孔特征的二维半MBB梁结构弹性成像迭代曲线
    Figure  10.  Iterative curve of elastography of two-dimensional half-MBB beam structure with porous features
    表  3  算例3中损伤模型和无损模型的特征频率(Hz)
    Table  3.  Eigen-frequencies of damaged model and undamaged model in the example 3 (Hz)
    Eigen-frequencyVibration order
    123456
    damaged model71.09202.34446.90633.191125.671209.26
    undamaged model71.24202.35449.28634.961122.411216.35
    下载: 导出CSV 
    | 显示表格
    图  11  具有多孔特征的二维半MBB梁结构成像结果
    Figure  11.  Imaging results of two-dimensional half-MBB beam structure with porous features

    为进一步验证本方法的通用性, 开展三维结构弹性成像的研究. 如图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  内置损伤的三维悬臂梁结构
    Figure  12.  Three-dimensional cantilever beam structure with built-in damaged

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

    图  13  内置缺损的三维悬臂梁结构弹性成像迭代曲线
    Figure  13.  Iterative curve of elastography of three-dimensional cantilever beam structure with built-in defect
    图  14  内置缺损的三维悬臂梁结构成像结果
    Figure  14.  Imaging results of three-dimensional cantilever beam structure with built-in defect
    图  15  内置结冰的三维悬臂梁结构弹性成像迭代曲线
    Figure  15.  Iterative curve of elastography of three-dimensional cantilever beam structure with built-in ice
    图  16  内置结冰的三维悬臂梁结构成像结果
    Figure  16.  Imaging results of three-dimensional cantilever beam structure with built-in ice
    图  17  内置积水的三维悬臂梁结构弹性成像迭代曲线
    Figure  17.  Iterative curve of elastography of three-dimensional cantilever beam structure with built-in water accumulation
    图  18  内置积水的三维悬臂梁结构成像结果
    Figure  18.  Imaging results of three-dimensional cantilever beam structure with built-in water accumulation

    本文提出了一种考虑刚度和质量耦合效应的结构弹性成像方法, 构建了缺损、夹杂、结冰和积水多种损伤类型的统一表征, 并构建了损伤表征、力学模型和特征值响应之间的映射关系, 建立并求解了结构弹性成像模型. 研究结论如下:

    (1)本方法可有效识别结构缺损、结冰、积水多种损伤类型的位置、数量和形状;

    (2)结构弹性成像结果不受特定边界条件限制, 无需先验信息可自动反演单/多个损伤;

    (3)本方法在二维和三维结构弹性成像问题中具有通用性. 本方法丰富了现有机械结构损伤识别方法的理论内涵, 有助于实现重大机械工程装备的服役安全评价.

  • 图  1   多种损伤类型表征

    Figure  1.   Characterization of damages

    图  2   考虑刚度和质量耦合效应的结构弹性成像流程

    Figure  2.   Structural elastography process considering the coupling effect of stiffness and mass

    图  3   内置缺损的二维半MBB梁结构

    Figure  3.   Two-dimensional half-MBB beam structure with built-in defect

    图  4   二维半MBB梁结构弹性成像迭代曲线

    Figure  4.   Iterative curve of elastography of two-dimensional half-MBB beam structure

    图  5   二维半MBB梁结构成像结果

    Figure  5.   Two-dimensional half-MBB beam structure imaging results

    图  6   内置结冰的二维悬臂梁结构

    Figure  6.   Two-dimensional cantilever beam structure with built-in ice

    图  7   具有多孔特征的二维悬臂梁结构弹性成像迭代曲线

    Figure  7.   Iterative curve of elastography of two-dimensional cantilever beam structure with porous features

    图  8   具有多孔特征的二维悬臂梁结构成像结果

    Figure  8.   Imaging results of two-dimensional cantilever beam structure with porous features

    图  9   内置积水的二维半MBB梁结构

    Figure  9.   Two-dimensional half-MBB beam structure with built-in water accumulation

    图  10   具有多孔特征的二维半MBB梁结构弹性成像迭代曲线

    Figure  10.   Iterative curve of elastography of two-dimensional half-MBB beam structure with porous features

    图  11   具有多孔特征的二维半MBB梁结构成像结果

    Figure  11.   Imaging results of two-dimensional half-MBB beam structure with porous features

    图  12   内置损伤的三维悬臂梁结构

    Figure  12.   Three-dimensional cantilever beam structure with built-in damaged

    图  13   内置缺损的三维悬臂梁结构弹性成像迭代曲线

    Figure  13.   Iterative curve of elastography of three-dimensional cantilever beam structure with built-in defect

    图  14   内置缺损的三维悬臂梁结构成像结果

    Figure  14.   Imaging results of three-dimensional cantilever beam structure with built-in defect

    图  15   内置结冰的三维悬臂梁结构弹性成像迭代曲线

    Figure  15.   Iterative curve of elastography of three-dimensional cantilever beam structure with built-in ice

    图  16   内置结冰的三维悬臂梁结构成像结果

    Figure  16.   Imaging results of three-dimensional cantilever beam structure with built-in ice

    图  17   内置积水的三维悬臂梁结构弹性成像迭代曲线

    Figure  17.   Iterative curve of elastography of three-dimensional cantilever beam structure with built-in water accumulation

    图  18   内置积水的三维悬臂梁结构成像结果

    Figure  18.   Imaging results of three-dimensional cantilever beam structure with built-in water accumulation

    表  1   算例1中损伤模型和无损模型的特征频率(Hz)

    Table  1   Eigen-frequencies of damaged model and undamaged model in the example 1 (Hz)

    Eigen-frequencyVibration order
    123456
    damaged model54.34197.25235.42371.65537.65616.35
    undamaged model54.81194.30235.36368.93548.82621.28
    下载: 导出CSV

    表  2   算例2中损伤模型和无损模型的特征频率(Hz)

    Table  2   Eigen-frequencies of damaged model and undamaged model in the example 2 (Hz)

    Eigen-frequencyVibration order
    123456
    damaged model35.84236.80411.77478.93883.40988.36
    undamaged model35.88236.06410.98479.23883.72987.73
    下载: 导出CSV

    表  3   算例3中损伤模型和无损模型的特征频率(Hz)

    Table  3   Eigen-frequencies of damaged model and undamaged model in the example 3 (Hz)

    Eigen-frequencyVibration order
    123456
    damaged model71.09202.34446.90633.191125.671209.26
    undamaged model71.24202.35449.28634.961122.411216.35
    下载: 导出CSV
  • [1]

    Ophir J, Céspedes I, Ponnekanti H, et al. Elastography: A quantitative method for imaging the elasticity of biological tissues. Ultrasonic Imaging, 1991, 13(2): 111-134 doi: 10.1177/016173469101300201

    [2]

    Karalilova R, Doykova K, Batalov Z, et al. Spleen elastography in patients with systemic sclerosis. Rheumatology International, 2021, 41(3): 633-641 doi: 10.1007/s00296-020-04772-5

    [3] 陈金龙, 潘航, 林祥龙等. 原位测量软材料弹性模量的光学相干弹性成像方法研究. 天津大学学报(自然科学与工程技术版), 2023, 56(2): 111-118 (Chen Jinlong, Pan Hang, Ling Xianglong, et al. Optical coherence elastography for the in situ measurement of the elastic modulus of soft materials. Journal of Tianjin University (Science and Technology), 2023, 56(2): 111-118 (in Chinese)

    Chen Jinlong, Pan Hang, Ling Xianglong, et al. Optical coherence elastography for the in situ measurement of the elastic modulus of soft materials. Journal of Tianjin University(Science and Technology), 2023, 56(2): 111-118(in Chinese))

    [4] 张鹏辉, 赵扬, 李鹏等. 超声成像检测技术研究进展综述. 激光与光电子学进展, 2022, 59(2): 1-15 (Zhang Penghui, Zhao Yang, Li Peng, et al. Research progress in ultrasonic imaging detection technology. Laser &Optoelectronics Progress, 2022, 59(2): 1-15 (in Chinese)

    Zhang Penghui, Zhao Yang, Li Peng, et al. Research progress in ultrasonic imaging detection technology. Laser & Optoelectronics Progress, 2022, 59(2): 1-15(in Chinese))

    [5] 傅健, 张昌盛, 朱国港等. X射线分层层析成像技术及在航空航天领域的应用. 航空制造技术, 2019, 62(14): 49-54 (Fu Jian, Zhang Changsheng, Zhu Guogang, et al. Development and application of x-ray computed laminography for aerospace. Aeronautical Manufacturing Technology, 2019, 62(14): 49-54 (in Chinese)

    Fu jian, Zhang Changsheng, Zhu Guogang, et al. Development and application of x-ray computed laminography for aerospace. Aeronautical Manufacturing Technology, 2019, 62(14): 49-54(in Chinese))

    [6] 翁顺, 朱宏平. 基于有限元模型修正的土木结构损伤识别方法. 工程力学, 2021, 38(3): 1-16 (Weng Shun, Zhu Hongping. Damage identification of civil structures based on finite element model updating. Engineering Mechanics, 2021, 38(3): 1-16 (in Chinese) doi: 10.6052/j.issn.1000-4750.2020.06.ST02

    Weng Shun, Zhu Hongping. Damage identification of civil structures based on finite element model updating. Engineering Mechanics, 2021, 38(3): 1-16(in Chinese)) doi: 10.6052/j.issn.1000-4750.2020.06.ST02

    [7]

    Minh H, Sang-To T, Wahab MA, et al. Structural damage identification in thin-shell structures using a new technique combining finite element model updating and improved cuckoo search algorithm. Advances in Engineering Software, 2022, 173: 103206 doi: 10.1016/j.advengsoft.2022.103206

    [8]

    Singh T, Sehgal S, Prakash C, et al. Real-time structural health monitoring and damage identification using frequency response functions along with finite element model updating technique. Sensors, 2022, 22(12): 4546 doi: 10.3390/s22124546

    [9]

    Simoen E, De Roeck G, Lombaert G. Dealing with uncertainty in model updating for damage assessment: a review. Mechanical Systems and Signal Processing, 2015, 56-57: 123-149 doi: 10.1016/j.ymssp.2014.11.001

    [10]

    Sigmund O, Maute K. Topology optimization approaches: a comparative review. Structural and Multidisciplinary Optimization, 2013, 48(6): 1031-1055 doi: 10.1007/s00158-013-0978-6

    [11]

    Zhang WH, Zhang ZD, Zhu JH, et al. Structural topology optimization: extensibility and attainability. Science China: Technological Sciences, 2014, 57(7): 1310-1321 doi: 10.1007/s11431-014-5580-7

    [12]

    Ibhadode O, Zhang ZD, Sixt J, et al. Topology optimization for metal additive manufacturing: current trends, challenges, and future outlook. Virtual and Physical Prototyping, 2023, 18(1): e2181192 doi: 10.1080/17452759.2023.2181192

    [13]

    Wang MY, Wang XM, Guo DM. A level set method for structural topology optimization. Computer Methods in Applied Mechanics and Engineering, 2003, 192(1-2): 227-246 doi: 10.1016/S0045-7825(02)00559-5

    [14]

    Wang YG, Yang HD, Kang Z. Velocity field level set method incorporating topological derivatives for topology optimization. Journal of Applied Mechanics, 2022, 89(6): 061002 doi: 10.1115/1.4053989

    [15]

    Guo X, Zhang WS, Zhong WL. Doing topology optimization explicitly and geometrically—a new moving morphable components based framework. Journal of Applied Mechanics, 2014, 81(8): 081009 doi: 10.1115/1.4027609

    [16] 李佳霖, 赵剑, 孙直等. 基于移动可变形组件法(MMC)的运载火箭传力机架结构的轻量化设计. 力学学报, 2022, 54(1): 244-251 (Li Jialin, Zhao Jian, Sun Zhi, et al. Lightweight design of transmission frame structures for launch vehicles based on moving morphable components (MMC) approach. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(1): 244-251 (in Chinese)

    Li Jialin, Zhao jian, Sun Zhi, et al. Lightweight design of transmission frame structures for launch vehicles based on moving morphable components (MMC) approach. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(1): 244-251(in Chinese))

    [17]

    Sigmund O. A 99 line topology optimization code written in matlab. Structural and Multidisciplinary Optimization, 2001, 21(2): 120-127 doi: 10.1007/s001580050176

    [18]

    Zhang LX, Yang G, Hu D, et al. An approach based on level set method for void identification of continuum structure with time-domain dynamic response. Applied Mathematical Modelling, 2019, 75(3): 446-480

    [19]

    Mei Y, Du ZL, Zhao DM, et al. Moving morphable inclusion approach: an explicit framework to solve inverse problem in elasticity. Journal of Applied Mechanics, 2021, 88(4): 041001 doi: 10.1115/1.4049142

    [20]

    Damghani F, Tavakkoli SM. Structural damage detection in plane stress problems by using time domain responses and topology optimization. International Journal of Optimization in Civil Engineering, 2023, 13(2): 257-274

    [21]

    Zhang LX, Yang G, Hu D. Identification of voids in structures based on level set method and fem. International Journal of Computational Methods, 2018, 15(1): 1850015

    [22]

    Zhang WS, Du ZL, Sun G, et al. A level set approach for damage identification of continuum structures based on dynamic responses. Journal of Sound and Vibration, 2017, 386: 100-115 doi: 10.1016/j.jsv.2016.06.014

    [23] 江守燕, 赵林鑫, 杜成斌. 基于频率和模态保证准则的结构内部多缺陷反演. 力学学报, 2019, 51(4): 1091-1100 (Jiang Shouyan, Zhao Linxin, Du Chengbin. Identification of multiple flaws in structures based on frequency and modal assurance criteria. Chinese Journal of Theoretical and Applied Mechanics, 2019, 51(4): 1091-1100 (in Chinese)

    Jiang Shouyan, Zhao Linxin, Du Chengbin. Identification of multiple flaws in structures based on frequency and modal assurance criteria. Chinese Journal of Theoretical and Applied Mechanics, 2019, 51(4): 1091-1100(in Chinese))

    [24]

    Lee JS, Kim JE, Kim YY. Damage detection by the topology design formulation using modal parameters. International Journal for Numerical Methods in Engineering, 2007, 69(7): 1480-1498 doi: 10.1002/nme.1817

    [25]

    Abdollahi F, Tavakkoli SM. Damage identification by using modal expansion and topology optimization in three dimensional elasticity problems. International Journal of Optimization in Civil Engineering, 2019, 9(4): 543-560

    [26]

    Saito A, Sugai R, Wang ZX, et al. Damage identification using noisy frequency response functions based on topology optimization. Journal of Sound and Vibration, 2023, 545: 117412 doi: 10.1016/j.jsv.2022.117412

    [27] 付君健, 李帅虎, 李好等. 基于拓扑优化的结构弹性成像方法. 力学学报, 2022, 54(5): 1331-1340 (Fu Junjian, Li Shuaihu, Li Hao, et al. Structural elastography method based on topology optimization. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(5): 1331-1340 (in Chinese)

    Fu Junjian, Li Shuaihu, Li Hao, et al. Structural elastography method based on topology optimization. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(5): 1331-1340(in Chinese))

    [28]

    Kim G, Sterkenburg R, Tsutsui W. Investigating the effects of fluid intrusion on nomex® honeycomb sandwich structures with carbon fiber facesheets. Composite Structures, 2018, 206: 535-549 doi: 10.1016/j.compstruct.2018.08.054

    [29]

    Song ZH, Luong S, Whisler D, et al. Honeycomb core failure mechanism of CFRP/nomex sandwich panel under multi-angle impact of hail ice. International Journal of Impact Engineering, 2021, 150: 103817 doi: 10.1016/j.ijimpeng.2021.103817

    [30]

    Abrate S. Soft impacts on aerospace structures. Progress in Aerospace Sciences, 2016, 81: 1-17 doi: 10.1016/j.paerosci.2015.11.005

    [31]

    Siddique SH, Hazell PJ, Wang HX, et al. Lessons from nature: 3 d printed bio-inspired porous structures for impact energy absorption—a review. Additive Manufacturing, 2022, 58: 103051 doi: 10.1016/j.addma.2022.103051

    [32]

    Fu JJ, Sun PF, Du YX, et al. Design and mechanical characterization of an s-based TPMS hollow isotropic cellular structure. Computer Modeling in Engineering & Sciences, 2022, 131(2): 695-713

    [33]

    Bendsøe MP, Sigmund O. Material interpolation schemes in topology optimization. Archive of Applied Mechanics, 1999, 69(9): 635-654

    [34] 胡骏, 亢战. 考虑可控性的压电作动器拓扑优化设计. 力学学报, 2019, 51(4): 1073-1081 (Hu Jun, Kang Zhan. Topology optimization of piezoelectric actuator considering controllability. Chinese Journal of Theoretical and Applied Mechanics, 2019, 51(4): 1073-1081 (in Chinese)

    Hu Jun, Kang Zhan. Topology optimization of piezoelectric actuator considering controllability. Chinese Journal of Theoretical and Applied Mechanics, 2019, 51(4): 1073-1081(in Chinese))

    [35]

    Fan W, Qiao PZ. Vibration-based damage identification methods: a review and comparative study. Structural Health Monitoring, 2011, 10(1): 83-111 doi: 10.1177/1475921710365419

    [36] 付君健, 张跃, 杜义贤等. 周期性多孔结构特征值拓扑优化. 振动与冲击, 2022, 41(3): 73-81 (Fu Junjian, Zhang Yue, Du Yixian, et al. Eigenvalue topology optimization of periodic cellular structures. Journal of Vibration and Shock, 2022, 41(3): 73-81 (in Chinese)

    Fu Junjian, Zhang Yue, Du Yixian, et al. Eigenvalue topology optimization of periodic cellular structures. Journal of Vibration and Shock, 2022, 41(3): 73-81(in Chinese))

图(18)  /  表(3)
计量
  • 文章访问数:  544
  • HTML全文浏览量:  179
  • PDF下载量:  87
  • 被引次数: 0
出版历程
  • 收稿日期:  2023-06-01
  • 录用日期:  2023-08-14
  • 网络出版日期:  2023-08-15
  • 刊出日期:  2023-09-12

目录

/

返回文章
返回