基于半解析有限元法的多层复合环面周向导波计算

摘要

周向超声导波在多层复合圆管的快速缺陷检测与定量材料表征方面具有优秀的应用潜力，结合传感器网络技术，可构建智能材料与结构。精确计算周向导波的频散特性是实现上述目标的关键。传统的频散曲线求解方法大多只适用于各向同性介质且单层环面的情形，而对于工程结构中多层复合圆管则难以甚至无法求解。从弹性波动理论出发，结合COMSOL Multiphysics商用分析平台，发展了半解析有限元（SAFE）方法的计算框架，实现了对周向导波频散关系的精确求解，并在两种典型工业复合圆管结构中得到应用验证。研究表明，所提出的方法可以计算任意材料属性、任意层数以及任意环向截面形状的复合圆管中周向导波的频散曲线，具有重要的工程应用价值。

关键词

周向超声导波; 半解析有限元法; 多层复合圆管; 频散曲线

多层复合圆管结构，相较于单层圆管，在比强度、比模量、线膨胀系数、耐腐蚀性、抗疲劳性和阻尼系数等方面具有显著优势，已发展成为航空航天、能源化工、船舶海洋等工业领域广泛应用的重要承力构件^［1］。因其在高温、超压、腐蚀等复杂苛刻的环境/条件下长时服役，复合圆管结构受疲劳、老化、应力集中、冲击载荷等因素影响，在结构内部及层间界面处易产生多尺度、多形态的缺陷损伤，如分层脱粘、变形、裂纹、腐蚀、凹坑等，进而严重削弱系统整体的承载能力和服役寿命。因此，针对多层复合圆管结构开展快速有效的无损检测尤为重要。超声导波是一种声场能量能够沿波导结构传播的弹性波，并在圆管类结构中可遍及整个壁厚；以超声导波为基础的无损检测与结构健康监测方法，因其具有传播距离长、检测效率高、激发/接收方式灵活、缺陷辨识能力强等优点，已被应用于圆管类结构的缺陷检测与表征。

超声导波检测技术的发展起始于对平板中传播导波的详细研究。Lamb最早在研究无限大平板中传播的正弦波时发现兰姆波现象^［2］，包括振动模态沿平板中性面对称分布的对称S模式波（symmetric mode）与振动模态沿中性面反对称分布的反对称A模式波（antisymmetric mode），并给出其严格的解析描述。Love^［3］在研究各向同性弹性半空间上覆盖另一种层状介质的结构中的波传播过程中发现水平剪切（SH）波现象，平板中SH波位移与兰姆波波场正交，其振动模态沿板厚呈现水平剪切特征。类似于板波（即Lamb波与SH波）^［4］，沿着圆管及圆柱状结构轴向传播的柱面导波亦被发现^［5］。Cawley等^［6］对管道中轴向导波的频散特性和缺陷散射规律开展深入研究，并将轴向导波技术成功应用于油气工业管道的长距离、快速检测。通过沿圆管轴向传播的柱面导波，可从单个超声换能器位置检测获取沿管道轴向约50 m（即正负方向各 25 m）范围内的裂纹、腐蚀等缺陷，并能将其与管道固有的焊接、法兰结构等有效区分^［7］。随着工业用复合圆管的尺寸、层数的不断增加以及管道所处周围介质 （如液体、包覆层等） 的影响，轴向导波在管中传播时会呈现较大的衰减系数，其检测范围和检测能力愈受限制^［8］，难以对圆管结构中的危险局部位置开展更为准确的缺陷检测与环面材料表征。而沿着圆管环面传播的周向导波亦具有扫查整个横截面的能力并对轴向裂纹、径向疲劳裂纹和分层缺陷等具有较高的检测灵敏度，可有效表征复合圆管局域内的多种类型缺陷，相关检测技术因此也获得越来越多的关注与研究。

厘清多层复合圆管中周向导波模式的频散特性（即频散曲线，描述导波相速度或能量速度随着频率改变的变化规律）是开展相关环面检测的首要前提。频散曲线可为选择最合适的导波模式和激励频率范围提供重要依据，以实现特定导波模式（对）的最佳激励和波传播条件；同时其亦能指导超声导波检测信号的分析与处理，实现缺陷的有效检测定位以及对其他结构特征的辨识。然而，工业用多层复合圆管呈现显著的各向异性，尤其是复合材料圆管中各层不同的纤维取向亦带来不可忽略的非均匀性。复杂的材料/结构特征将使得多层复合环面中周向导波问题的解析求解尤为艰巨。

文献中关于管中周向导波问题研究及频散曲线提取工作，相较于轴向导波少，起始于Viktorov明晰圆柱状结构中周向导波与平板结构中导波传播特性的主要差别在于圆柱状结构凸面可支撑表面波的传播，其凹面则不可支撑，而平板的两个表面均可支撑表面波的传播，并可由此推导出单层环面中类兰姆波的特征方程 ^［9］。与平板中的导波传播类似，在均匀各向同性介质环面中，周向导波亦可分为周向水平剪切（CSH）波和周向类兰姆（CLT）波。Liu等^［10］基于频散方程的解析表达式给出不同厚度/半径比下的各向同性介质环面中CLT模式的频散曲线及位移场的数值结果。此外，Zhao等^［11］利用贝塞尔函数扩展位移场解析推导出CSH模式的频散方程，计算了CSH模式的频散曲线以及沿管壁厚度方向的位移和应力场。然而，该方法在较大环半径（即贝塞尔函数阶数较高）情形下较难适用，高阶数、小变量的贝塞尔函数求解时会出现计算效率低、数值上不稳定、可计算频散空间有限等问题^［12］。为了克服上述寻根（root-finding）方法的局限性，Gridin等^［13］在研究圆形钢环中周向导波的传播问题时，引入多种渐近方法替代控制特征方程中的贝塞尔函数以推导出周向导波模式的频散关系，但该类高频渐进分析方法仅仅适用于单层各向同性介质环面。

为进一步研究各向异性材料多层环面中周向导波的频散曲线， Lowe^［14］总结传递矩阵法（transfer matrix method）和全局矩阵法（global matrix method），构建并求解多层结构导波问题的频散方程。然而这两种方法均有局限性，即当结构层数/厚度逐渐增大时，利用传递矩阵法求解超越方程变得既困难又费时，将会出现‘large $f\cdot d$’问题（$f$为频率，$d$为结构厚度）。该类问题是指在导波求解过程中出现频厚积（$f\cdot d$）较大的情形时，会出现矩阵的病态现象，导致计算不稳定且误差增大。全局矩阵法对传递矩阵法进行了修正，可消除传递矩阵法中出现大频厚积情形计算精度损失的问题，但对于层数较多的情况，全局矩阵方法仍需处理较大的矩阵，计算量显著增大。

采用上述两种矩阵方法求解各向异性多层复合圆管中的周向导波频散问题已有相关研究。Towfighi等^［15］引入傅里叶级数扩展未知位移函数，以求解横观各向同性圆管中周向导波的频散曲线。Yu等^［16］基于线性三维弹性，利用Legendre正交多项式级数展开法确定均质正交各向异性圆管中周向导波的频散曲线。上述方法求解各向异性多层圆管结构中周向导波问题通常计算量较大。为进一步提升计算效率，Lugovtsova等^［17］基于比例边界有限元法（scaled boundary finite element method， SBFEM）计算得到铝‒碳纤维多层复合圆管中周向导波的频散曲线及模态。Lin等^［18］将Floquet边界条件方法与扫频有限元方法（sweeping frequency finite element modeling， SFFEM）相结合，高效求解了正交各向异性多层复合材料圆管中的周向导波频散曲线并加以实验验证。

半解析有限元（ semi-analytical finite element， SAFE ）方法被广泛用于求解描述导波传播的特征方程，通过采用有限元网格对波导结构的横截面进行离散，而在波传播方向上（即与波导截面垂直的方向）采用解析描述^［19］，将三维波导问题转化为二维模型进行求解（不改变三维线弹性假设），进而大幅减小计算成本。SAFE方法最初被用于替代一些较为传统的解析求解过程（全局矩阵法、传递矩阵法等）以解决经典的‘large $f\cdot d$’导波问题。Van Velsor^［20］提出SAFE方法数学框架以求解含粘弹性材料的多层各向同性圆管中周向导波的频散曲线。Matuszyk等^［21］提出一种将SAFE方法与完全匹配层（PML）技术相结合的手段，求解了埋地圆管中的周向导波频散曲线。Lin等^［22］基于所提出的场变分方法，改进了已应用于直波导的半解析技术以研究各向异性圆管中周向导波传播问题，并计算了复杂周向截面结构中周向导波的频散曲线。文献报道中，通过独立编程手段执行SAFE方法，已成功求解各向异性^［23］、埋地^［24］或水浸^［25］的圆管结构中周向导波的频散曲线。而本课题组与合作者以及多个国内外研究组则借助于商用有限元分析平台COMSOL Multiphysics在笛卡尔坐标系下执行SAFE方法，已成功求解具有任意材料属性^［26］、任意截面形状^［27］以及任意预应力边界条件^［28］的波导问题，并能求解固埋^［29］、液浸^［30］等情形下的频散曲线。其可望精确计算复杂圆管结构中周向导波的频散曲线，对多层复合环面更具优势。

本文进一步借助COMSOL Multiphysics商用平台在柱坐标系下执行SAFE方法，通过控制特征方程数学推导与弱形式偏微分方程求解模块，提出多层环面的SAFE计算模型，无需编写冗长程序代码，即可高效、准确计算具有任意材料属性、任意曲率、任意层数的多层复合圆管结构中周向导波的频散曲线。通过与文献中已有结果比较，交互验证该数值计算方法的准确性。本文选取工业中两种典型复合圆管结构，即金属‒碳纤维双层环面和多层碳纤维树脂基复合材料环面，计算获取其中传播的所有周向导波模式的频散曲线，为后续对该类结构开展周向检测与表征提供有效基础。

1 理论框架及半解析有限元法

1.1　周向导波的控制方程

对具有任意截面形状波导结构中沿轴向传播的导波问题通常采用笛卡尔坐标系进行描述，而针对具有任意曲率板类结构中周向导波传播问题，需采用柱坐标系进行描述。在柱坐标系（$$r$$，$$\theta$$，$$z$$）中，弹性体应力矩阵的正应力分量为$${\sigma }_{rr}$$，$${\sigma }_{\theta \theta }$$，$${\sigma }_{zz}$$，切应力分量为$${\sigma }_{\theta z}$$，$${\sigma }_{rz}$$，$${\sigma }_{r\theta }$$；而应变包括法向应力长度方向的正应变$${\epsilon }_{rr}$$，$${\epsilon }_{\theta \theta }$$，$${\epsilon }_{zz}$$以及切应力长度方向的切应变$${\epsilon }_{\theta z}$$，$${\epsilon }_{rz}$$，$${\epsilon }_{r\theta }$$。在柱坐标系中，$$u_r$$，$$u_{\theta }$$和$$u_z$$分别表示$$r$$，$$\theta$$和$$z$$方向的变形，由于$$z$$轴方向的分量不随$$z$$变化，其应变张量分量为 $$\begin{array}{l}{\epsilon }_{\mathrm{11}}={\epsilon }_{rr}=\frac{\partial u_r}{\partial r}\\ {\epsilon }_{\mathrm{22}}={\epsilon }_{\theta \theta }=\frac{u_r}{r}+\frac{\mathrm{1}}{r}\frac{\partial u_{\theta }}{\partial \theta }\\ {\epsilon }_{\mathrm{33}}={\epsilon }_{zz}=\mathrm{0}\\ {\epsilon }_{\mathrm{23}}={\epsilon }_{\theta z}=\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathrm{1}}{r}\frac{\partial u_z}{\partial \theta }\\ {\epsilon }_{\mathrm{13}}={\epsilon }_{rz}=\frac{\mathrm{1}}{\mathrm{2}}\frac{\partial u_z}{\partial r}\\ {\epsilon }_{\mathrm{12}}={\epsilon }_{r\theta }=\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\mathrm{1}}{r}\frac{\partial u_r}{\partial \theta }+\frac{\partial u_{\theta }}{\partial r}-\frac{u_{\theta }}{r}\right)\end{array}$$

（1）

取弹性波导中一个无限小的体积微元，由力和位移表示的平衡方程以及本构方程分别为

$$\nabla {\sigma }_{ik}+f_i=\rho \frac{{\partial }^{\mathrm{2}}{u}_i}{\partial t^{\mathrm{2}}},$$

（2）

$${\sigma }_{ik}=C_{ikjl}{\epsilon }_{jl},$$

（3）

式中：$\rho$为介质密度； $\mathit{f}$为体力； $t$为时间变量；$C_{ikjl}$为四阶弹性系数张量。如不计体力，即$f_i=\mathrm{0}$，式（2）可表示为$r$，$\theta$和$z$三个方向上应力分量的运动方程。将式（1）和式（3）代入$r$，$\theta$和$z$三个方向上应力分量的运动方程中可分别得到周向导波沿$r$，$\theta$和$z$方向上位移分量的运动方程。

周向导波是在圆管中沿着圆周$\theta$方向传播的导波，其传播路径闭合；周向导波控制方程的解$u_r$，$u_{\theta }$，$u_z$分别为$r$轴，$\theta$轴和$z$轴三个方向上的位移分量。在均匀各向同性环面中，周向导波分周向水平剪切（CSH）波和周向类兰姆（CLT）波两种。图1显示了圆管中周向导波的质点位移，其中CSH模式的质点运动方向（沿$z$轴）与波传播$\theta$方向垂直，且在$r$轴和$\theta$轴方向的位移分量均为零，即$u_r=u_{\theta }=\mathrm{0}$；CLT模式质点位移在$r-\theta$面内，位移在$z$轴方向的位移分量为零，即$u_z=\mathrm{0}$。

图1  周向导波模式的质点位移示意图

Fig. 1  Schematics of particle displacement for circumferential guided waves modes

1.2　周向导波的半解析有限元（SAFE）求解方法

本文基于COMSOL Multiphysics商用有限元分析平台执行周向导波的SAFE计算方法，通过其中的PDE weak form功能模块求解。对于周向传播的导波，弹性介质中的波场不依赖于$z$方向，其位移矢量$u$只是$r$和$\theta$的函数。在广义平面应变假设下，$z-\theta$平面中的$z$向延伸到无穷大；波传播沿$\theta$方向，并假设其作简谐波形式${\mathrm{e}}^{-\mathrm{i}\omega t}$，由函数${\mathrm{e}}^{\mathrm{i}p\theta }$描述。波导中位移矢量则可描述为

$$u_i\left(r,\theta ,t\right)=U_i\left(r\right){\mathrm{e}}^{\mathrm{i}\left(p\theta -\omega t\right)}$$

（4）

式中：$\mathrm{i}=\sqrt[]{-\mathrm{1}}$是虚数单位；$p$是角波数，最早由Viktorov定义 ^［9］，为环向波数$k$与半径$R$的乘积（$p=kR$，量纲一），在半径$R$处可确定线性相速度值；$\omega =\mathrm{2}\pi f$为角频率；下标$i=r,\theta ,z=\mathrm{1,2},\mathrm{3}$对应全局坐标系的自由度；$U_r(r)$，$U_{\theta }(r)$和$U_z(r)$分别表示各方向上的波振幅。

由于符合广义平面应变假设，$z$轴方向的分量不随$z$变化，即$\partial /\partial z=\mathrm{0}$。将波导中导波模式质点位移分别对柱坐标系下三个方向求偏导可得：

$$\frac{\partial u_i}{\partial r}=\frac{\partial U_i}{\partial r}{\mathrm{e}}^{\mathrm{i}\left(p\theta -\omega t\right)};$$

（5a）

$$\frac{\partial u_i}{\partial \theta }=U_i(\mathrm{i}p){\mathrm{e}}^{\mathrm{i}\left(p\theta -\omega t\right)};$$

（5b）

$$\frac{\partial u_i}{\partial z}=\mathrm{0}.$$

（5c）

将式（5）代入$r$、$\theta$和$z$方向的运动方程，等式两边同除${\mathrm{e}}^{\mathrm{i}\left(p\theta -\omega t\right)}$项，化简可分别得到三个方向上控制微分方程。计算真空环境下环面中周向导波的频散曲线，圆管的内外表面须满足无应力边界条件，即在$r=a$和$r=b$处，${\sigma }_{rr}=\mathrm{0}$，${\sigma }_{r\theta }=\mathrm{0}$，${\sigma }_{rz}=\mathrm{0}$。

控制微分方程及边界条件对应COMSOL Multiphysics中特征值问题的一般输入形式为

$$-{\lambda }^{\mathrm{2}}{\mathit{e}}_a\mathit{u}+\lambda {\mathit{d}}_a\mathit{u}+\nabla \left(\mathit{c}\nabla \mathit{u}+\mathit{\alpha }\mathit{u}\right)-\mathit{\beta }\nabla \mathit{u}-\mathit{a}\mathit{u}=\mathrm{0}$$

（6）

式中：$\lambda$为待求解的特征值，即周向导波模式对应的角波数；$\mathit{u}$为对应的特征向量；系数矩阵${\mathit{e}}_a$、${\mathit{d}}_a$、$\mathit{c}$、$\mathit{\alpha }$、$\mathit{\beta }$、$\mathit{a}$与控制方程参数对应。将控制偏微分方程转换为特征值问题的一般表达式，以矩阵形式表示为

$$\begin{array}{l}\left[\begin{array}{ccc}{C}_{\mathrm{1111}}& C_{\mathrm{1112}}& C_{\mathrm{11113}}\\ C_{\mathrm{1211}}& C_{\mathrm{1212}}& C_{\mathrm{1213}}\\ C_{\mathrm{1311}}& C_{\mathrm{1312}}& C_{\mathrm{1313}}\end{array}\right]\left[\begin{array}{c}\frac{{\partial }^{\mathrm{2}}{U}_r}{\partial r\partial r}\\ \frac{{\partial }^{\mathrm{2}}{U}_{\theta }}{\partial r\partial r}\\ \frac{{\partial }^{\mathrm{2}}{U}_z}{\partial r\partial r}\end{array}\right]+\left(p\right)\frac{\mathrm{i}}{r}\left[\begin{array}{ccc}{C}_{\mathrm{1112}}+C_{\mathrm{1211}}& C_{\mathrm{1122}}+C_{\mathrm{1212}}& C_{\mathrm{1123}}+C_{\mathrm{1213}}\\ C_{\mathrm{1212}}+C_{\mathrm{2211}}& C_{\mathrm{1222}}+C_{\mathrm{2212}}& C_{\mathrm{1223}}+C_{\mathrm{2213}}\\ C_{\mathrm{1312}}+C_{\mathrm{2311}}& C_{\mathrm{1322}}+C_{\mathrm{2312}}& C_{\mathrm{1323}}+C_{\mathrm{2313}}\end{array}\right]\left[\begin{array}{c}\frac{\partial U_r}{\partial r}\\ \frac{\partial U_{\theta }}{\partial r}\\ \frac{\partial U_z}{\partial r}\end{array}\right]+\\ \frac{\mathrm{1}}{r}\left[\begin{array}{ccc}{C}_{\mathrm{1122}}& -C_{\mathrm{1112}}& \mathrm{0}\\ C_{\mathrm{1222}}& -C_{\mathrm{1212}}& \mathrm{0}\\ C_{\mathrm{1322}}& -C_{\mathrm{1312}}& \mathrm{0}\end{array}\right]\left[\begin{array}{c}\frac{\partial U_r}{\partial r}\\ \frac{\partial U_{\theta }}{\partial r}\\ \frac{\partial U_z}{\partial r}\end{array}\right]+\left(-p^{\mathrm{2}}\right)\frac{\mathrm{1}}{r^{\mathrm{2}}}\left[\begin{array}{ccc}{C}_{\mathrm{1212}}& C_{\mathrm{1222}}& C_{\mathrm{1223}}\\ C_{\mathrm{2212}}& C_{\mathrm{2222}}& C_{\mathrm{2223}}\\ C_{\mathrm{2312}}& C_{\mathrm{2322}}& C_{\mathrm{2323}}\end{array}\right]\left[\begin{array}{c}{U}_r\\ U_{\theta }\\ U_z\end{array}\right]+\\ \frac{\mathrm{i}p}{r^{\mathrm{2}}}\left[\begin{array}{ccc}-C_{\mathrm{2212}}+C_{\mathrm{1222}}& -C_{\mathrm{2222}}-C_{\mathrm{1212}}& -C_{\mathrm{2223}}\\ C_{\mathrm{1212}}+C_{\mathrm{2222}}& C_{\mathrm{1222}}-C_{\mathrm{2212}}& C_{\mathrm{1223}}\\ C_{\mathrm{2322}}& -C_{\mathrm{2312}}& \mathrm{0}\end{array}\right]\left[\begin{array}{c}{U}_r\\ U_{\theta }\\ U_z\end{array}\right]+\frac{\mathrm{1}}{r^{\mathrm{2}}}\left[\begin{array}{ccc}-C_{\mathrm{2222}}& C_{\mathrm{2212}}& \mathrm{0}\\ C_{\mathrm{1222}}& -C_{\mathrm{1212}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]\left[\begin{array}{c}{U}_r\\ U_{\theta }\\ U_z\end{array}\right]+\\ \frac{\mathrm{1}}{r}\left[\begin{array}{ccc}{C}_{\mathrm{1111}}-C_{\mathrm{2211}}& C_{\mathrm{1112}}-C_{\mathrm{2212}}& C_{\mathrm{1113}}-C_{\mathrm{2213}}\\ \mathrm{2}{C}_{\mathrm{1211}}& \mathrm{2}{C}_{\mathrm{1212}}& \mathrm{2}{C}_{\mathrm{1213}}\\ C_{\mathrm{1311}}& C_{\mathrm{1312}}& C_{\mathrm{1313}}\end{array}\right]\left[\begin{array}{c}\frac{\partial U_r}{\partial r}\\ \frac{\partial U_{\theta }}{\partial r}\\ \frac{\partial U_z}{\partial r}\end{array}\right]+\rho {\omega }^{\mathrm{2}}{U}_j{\delta }_{ij}=\mathrm{0}\end{array}$$

（7）

式中：${\delta }_{ij}$为Kronecker符号。式（7）为二次特征值问题，为将其转换为线性形式进行简化求解，需引入一个新变量$V_j=p{U}_j$；式（6）待求解的特征向量定义为$\mathit{u}={\left[\begin{array}{llllll}{U}_{\mathrm{1}}& U_{\mathrm{2}}& U_{\mathrm{3}}& V_{\mathrm{1}}& V_{\mathrm{2}}& V_{\mathrm{3}}\end{array}\right]}^{\mathrm{T}}$。根据式（7）可确定COMSOL Multiphysics中待求解的特征值$\lambda =p$，以及相应的系数矩阵。刚度矩阵采用Voigt标注形式。对该特征值问题在频域下求解时，通过指定单个频率$\omega$，可计算获得所研究波导结构中所有导波模式的特征值波数$p$及其特征向量$\mathit{u}={\left[\begin{array}{llllll}{U}_{\mathrm{1}}& U_{\mathrm{2}}& U_{\mathrm{3}}& V_{\mathrm{1}}& V_{\mathrm{2}}& V_{\mathrm{3}}\end{array}\right]}^{\mathrm{T}}$；在多个频率下重复上述计算即可获得各导波模式在所研究频率范围内的频散曲线。

本文中研究多层复合环面波导问题的SAFE计算模型如图2所示。在COMSOL Multiphysics平台中选择二维轴对称组件，采用系数型偏微分方程（coefficient PDE）模块；模型构建采用柱坐标系，坐标变量为$r$，$\theta$，$z$；圆管内径为$a$，外径为$b$，管壁厚度即为$b-a$。设置圆管内外表面满足无应力边界条件；圆管在$z$轴无限延伸，则沿$z$向设置波导截面的周期性边界条件。周向导波沿$\theta$方向传播，波函数用$e^{\mathrm{i}p\theta }$表示，其中角波数$p$与环向波数$k$的关系为$p=kb$。该圆管的SAFE模型采用二阶拉格朗日三角形网格进行截面划分，在$z=\mathrm{0}$和$z=h$边界上设置周期性边界条件以模拟$z$向无限大圆管，如式（8）所定义。在该模型设置中，$h$可选取相对壁厚较小的高度值，减少划分网格数量，并有效避免导波振动模态在较大$h$值情形下的重复分布^［31］。

$$\begin{array}{l}{\mathit{n}}_{z=\mathrm{0}}\cdot {\left(c\nabla \mathit{u}+\alpha \mathit{u}-\gamma \right)}_{z=\mathrm{0}}=\\ {\mathit{n}}_{z=h}\cdot {\left(c\nabla \mathit{u}+\alpha \mathit{u}-\gamma \right)}_{z=h},{\mathit{u}}_{z=\mathrm{0}}={\mathit{u}}_{z=h}.\end{array}$$

（8）

图2  COMSOL Multiphysics平台中建立的SAFE计算模型

Fig. 2  Diagrams for SAFE modelling via COMSOL Multiphysics

通过环面SAFE模型计算提取周向导波的频散关系时，先根据上述求解方法在频率‒波数域求解特征方程的根$p$，由角波数$p$计算出在任意半径$R$处的环向波数$k$（$k=p/R$），然后分别求出多层结构中任意半径处的相速度（$c_{\mathrm{p}}=\omega /k=\omega R/p$）和群速度（$c_{\mathrm{g}}=\partial \omega /\partial k$），以获得其频散曲线。其中，相速度$c_{\mathrm{p}}$从环的内径到外径线性增大。这种相速度随半径增大而增大的现象是在环面沿厚度方向保持相位不变的必要条件，且属于波沿着曲面传播时的独特现象。

通过周向SAFE计算，对环面周向导波的模态特性进行研究，旨在识别可用于筛选的适当模式‒频率组合，并预测其对多层复合环面不同区域缺陷的检测灵敏度。为了选择感兴趣的解决方案，在网格的每个节点位置使用计算模式的时间平均功率流。式中，$P_{{\xi }_i}$为坡印亭（Poynting）矢量在${\xi }_i$方向上分量，为一段时间内的时间平均值；$v_i=-\mathrm{i}\omega U_i$为速度分量；${\tau }_{ji}$为应力张量；下标$i,j=\mathrm{1,2},\mathrm{3}$分别表示$r,\theta ,z$方向。

2 周向导波SAFE计算模型与验证

本节基于1.2节周向SAFE方法求解单层各向同性圆管中周向导波的频散曲线，计算结果与公开文献中所报道的频散曲线结果比较，以验证周向SAFE方法及其计算执行的准确性。所用验证算例中各向同性圆管的几何与材料参数见表1。用于计算环向波数$k$的半径对应空心圆管的外径$b$；圆管的内外径比定义为$\eta =IR/OR=a/b$。

表1  公开文献报道算例中空心圆管的几何与材料参数

Tab. 1  Geometric and material properties for circular tubes reported in the literature

验证算例	内外径比$\eta$	密度$\rho$/（$\mathrm{k}\mathrm{g}\cdot {\mathrm{m}}^{-3}$）	体纵波波速$c_L$/（$\mathrm{m}\cdot {\mathrm{s}}^{-1}$）	体横波波速$c_T$/（$\mathrm{m}\cdot {\mathrm{s}}^{-1}$）
铝制空心圆管	0.1	2 700	6 420	3 020
钢制空心圆管	0.5	7 850	5 960	3 230

首先将周向SAFE方法计算得到的$\eta =\mathrm{0.1}$铝管中CLT模式的频散曲线与Towfighi等论文^［15］中报道的$\stackrel{-}{k}$–$\stackrel{-}{\omega }$频散曲线进行比较，结果如图3所示。其中，量纲一化波数为$\stackrel{-}{k}$（$\stackrel{-}{k}=k(b-a)$），量纲一化角频率为$\stackrel{-}{\omega }$（$\stackrel{-}{\omega }=\omega (b-a)c_{\mathrm{T}}=\omega (b-a)\sqrt[]{\rho /\mu }$）。式中，$c_{\mathrm{T}}$为各向同性介质中传播的体横波波速；$\mu$为介质的拉梅常数。由图3中两种方法计算结果的比较可知，周向SAFE方法计算得到的频散曲线与文献中解析方法计算得到的频散曲线结果高度吻合，在所研究的频率范围内各CLT模式波数对应的最大相对误差均小于1.7 %。

图3  内外径比为0.1铝制空心圆管中CLT导波模式的量纲一角频率与量纲一波数的频散曲线

Fig. 3  Dispersion curves of relation between dimensionless frequency and dimensionless wavenumber for an aluminum circular tube with a ratio of inner diameter to outer diameter of 0.1

根据上述$\stackrel{-}{k}$–$\stackrel{-}{\omega }$频散关系，即可计算周向导波在圆管任意半径R处的相速度（$c(R)=\stackrel{-}{\omega }\cdot R\cdot c_{\mathrm{T}}/\left(\stackrel{-}{k}b\right)$）和群速度（$c_{\mathrm{g}}=\partial \omega /\partial k=\partial \stackrel{-}{\omega }\cdot c_T/\partial \stackrel{-}{k}$）。

类似地，本节采用周向SAFE方法进一步计算了内外径比$\eta =\mathrm{0.5}$的铁制空心圆管中CLT和CSH模式导波的频散曲线，并与Rose等^［12］在公开文献中报道的对应$kd$–$\omega$频散曲线进行比较，结果如图4所示。CSH及CLT导波模式的数值与解析计算结果高度吻合，最大相对误差为2.4 %，有效验证了周向SAFE计算方法的准确性。相较于平板中的兰姆波，环面中CLT波的角频率‒波数频散曲线图原点处出现新的导波模式，即CLT0模式，如图4 b所示。该模式始于原点，然后迅速到达一个峰值，之后降低并与$kd$轴相交于点（$kd=\mathrm{1}-\eta$），这与Rose公开文献^［12］中观察到CLT波的特征一致。CLT模式的产生是由于环的内表面不能支撑表面波，而外表面可支撑表面波传播。因此，任何在环面内表面产生的表面波都会迅速泄漏到环内，然后反射到外表面，从而形成新的导波模式。在平板中，由于板的两个表面都能支持表面波传播，因而不存在该种模式^［10］。当环的内外径比值逐渐减小时，类似于S0板波的环面CLT模式将与其逐渐产生显著差异，该现象说明在具有较小内外径比值的环面上不能支持对称的波位移剖面。

图4  内外径比为0.5钢制空心圆管中周向导波模式的角频率与波数-厚度积的频散曲线

Fig. 4  Dispersion curves of relation between angular frequency and wavenumber-thickness product for a steel circular tube with a ratio of inner diameter to outer diameter of 0.5

3 工业用多层复合圆管结构中的周向导波频散特性

相较于单层各向同性圆管，工业中常用多层复合圆管结构，其在力学性能方面具有更大的裁剪自由度，但在材料组分与结构构型方面亦更为复杂；因此对该类多层复合圆管的加工制造、质量监控及服役性能评估任务均更为艰巨。本文选取了工业中两种典型复合圆管构件，即钢‒碳纤维双层复合圆管和多层碳纤维树脂基复合材料圆管。其中，钢‒碳纤维复合圆管常作为III型复合材料压力容器以用于储存氢气^［17］；因其具有重量轻、可提高压缩气体储存效率等优点，正逐渐代替I型全金属压力容器而应用。而多层碳纤维树脂基复合材料圆管构件具有强度高，寿命长、耐腐蚀，质量轻、密度低等优点，已被用作无人机、机器人、卫星支承等工业设备中的关键承载结构件。沿该类多层复合圆管中传播的周向导波将携带多层环面中的材料/结构/缺陷信息；深入研究并厘清其中周向导波的频散特性则可望形成高效、准确的无损检测与评估技术。在实际检测过程中，通常易测量导波模式沿空心圆管表面传播的相速度或群速度，因此本节选取相速度与频率的关系表征周向导波模式的频散曲线。

基于第2节验证的周向SAFE计算方法，进一步建模计算钢‒碳纤维双层复合圆管中的周向导波。该双层复合圆管由各向同性钢制环面与各向异性碳纤维复合材料环面叠压而成。钢‒碳纤维双层复合圆管的几何与材料参数：钢层的内径为30 mm，外径为32 mm，密度$\rho =$7 850 $\mathrm{k}\mathrm{g}\cdot {\mathrm{m}}^{-\mathrm{3}}$，体纵波波速$c_{\mathrm{L}}=$5 960 $\mathrm{m}\cdot {\mathrm{s}}^{-\mathrm{1}}$，体横波波速$c_{\mathrm{T}}=$3 230 $\mathrm{m}\cdot {\mathrm{s}}^{-\mathrm{1}}$；碳纤维层的内径为32 mm，外径为36 mm，密度$\rho =$1 600 $\mathrm{k}\mathrm{g}\cdot {\mathrm{m}}^{-\mathrm{3}}$。纤维方向沿圆周$\theta$方向，表2给出表征碳纤维复合材料层正交各向异性对称性所需的9个独立弹性常数，然后可依据正交各向异性弹性常数转换得到碳纤维层的刚度矩阵。值得注意的是，复合圆管的外层碳纤维树脂基环面介质采用均匀（homogeneous）假设，认为其正交各向异性刚度矩阵在该层内各个位置是保持一致的。

表2  碳纤维正交各向异性复合材料的弹性常数

Tab.2  Elastic properties for carbon fiber orthotropic composites

弹性模量/$\mathrm{G}\mathrm{P}\mathrm{a}$	泊松比	剪切模量/$\mathrm{G}\mathrm{P}\mathrm{a}$
$E_x=154.1$	${\nu }_{xy}=0.33$	$G_{xy}=7.092$
$E_y=11.4$	${\nu }_{yz}=0.49$	$G_{yz}=3.792$
$E_z=11.4$	${\nu }_{xz}=0.15$	$G_{zx}=7.092$

通过周向SAFE模型计算其特征值波束，图5则给出该双层复合环面中周向导波的相速度频散曲线，包含所研究频率范围内所有CLT模式与CSH模式。其中环向波数及相速度对应的半径选取该环面内层的外径值，图中虚线表示周向水平剪切（CSH）波的频散曲线，实线表示周向类兰姆（CLT）波的频散曲线。不同于平板中SH0模式，环面中CSH0模式出现频散现象。其他高阶CSH模式亦发生频散，随着频率的增加频散降低，均存在截止频率。CLT1模式随着频率增加频散逐渐降低；其他高阶CLT模式亦均存在截止频率。针对频散较低的周向导波模式‒频率组合，可激励对应导波模式，通过缺陷反射信号到达时间以定位环面缺陷，并通过导波透射或反射系数推测缺陷尺寸。针对频散较高的周向导波模式‒频率组合，可利用频厚积微小变化对应显著速度差的关系，结合导波层析成像反演方法^［32］ ，可定量表征环面上的腐蚀减薄特征或应用于其环面材料表征^［33］。

图5  钢‒碳纤维双层复合环面中周向导波的相速度频散曲线

Fig. 5  Phase velocity dispersion curves of circumferential guided waves in steel-CFRP double-layered composite circular tube

通过周向SAFE计算进一步提取特征波数对应的特征向量，可获取对应周向导波的传播模态。典型地，图6分别显示100 kHz频率下钢‒碳纤维双层复合环面中CLT1模式、CSH0模式和CLT2模式的波结构及坡印亭矢量分布情形，其中颜色幅值表示导波模式能流密度的相对幅值（低至高：黑色至白色）；箭头表示导波模式的粒子位移（幅值与方向）；振动模态沿径向壁厚提取$r,\theta ,z$方向的位移分量与各应力分量（下同）。由结果可知，不同周向导波模式的声场能量汇聚于钢‒碳纤维复合圆管的不同区域，可依据能量分布以预测周向导波模式对钢‒碳纤维双层复合环面中不同区域的检测能力。例如，其中CSH0模式的能量主要集中于内环（钢层），可能对界面状态改变不灵敏；而CLT1模式和CLT2模式部分能量局域于层间界面处，对界面处的分层脱粘缺陷将呈现较高的检测灵敏度。同时，沿双层环面的径向壁厚提取了各应力分量的分布情形，亦可为各导波模式的缺陷检测灵敏度提供参考依据。

图6  钢-碳纤维双层复合圆管中周向导波模式在100 kHz频率下典型波结构及坡印亭矢量分布示意图

Fig. 6  Diagrams of typical wave structure for varied circumferential guided wave modes in a steel-CFRP compound circular annuli, calculated by the SAFE method at 100 kHz

为观察不同模式在$r$，$\theta$，$z$方向上的位移分量占比情形，可沿圆管径向厚度（即虚线截取线）提取声场位移分量，可知CLT1模式的径向位移$u_r$占主导，内层（钢层）径向位移$u_r$振幅略高于外层（碳纤维层）径向位移幅值，周向位移$u_{\theta }$和轴向位移$u_z$几乎为0，这与平板中反对称式（A）导波模式的模态类似。而CSH0模式的轴向剪切位移$u_z$占主导，内层（钢层）轴向位移$u_z$幅值较低于外层（碳纤维复合层）轴向位移幅值，其径向位移$u_r$和周向位移$u_{\theta }$几乎为0，这与平板中水平剪切型（SH）导波模式的振动模态相似。CLT2模式的周向位移$u_{\theta }$占主导，其轴向位移$u_z$为0，模态与平板中对称式（S）导波模式相似。

第二类经典环面结构选取12层T800/924碳纤维增强树脂基复合材料圆管结构，其纤维铺层顺序为［（+45°/-45°）3］ s，其中s表示为纤维铺层对称，每层厚度为1 mm，叠压成形后的环面内径和外径分别为 138 mm和150 mm，密度$\rho =$1 500 $\mathrm{k}\mathrm{g}\cdot {\mathrm{m}}^{-\mathrm{3}}$。通过周向SAFE方法可计算该多层复合环面中周向导波的频散曲线，这里注意模型采用非均匀（heterogeneous）假设，即各层定义随纤维取向改变的单方向预浸料的刚度矩阵，其中纤维沿主轴0°方向定义刚度矩阵如表3所示。

表3  单方向碳纤维增强树脂基T800/924复合材料的刚度系数

Tab. 3  Elastic constants of unidirectional T800/924 CFRP composites ( GPa )

$C_{11}$	$C_{12}$	$C_{13}$	$C_{22}$	$C_{23}$	$C_{33}$	$C_{44}$	$C_{55}$	$C_{66}$
11.30	5.45	4.74	164.71	5.45	11.30	6.0	3.28	6.0

借助于商用有限元平台COMSOL Multiphysics，基于所提出的周向SAFE方法获得T800/924多层复合材料环面中周向导波的频散曲线，如图7所示，其中环向波数、相速度的半径选取第四层环面的外径值。图7包含所研究频率范围内所有CLT模式与CSH模式。环面中CSH波的频散特性及振动模态均与板中SH波相似。这里CSH0模式与SH0板波类似，均具有非频散特性。工业中非频散的导波模式能够克服导波频散特性给信号处理带来的困难，最易直接应用于无损检测。高阶CSH模式均发生频散，且随着频率增加，其频散程度降低，均存在对应的截止频率。所研究频率范围内，各CLT模式频散曲线亦给出，其频散曲线并无相交。该现象在多层平板中亦然^［34］。

图7  T800/924多层复合环面中周向导波的相速度频散曲线（虚线表示周向水平剪切（CSH）波，实线表示周向类兰姆（CLT）波）

Fig. 7  Phase velocity dispersion curves of circumferential guided waves in T800/924 multilayered composite circular tube(the dashed line representing the circumferential horizontal shear (CSH) wave, and the solid line representing the circumferential lamb-like (CLT) wave)

平板中水平剪切（SH）模式位移与兰姆波波场正交，其振动模态沿板厚呈现水平剪切特征。类似于平板中的SH模式，周向水平剪切（CSH）模式位移与周向类兰姆（CLT）模式波场正交，其振动模态沿环面厚度呈现水平剪切特征，其粒子位移振动方向（沿$z$轴）与波传播$\theta$方向垂直，且在$r$轴和$\theta$轴方向的位移分量均为零。图8显示了在300 kHz下由周向SAFE方法计算得到的一系列高阶CSH模式的波结构及坡印亭矢量分布情况。CSH模式的波能量集中于多层复合材料环面中的不同区域及层间界面且分布连续，其沿环面传播过程沿厚度方向传递剪切应力，可望对各界面处的分层脱粘缺陷具有较高的检测灵敏度。此外，$\mathrm{C}\mathrm{S}{\mathrm{H}}_{\mathrm{n}}$导波模式的顺序根据其从低频到高频出现在频散曲线中的顺序进行索引，亦与圆管厚度方向上水平剪切振动的节点数量对应。

图8  T800/924多层复合圆管中周向导波模式在300kHz频率下典型的高阶CSH模式的波结构及坡印亭矢量分布

Fig. 8  Diagrams of typical wave structure for varied circumferential guided wave modes in a T800/924 multilayered composite annuli, calculated by the SAFE method at 300 kHz

4 结论

本文提出了一种基于商用有限元平台COMSOL Multiphysics的半解析有限元（SAFE）执行方法，可高效准确地计算多层复合环面中周向导波的频散特性。文中首先从各向异性弹性波导的平衡方程和几何关系出发，结合本构方程，推导了柱坐标系下圆管结构中周向导波的控制微分方程；基于周向SAFE法及特征值问题数学框架，借助COMSOL Multiphysics中弱形式系数偏微分方程模块，实现了对任意圆管中周向导波模式频散关系的精确求解，并交互验证该数值计算方法的准确性。在此基础上，本文将所提出的周向SAFE方法应用于工业中两种典型复合圆管结构，即金属‒碳纤维双层环面和多层碳纤维树脂基复合材料环面，有效计算获取其中传播的所有周向导波模式的频散曲线，并评估了不同周向导波模式应用于环面缺陷检测与材料表征的潜力。

本文所提出的周向SAFE方法可计算任意材料属性、任意曲率及任意层数的复合圆管中周向导波的频散特性，并可推广至具有任意环向截面形状的复杂结构中的导波问题求解。该方法的提出为获取复杂材料与结构中周向导波的频散曲线提供了高效且便捷的手段。在此工作的基础上，可以进一步开展周向导波检测技术研究，并结合超声换能器/传感器布置策略构建自诊断的智能材料与结构，实现对该类工程多层圆管结构的质量控制与健康状态监测。

作者贡献声明

余旭东：提出研究选题、设计研究方案、修订论文、获取研究经费、指导性支持。

秦荣：实施研究过程、理论推导与计算、统计分析、起草论文。

沈海：实施研究过程、数值计算、调研整理文献。

左鹏：设计研究方案、修订论文、指导性支持。

邵照宇：技术与材料支持。