引 言

基于有关实际测试数据，建立有效的反演分析模型，检测识别出结构内部缺陷的位置、尺寸和类型等参数，一般包括正分析和目标函数极小化迭代过程两部分.早期，有限元法和边界元法被广泛地用于正分析数值计算，但在每次迭代优化中，有限元法都需要进行网格重剖分，在成百上千次的迭代过程中其计算工作量巨大.边界元法尽管只需对边界单元进行重划分，但其工作量也不容小觑，计算效率同样受到很大的影响.1999年，美国西北大学的学者基于单位分解理论^[1,2]提出的扩展有限元法(XFEM)^[3]弥补了这些不足，它是在标准有限元框架下，连续区域采用标准有限元法，在包含不连续边界的局部区域内，修正有限元的位移场近似函数，并增加了对不连续边界的描述方法，扩展有限元法为解决结构包含裂纹、孔洞等缺陷的数值仿真模拟提供了新的有效途径^[4,5]，扩展有限元法与水平集法^[6]的结合避免了反演分析迭代过程中网格的重剖分.

在国外已有的相关研究中，有人基于扩展有限元法和遗传算法提出了一种裂纹反演识别的数值方法，并利用该方法识别出静、动载下平板薄膜结构内部的裂纹^[7,8]；还有人运用扩展有限元法和遗传算法反演识别出线弹性结构内各种类型的缺陷，包括直裂纹、圆形孔洞及不规则形孔洞等^[9]；有学者基于扩展有限元法和遗传算法采用椭圆形孔洞近似识别任意形状的缺陷，但这些研究中采用的遗传算法效率较低且容易陷入局部最优解^[10].除了遗传算法，一些其他的智能优化算法(如梯度优化算法^[11]、人工蜂群算法^[12,13]、多级坐标搜索算法^[14]等)也被较为广泛地用于反分析问题中.其中，人工蜂群算法在每次迭代中都采用全局和局部搜索，找到最优解的概率大大增加并可很好地避免局部最优，而且收敛速度较快.人工蜂群算法与扩展有限元法的结合可有效地减少反演分析的计算工作量.在国内已有的研究中，主要应用扩展有限元法进行正分析问题的研究，如：断裂力学问题^{[15,16,17,18]}，扩展有限元法在反分析问题中的应用研究尚未见报道.

文中结合扩展有限元法和人工蜂群智能优化算法的优点，重点研究结构内部缺陷(夹杂)的反演分析模型，阐述利用扩展有限元法进行结构内部圆形或椭圆形缺陷(夹杂)的反演分析时扩展有限元位移模式的构建以及水平集函数的表征，给出基于扩展有限元法和人工蜂群智能优化算法的结构内部缺陷(夹杂)的反演分析流程.最后，通过若干算例验证建立的反演分析模型能够准确地探测结构内部的单个缺陷(夹杂).

1 反分析问题

在弹性力学正分析问题中，一般是已知结构的平衡方程、几何方程、本构关系以及边界条件等物理特性，求解结构的响应量(应变、位移等). 而在反分析问题中，往往是已知结构某些关键点的响应量(如：位移)来反演结构一些未知的物理特性. 结构内部缺陷(夹杂)的反演是一类典型的反分析问题.

图1为含缺陷/夹杂的结构体及其响应测试示意图.为准确探测结构内部存在的缺陷(孔洞、裂缝)及夹杂的位置及尺寸，需要结合一些优化算法寻找合适的参数模型估计结构关键点的响应，并与结构真实响应进行比较，使得二者之间误差最小化.

参数模型通过下列一组向量表示

$$ { \alpha} = \left\{ {\alpha _1 ,} {\alpha _2 ,} { \cdots ,} {\alpha _n } \right\} \in { R}^n $$

(1)

其中，$n$为待反演的参数个数.

反演单个圆形夹杂(孔洞)时，参数为

$$ { \alpha}_{{\rm c }} = \left\{ {x_{\rm c} ,} {y_{\rm c} ,} {r_{\rm c} } \right\} $$

(2)

其中，$\left( x_{\rm c} ,y_{\rm c} \right)$为圆心坐标，$r_{\rm c} $为圆的半径.

反演单个椭圆形夹杂(孔洞)时，参数为

$$ { \alpha}_{{\rm e }} = \left\{ {x_{\rm e} ,} {y_{\rm e} ,} {a,} {b,} \ \beta \right\} $$

(3)

其中，$\left( {x_{\rm e} ,y_{\rm e} } \right)$为椭圆中心点坐标；$a$为椭圆半长轴，$b$为椭圆半短轴；$\beta $为方位角，即整体坐标系与局部坐标系之间的夹角.

反演单个裂纹时，参数为

$$ { \alpha}_{{\rm crack}} = \left\{ {x_{{\rm cb}} ,} {y_{{\rm cb}} ,} {x_{{\rm ce}} ,} {y_{{\rm ce}} } \right\} $$

(4)

其中，$\left(x_{{\rm cb}} ,y_{{\rm cb}} \right)$为裂纹起点坐标，$\left( x_{{\rm ce}} ,y_{{\rm ce}} \right)$为裂纹终点坐标.

反演分析的目标函数为

$$ O\left( { \alpha} \right) = \sum {\dfrac{\left\| { { u }^{\rm c} \left( { \alpha} \right) - { u}^{\rm m} } \right\|}{\left\| {{ u }^{\rm m} } \right\|}} $$

(5)

式中，${ u}^{\rm c} \left( { \alpha} \right)$为结构关键点响应量的估计值，${ u}^{\rm m} $为结构关键点响应量的真实值.

因此，反分析问题可描述为结合智能算法搜索寻找最佳的参数向量

$$ \mathop{\alpha}\limits^\frown = \left\{ {\mathop{\alpha}\limits^\frown}_1 ,\ {\mathop{\alpha}\limits^\frown}_2 ,\cdots ,\ {\mathop{\alpha}\limits^\frown}_n \right\} \in { R}^n $$

(6)

且

$$ { \alpha}_{\min } \leqslant \mathop{\alpha}\limits^\frown \leqslant { \alpha}_{\max } $$

(7)

使得目标函数式(5)最小化. 其中，${ \alpha}_{\max } $和${\alpha }_{\min } $分别为参数的上限值和下限值.

2 扩展有限元法

采用扩展有限元法计算含缺陷(夹杂)结构的响应量估计值. 扩展有限元法通过引入非连续位移模式，使得不连续位移场的描述独立于网格划分，可以在不重新划分网格的前提下，通过改变水平集函数反映缺陷(夹杂)的位置及大小.

扩展有限元法的位移模式可表示为

$$ { u}^{\rm h} ({ x})=\underbrace{\sum_{i\in I}N_i({ x}){ u}_i}_{\hbox{ 常规FEM部分}}+ \underbrace{\sum_{i\in I^*}N^*_i({ x}){\psi}_i({ x}){ a}_i}_{\hbox{ 改进部分}} $$

(8)

式中，$N_i ({ x})$为结点$i$对应的常规有限元法的插值形函数；${ u }_i $为常规有限元法结点$i$处的位移未知量；$\psi _i ({x })$为扩展有限元法的改进函数；${ a}_i $为结点$i$的改进位移未知量；$N_i^\ast ({ x})$为单位分解函数，其形式可以和$N_i ({ x})$相同，也可以不同，文中取二者相同.

如图2，对于被裂纹完全切割的单元，这些单元的结点(结点子集：$I_{{\rm abs}}^\ast $)通过阶跃函数$H({ x})$进行改进^[19]

$$H(x) = \left\{ {\matrix{ { + 1,(x - {x^ * }) \cdot n > 0} \hfill \cr { - 1,(x - {x^ * }) \cdot n ＜ 0} \hfill \cr } } \right.$$

(9)

式中，${ x}^\ast $是考察点${ x}$在裂纹面上的投影点；${ n} $为裂纹面的单位外法向.

对于被裂纹部分切割的单元，这些单元的结点(结点子集：$I_{{\rm br}}^{\rm \ast } $)通过裂尖改进函数$F_j ({ x})$进行改进^[19]

$$\left. \matrix{ {F_1}(r,\theta ) = \sqrt r \sin (\theta /2) \hfill \cr {F_2}(r,\theta ) = \sqrt r \cos (\theta /2) \hfill \cr {F_3}(r,\theta ) = \sqrt r \sin \theta \sin (\theta /2) \hfill \cr {F_4}(r,\theta ) = \sqrt r \sin \theta \cos (\theta /2) \hfill \cr} \right\}$$

(10)

式中，$r$和$\theta $为裂尖局部坐标系.

对于被孔洞边界切割的单元，这些单元的结点(结点子集：$I_{{\rm void}}^\ast$)通过函数$V({ x})$进行改进，若结点位于孔洞内，$V({ x}) =0$，否则$V({ x}) = 1$.

对于被夹杂边界切割的单元，这些单元的结点(结点子集：$I_{{\rm inc}}^\ast $)通过函数$\varphi ({x})$进行改进^[20]

$$ \varphi ({ x}) = \sum_{i \in I^\ast } {N_i ({ x} )\left| {\phi _i } \right|} - \left| {\sum_{i \in I^\ast } {N_i ({ x})\phi _i } } \right| $$

(11)

式中，$N_i ({ x})$为常规有限元插值形函数；$\phi _i$为结点水平集值.

对于实际问题，结构内部不规则的孔洞(夹杂)在反演分析时可采用规则的圆形或椭圆形孔洞(夹杂)进行最佳拟合.

圆形夹杂(孔洞)的水平集函数可表示成^[20]

$$ \phi ({ x}) = \left\| {{ x} - { x}_{\rm c} } \right\| - r_{\rm c} $$

(12)

式中，${ x}_{\rm c} $为圆心点的坐标向量；$r_{\rm c} $为半径.

椭圆形夹杂(孔洞)的水平集函数可表示成^[21]

$$\phi ({\rm }x) = {{{{\left( {\left( {x - {x_{\rm {c}}}} \right)\cos \theta + \left( {y - {y_{\rm {c}}}} \right)\sin \theta } \right)}^2}} \over {{a^2}}} + {{{{\left( {\left( {y - {y_{\rm {c}}}} \right)\cos \theta - \left( {x - {x_{\rm {c}}}} \right)\sin \theta } \right)}^2}} \over {{b^2}}} - 1$$

(13)

式中，$\left( {x_{\rm c} ,y_{\rm c} } \right)$为椭圆中心点坐标；$a$为椭圆的半长轴；$b$为椭圆的半短轴；$\theta $为椭圆的方向角，以逆时针为正，$\theta \in \left[{ - {\rm { }}\pi ,{\rm { }}\pi } \right]$.

如图3所示，裂纹采用两个水平集$\psi ({ x},t)$和$\phi ^k({ x},t)(k = 1,2)$描述，图中简写为$\psi$和$\phi ^k$，$\psi ({ x},t)$的零水平集$\psi ({ x},t) = 0$为裂纹面，波前水平集$\phi ^k({x},t)$与$\psi ({ x},t)$是正交的. $\psi ({ x},t)$取为符号距离函数，可表达为^[22]

$$ \psi ({ x},t) = \left\| {{ x} - { x}^\ast } \right\|{\rm sign}[({ x} - { x}^\ast ) \cdot { n}] $$

(14)

波前水平集$\phi ^k({ x},t)$为 ^[22]

$$ \phi^k({ x},t) = ({ x} - { x}_k ) \cdot { t} $$

(15)

式中，${ x}$为考察点$P$的坐标；${ x}^\ast$为考察点$P$在裂纹面上的投影点坐标；${n} $为裂纹面的单位外法向；${\rm sign}(x)$为符号函数，$x >0$时${\rm sign}(x) = 1$，$x = 0$时${\rm sign}(x) = 0$，$x ＜0$时${\rm sign}(x) = - 1$；${ x}_k $为第$k$个裂尖坐标；${t} $为第$k$个裂纹尖端处的单位切向矢量.

3 人工蜂群算法

在扩展有限元法中扩充人工蜂群算法反演最佳的模型参数. 人工蜂群算法是一种将求解优化问题的过程转化为蜜蜂群寻找优良蜜源过程的仿生智能计算方法，在基本的人工蜂群算法中主要有3个控制参数：蜜源的数量($N_{\rm Food}$)、优化极值不变的最大搜索限定次数($N_{\rm Limit}$)和最大迭代次数($N_{\rm Iter}$)，这使得算法的空间复杂度较低.在该算法中，蜜源的位置对应优化问题的可行解(缺陷的位置、大小)，蜜源的花蜜量对应可行解的质量(适应度函数值)，寻找花蜜速度对应可行解的优化速度，最大花蜜量对应优化问题最优解.适应度函数值$f\left( { \alpha} \right)$与目标函数$O\left( { \alpha} \right)$的关系为

$$ f\left( { \alpha} \right) = 1 / O\left( { \alpha} \right) = \sum \dfrac{\left\| {{ u}^{\rm m} } \right\|}{\left\| {{ u}^{\rm c} \left( { \alpha} \right) - { u}^{\rm m} } \right\|} $$

(16)

基于人工蜂群算法随机搜索蜜源(式(2),(3)或式(4)中水平集函数的参数)，并根据蜜源情况(确定的水平集参数)运用扩展有限元法求解结构关键点响应量的估计值，与结构关键点响应量的真实值比较，得出蜜源花蜜量(式(16)的函数值)，即蜜源的``收益度''. 当根据搜索的蜜源情况得出的蜜源花蜜量小于优化误差时，即可确定缺陷(夹杂)的具体位置及大小. 人工蜂群算法的求解过程为：

(1) 初始化求解参数. 初始时刻随机生成$N$个可行解矩阵$\left[ { \alpha}^1, { \alpha}^2,\\cdots , { \alpha}^{ N} \right]^{\rm T}$，某个随机可行解参数$\alpha _i^j $为^[23]

$$ \alpha _i^j = \alpha^{j}_{i\min } + \lambda _i^j \left( \alpha^{j}_{i\max } - \alpha^{j}_{i \min } \right) $$

(17)

式中，$i$取值于$\left\{ {1,2,\cdots ,n} \right\}$，为$n$维解向量的某个分量；$j$取值于$\left\{ {1,2,\cdots,N} \right\}$；$\lambda _i^j $为$\left( {0,1} \right)$之间的随机数；$\alpha^{j}_{i \max } $和$\alpha^{j}_{i \min } $为随机可行解参数$\alpha _i^j $的上限值和下限值.利用扩展有限元法计算各随机向量的适应度函数值，将排名前一半的作为初始的采蜜蜂种群.

(2) 采蜜蜂搜索阶段. 在这一阶段，采蜜蜂通过邻域搜索得到优化解向量$\beta _i^j$，搜索公式为^[23]

$$ \beta _i^j = \alpha _i^j + \delta _i^j \left( {\alpha _i^j - \alpha _i^k } \right) $$

(18)

式中，$k$取值于$\left\{ {1,2,\cdots ,n} \right\}$，且$k \ne i$，$k$和$i$均随机生成；$\delta _i^j $为$\left( { - 1,1}\right)$之间的随机数. 根据采蜜蜂搜索到的新位置$\beta _i^j $和原位置$\alpha_i^j $利用扩展有限元法计算适应度函数值$f\left( { \alpha} \right)$，若根据新位置计算的$f\left( { \alpha}\right)$更优，则该采蜜蜂的位置更新到新位置，否则位置不变.

(3)观察蜂搜索阶段. 观察蜂根据采蜜蜂种群适应度值大小选择一个采蜜蜂，并在其邻域内同样进行新位置的搜索，并按照与采蜜蜂相同的规则更新其位置. 蜜源由观察蜂选择的概率计算公式为^[23]

$$ P_i = f\left( { \alpha}^i \right) \Bigg/ \sum_{k = 1}^{ NF} f\left( { \alpha}^k \right) $$

(19)

(4)侦察蜂搜索阶段. 在某只采蜜蜂的位置周围搜索次数达到一定限值($N_{\rm Limit}$)，仍然没有搜索到更优的新位置时，采蜜蜂放弃当前蜜源成为侦察蜂，然后根据式(17)在解空间随机产生新蜜源.

4 反演分析模型

人工蜂群算法与扩展有限元法的结合可有效地减少反演分析的计算工作量，并能快速探测结构内部存在的缺陷(夹杂). 在笔者前期研究工作中^[5,24]，对应用扩展有限元法求解夹杂、孔洞及裂纹问题的计算结果的准确性进行了验证，结果表明：扩展有限元法在求解这些问题时的精度是有保障的.

如第3节所述，人工蜂群算法在解空间生成大量的随机样本(可能解)，通过蜜蜂的搜索行为不断地更新可能解，最终寻找到最优解.在最优解的搜索过程中，利用扩展有限元法计算适应度函数值来衡量可能解的质量.文中基于“Fortran”算法语言自主研制了相应的程序，可用于识别结构内部存在的夹杂以及缺陷(裂纹、孔洞).反演分析过程如图4所示.

5 数值算例 5.1 算例1：单个圆形孔洞(夹杂)的反演

如图5(a)所示，含圆形孔洞(夹杂)的方形板两边受水平向拉伸载荷$P=1.0$ MPa的作用，板的边长$l=2$ m.数值计算时，假设 板处于平面应变状态，并在板的左边界施加合适的约束条件消除刚体位移，板被离散成39$\times$39的均匀网格(图5(b)). 基体的弹性模量$E=22$ GPa，泊松比$\nu =0.167$；夹杂的弹性模量$E=55$ GPa，泊松比$\nu=0.3$.为获取该结构的真实响应量，在板的上、下边缘处分别布置4个传感器(图5(b))，结构的真实响应通过真实结构的扩展有限元法数值解得到.

本算例中，圆形孔洞(夹杂)有3个待反演参数，圆心坐标$\left( {x_{\rm c} ,y_{\rm c} } \right)$和圆的半径$r_{\rm c}$，如表1所示. 反演分析时，蜜源的数量$N_{\rm Food}=100$，优化极值不变的最大搜索限定次数$N_{\rm Limit}=20$，最大迭代次数$N_{\rm Iter}=100$. 待反演参数的限值为$x_{\rm c} \in \left[ {0.6},{1.8} \right]$，$y_{\rm c} \in \left[ {0.6},{1.8} \right]$，$r_{\rm c} \in \left[ {0.1},{0.5} \right]$.

从表1可以看出，在当前的计算参数下，反演值与真实值非常一致. 反演孔洞时，迭代至99步，目标函数值为4.54$\times$10$^{ - 6}$； 反演夹杂时，迭代至85步，目标函数值为2.70$\times $10$^{-8}$.图6给出了单个圆形孔洞(夹杂)的反演过程.

5.2 算例2：单个椭圆形孔洞(夹杂)的反演

如图7所示，方形板的尺寸、材料参数、载荷条件、边界条件、网格剖分以及传感器布置同5.1节，但这里方形板含一椭圆形的孔洞(夹杂). 结构的真实响应仍通过真实结构的扩展有限元法数值解得到.

本算例中，椭圆形孔洞(夹杂)有5个待反演参数，椭圆中心点坐标$\left( {x_{\rm e} ,y_{\rm e} }\right)$，椭圆半长轴$a$，椭圆半短轴$b$，椭圆方向角$\beta $，如表2所示. 反演分析时，蜜源的数量$N_{\rm Food}=100$，优化极值不变的最大搜索限定次数$N_{\rm Limit} =20$，最大迭代次数$N_{\rm Iter} =300$.待反演参数的限值为$x_{\rm e} \in \left[{0.6},{1.8} \right]$，$y_{\rm e} \in \left[ {0.6},{1.8} \right]$，$a \in \left[ {0.1} ,{0.5} \right]$,$b \in \left[ {0.1},{0.5} \right]$，根据椭圆的对称性，取$\beta \in \left[0 , {90} \right]$.

从表2可以看出，在当前的计算参数下，反演值与真实值吻合较好.反演孔洞时，迭代至291步，目标函数值为5.49$\times $10$^{ -4}$；反演夹杂时，迭代至272步，目标函数值为1.64$\times $10$^{ - 6}$.图8给出了单个椭圆形孔洞(夹杂)的反演过程.相比圆形孔洞(夹杂)的反演分析参数(3个参数)，椭圆形孔洞(夹杂)的反演分析参数(5个参数)较多，其迭代收敛过程也相对慢一些，但椭圆形孔洞(夹杂)更易于适应一些不规则形状的孔洞(夹杂).

5.3 算例3：单个裂纹的反演

如图9所示，方形板的尺寸、材料参数、载荷条件、边界条件、网格剖分以及传感器布置同5.1节，但这里方形板含一裂纹.

本算例中，裂纹有4个待反演参数，裂纹起点坐标$\left( {x_{{\rm cb}} ,y_{{\rm cb}} } \right)$和裂纹终点坐标$\left( {x_{{\rm ce}} ,y_{{\rm ce}} }\right)$，如表3所示. 反演时，蜜源的数量$N_{\rm Food} =100$，优化极值不变的最大搜索限定次数$N_{\rm Limit}=20$，最大迭代次数$N_{\rm Iter} =100$. 待反演参数的限值为$x_{{\rm cb}} \in \left[{0.6},{1.8}\right]$，$y_{{\rm cb}} \in \left[ {0.6},{1.8} \right]$，$x_{{\rm ce}} \in \left[ {0.1},{1.8} \right]$，$y_{{\rm ce}} \in \left[ {0.1},{1.8} \right]$.

从表3可以看出，在当前的计算参数下，反演值与真实值吻合较好. 迭代至第65步时，目标函数值为1.36$\times $10$^{ -4}$.图10给出了内部裂纹的反演过程.

6 结 论

结合扩展有限元法和人工蜂群智能优化算法的优点，建立了结构内部缺陷(夹杂)的反演分析模型，为结构的无损检测技术提供了一条新的途径. 扩展有限元法通过改变水平集函数反映缺陷(夹杂)的位置及大小，避免了反演分析每次迭代过程中的网格重剖分，人工蜂群算法在每次迭代中都采用全局和局部搜索，找到最优解的概率大大增加并可很好地避免局部最优. 扩展有限元法与人工蜂群智能优化算法的结合有效地减少反演分析的计算工作量，能够快速探测结构内部的缺陷(夹杂). 通过若干算例的分析表明：建立的反演分析模型能准确地探测结构内部存在的单个缺陷(夹杂).

文中建立的分析模型能够较易拓展到多个缺陷(夹杂)的反演，为非均质材料中缺陷的无损检测提供了一条新的途径，这也是本文后续的研究重点.