DETERMINATION OF WAVE SPEED OF LEAK SOLID SOUND IN FLUID-FILLED PIPELINE UPON LOW-FREQUENCY SIMPLIFIED DISPERSION EQUATIONS
-
摘要:
由充流管道泄漏引致并在管壁上传播的声波被称为泄漏固体声波，其中波速的确定是实现管道泄漏定位的关键。针对适合泄漏定位的低频段的T(0, 1)模态、L(0, 1)模态和F(1, 1)模态，该文提出了一种基于低频简化频散方程的声波波速计算方法。该方法将低频简化的流体荷载和土体荷载施加于管道，从而建立了流体-管道-土体的耦合方程。不同模态的波速可以通过求解耦合方程来得到。以钢管、HDPE管和球墨铸铁管为对象进行算例研究，得到了三种管道在不同流体和土体工况下的频散曲线。结果表明，流体性质对波速的影响较小，而土体性质会对波速产生显著且复杂的影响。该文通过标定实验获取了埋地充水DN200球墨铸铁管的频散曲线。标定结果与所提出的计算方法的结果较为一致，这支撑了所提出方法的工程应用价值。
Abstract:Leak solid acoustic waves are acoustic waves that are initiated by pipe leaks and propagate along pipe walls, and their wave velocities are critical to locating a pipe leak. For the mode T(0, 1), mode L(0, 1) and mode F(1, 1) in the low frequency band suitable for leak location, a method is proposed for calculating the acoustic wave velocity based on the simplified low-frequency dispersion equations. The method applies simplified low-frequency fluid loads and soil loads to pipes, thus establishing the fluid-pipe-soil coupling equations. The wave velocities of the different modes can be obtained by solving the coupled equations. The dispersion curves of three pipe materials, namely steel, HDPE and, ductile iron, are calculated for different fluid and soil conditions. The results show that: fluid properties have a small effect on wave velocities, while soil properties can have a significant and complex effect on wave velocities. The dispersion curve of a buried water-filled DN200 ductile iron pipe is obtained through a calibration experiment. The calibration result is in a good agreement with the result of the calculation method proposed, which supports the engineering application value of the method proposed.
-
管道在城市供水、供燃气、供热和工业生产物料输送中均承担着重要作用。2020年，我国供水管道的平均漏损高达13.4%，燃气管道的漏损率为2.1%[1]。我国在2019年−2021年间共发生管道泄漏事故2049起，占所有管道事故的59.39%[2]，管道泄漏已经成为管线安全运行的重大威胁。
基于管壁固体声波的泄漏识别和定位方法近年来得到了深入的研究与发展。其中，基于固体声波到达不同位置传感器的时间差来确定泄漏位置的时延估计法(time difference of arrival, TDOA)是应用最为广泛的定位算法[3-6]。泄漏声波沿管壁传播的波速是决定该算法定位精度的关键参数。然而，对于实际场景下的复杂工况，难以快速、准确地获取沿管壁传播的声波波速成为了提升管道泄漏定位精度的重大阻碍[7]。
沿管壁传播的声波波速具有频散性和多模态性两大特性。管壁波速的研究工况是由简单到复杂的，即遵循了“非埋地无充流管道—非埋地充流管道—埋地充流管道”的发展历程。由于高频泄漏声波在传播中衰减极快，环频率以下的低频区域是管道泄漏检测重点关注的频段[8]。在环频率以下的低频段，共有T(0, 1)、L(0, 1)、F(1, 1)和S1四种模态，其中，L(0, 1)和F(1, 1)是泄漏定位最为常用的模态。
其论求解方法是常见的沿管壁传播波速的求解方法，可以方便地计算复杂工况下的波速。理论求解方法可以分为圆柱壳波法[9-12]、数值模拟法[13-16]、局部体积波叠加法[17-21]及低频近似求解法[22-27]四个类别。早期的研究将研究对象简化为自由边界的弹性各向同性管道，发展出了针对非埋地无充流管道的圆柱壳波法。为进一步贴合管道的实际场景，研究者们开始进一步考虑管道内流体和管道外土体对沿管壁传播的声波波速的影响。为了解决流体-管壁-土体的复杂耦合系统的管壁波速求解问题，大多数研究采用时空谐波函数ei(kx−ωt)来近似地表达声波沿管道轴向的传播，从而将三维问题简化为二维问题。在这一假设的基础上，以半解析有限元法[13]为代表的数值模拟法得到了较为充分的研究和发展。包含了转移矩阵法和全局矩阵法的局部体积波叠加法，同样采用了时空简谐函数的假设。由于管道泄漏声波的高频信息在传播过程中衰减剧烈，管壁波速研究开始聚焦于低频区域。区别于传统的全频段求解方法，低频近似求解法[22-27]在求解过程中引入低频假设简化了方程，从而获得了仅适用于低频段的管壁波速近似解。表1对比了现有四类理论求解方法。可以发现，现有的四类求解方法均没有能力求解各类工况下的全部低频模态波速，甚至无法很好地求解更受关注的T(0, 1)、L(0, 1)和F(1, 1)三种模态。现在的理论求解方法无法满足管道泄漏定位的需要。
表 1 理论求解方法对比Table 1. Theoretical solution method comparison方法类别 计算工况 计算模态 计算
准确性多层
管壁管外
土体管内
流体T(0, 1) L(0, 1) F(1, 1) S1 圆柱壳波法 是 否 否 是 是 是 否 好 数值模拟法 是 否 否 是 是 是 否 差 局部体积波叠加法 是 是 否 是 是 是 否 中 低频近似求解法 否 是 是 是 是 否 是 好 考虑到当前理论求解方法的不足，本文提出了基于低频简化频散方程的固体声波波速计算方法，实现了T(0, 1)、L(0, 1)和F(1, 1)三种低频模态的波速求解。该方法参考了局部体积波叠加法和低频近似求解法的优点：通过全局矩阵法的局部体积波叠加来描述管壁；采用低频简化荷载来描述流体和土体对传播波速的影响。最后，通过网格化寻根实现了不同模态波速的求解。第1节和第2节分别阐述了低频简化频散方程的构建与求解；第3节以不同流体和土体工况下的钢管、HDPE管、球墨铸铁管为对象进行了算例研究；第4节开展了现场实验验证所提出的方法。
1 低频简化频散方程
图1给出了低频简化频散方程的建立流程。本文首先结合管壁的局部体积波方程和所求解的模态建立了基础方程。然后，在低频假设的基础上推导了流体和土体的应力-位移关系式，使得流体和土体的影响可以以低频简化荷载的形式施加到基础方程上。最终，本文建立了描述流体-管壁-土体耦合系统的低频简化频散方程，可以用于T(0, 1)、L(0, 1)和F(1, 1)三种低频模态的波速求解。该建立流程能够使得建立的低频简化频散方程的矩阵维度不过大，并且保证了矩阵的行数和列数是相等的。
![]()
图 1 低频简化频散方程的建立流程Figure 1. Process for establishing the low frequency simplified dispersion equation埋地充流管道的坐标系与位移符号规定如图2所示：r、x和θ分别为流体-管壁-土体耦合系统的径向、轴向和周向的坐标轴；u、v和w分别为管壁在径向、轴向和周向的位移；ur、ux和uθ分别为土体在径向、轴向和周向的位移。
流体-管壁-土体系统采用了忽略周向耦合的二维耦合。耦合面上的应力符号规定如图3所示：在流体-管壁耦合面，流体对管壁的径向和轴向应力分别定义为pf和τrx；在土体-管壁耦合面，土体对管壁的径向和轴向应力分别被定义为˜σrr和˜σrx。
1.1 基础方程的建立
低频简化频散方程的基础方程是在管壁的局部体积波方程的基础上建立的。管壁局部体积波方程则参考了全局矩阵法的介质层矩阵的建立方法[28-29]。全局矩阵法的核心是，构建并求解一个描述了结构所有界面的位移和应力耦合情况的全局矩阵。对于每一层介质中，该方法假设：
1) 所求解结构中的每一层介质的存在6种局部体积波，其中，向外传播的局部体积波被命名为L+波、SV+波和SH+波，向内传播的局部体积波被命名为L−波、SV−波和SH−波。
2) 6种局部体积波将分别在该层介质的内表面或外表面上产生不同大小的应力和位移。
3) 在内表面或外表面上，6种局部体积波产生的应力和位移的大小是同时受波速、频率、材料性质和体积波强度等多种因素影响的。
本文的管壁局部体积波方程同样采用上述假设，对于单层管壁介质的内表面和外表面，它们的位移场和应力场分别可以表示为由波速、频率和材料性质决定的介质矩阵\boldsymbol{D}与6种体积波的强度向量\bf{Amp}乘积的形式，即：
\left[ \begin{matrix} u \\ v \\ w \\ {{\sigma _{rr}}} \\ {{\sigma _{rx}}} \\ {{\sigma _{r\theta }}} \end{matrix} \right] = \boldsymbol{D}\left[ {\begin{matrix} {{\rm{Amp}}_{\rm{L + }}} \\ {{\rm{Amp}}_{\rm{L - }}} \\ {{\rm{Amp}}_{\rm{SV + }}} \\ {{\rm{Amp}}_{\rm{SV - }}} \\ {{\rm{Amp}}_{\rm{SH + }}} \\ {{\rm{Amp}}_{\rm{SH - }}} \end{matrix}} \right] (1) 式中： u 、 v 和 w 分别为径向位移、轴向位移和周向位移； {\sigma _{rr}} 、 {\sigma _{rx}} 和 {\sigma _{r\theta }} 分别为径向应力、轴向应力和周向应力；\boldsymbol{D}为一个 6 \times 6 维度的矩阵，与全局矩阵法中的介质层矩阵一致，其推导过程和矩阵元素可参见[29]；{\bf{Amp}}为一个 6 \times 1 维度的矩阵，其中的6个元素代表了6种体积波的强度，共同决定了管壁在径向、轴向和周向上的振动强度。由于泄漏定位依赖于沿管壁传播的声波波速，本文无需求解\bf{Amp}中的体积波强度。
当求解不同的模态时，内表面的应力场\boldsymbol{\sigma}_{\rm int}与外表面的应力场\boldsymbol{\sigma}_{\rm ext}需要进行调整和组合来构建基础方程。当求解F(1, 1)模态时，基础方程为：
\left[ {\begin{matrix} {{\sigma _{\rm int\_rr}}} \\ {{\sigma _{\rm int\_rx}}} \\ {{\sigma _{\rm int\_r\theta }}} \\ {{\sigma _{\rm ext\_rr}}} \\ {{\sigma _{\rm ext\_rx}}} \\ {{\sigma _{\rm ext\_r\theta }}} \end{matrix}} \right] = \boldsymbol{D_{\rm F}}\left[ {\begin{matrix} {{\rm{Amp}}_{\rm{L + }}} \\ {{\rm{Amp}}_{\rm{L - }}} \\ {{\rm{Amp}}_{\rm{SV + }}} \\ {{\rm{Amp}}_{\rm{SV - }}} \\ {{\rm{Amp}}_{\rm{SH + }}} \\ {{\rm{Amp}}_{\rm{SH - }}} \end{matrix}} \right] (2) 式中，{\boldsymbol D_{\rm F}}是一个 6 \times 6 的矩阵，它是由外表面介质矩阵{\boldsymbol D_{\rm ext}}的第4行~第6行和内表面介质矩阵\boldsymbol{D}_{\rm int}的第4行~第6行组成的。
在求解L(0, 1)轴向模态时，内外表面应力的径向和轴向分量被保留。体积波则对应地保留了L+、L−、SV+和SV−四种类型的波。最终形成的基础方程为：
\left[ {\begin{matrix} {{\sigma _{\rm int\_rr}}} \\ {{\sigma _{\rm int\_rx}}} \\ {{\sigma _{\rm ext\_rr}}} \\ {{\sigma _{\rm ext\_rx}}} \end{matrix}} \right] = {\boldsymbol D_{\rm L}}\left[ {\begin{matrix} {{\rm{Amp}}_{\rm{L + }}} \\ {{\rm{Amp}}_{\rm{L - }}} \\ {{\rm{Amp}}_{\rm{SV + }}} \\ {{\rm{Amp}}_{\rm{SV - }}} \end{matrix}} \right] (3) 式中，{\boldsymbol D_{\rm L}}为一个 4 \times 4 的矩阵。
在求解T(0, 1)扭转模态时，内外表面应力的周向分量被保留。体积波则仅保留了SH+和SH−两种类型。最终形成的基础方程为：
\left[ {\begin{matrix} {{\sigma _{{\rm{int}}\_{\rm{r}}{\text{θ}} }}} \\ {{\sigma _{{\rm{ext}}\_{\rm{r}}{\text{θ}} }}} \end{matrix}} \right] = \boldsymbol{D_{\rm T}}\left[ {\begin{matrix} {{\rm{Amp}}_{\rm{SH + }}} \\ {{\rm{Amp}}_{\rm{SH - }}} \end{matrix}} \right] (4) 式中，\boldsymbol{D_{\rm T}}为一个 2 \times 2 的矩阵。
1.2 低频简化的流体荷载与土体荷载
计算低频简化的流体荷载和土体荷载实质上是在探讨在耦合面上的应力-位移关系。由于流体不能承受剪切波，低频简化的流体荷载的计算较为简单。假设管壁的粗糙度为0，流体-管壁耦合面的轴向荷载 {\tau _{rx}} 为0，则可以通过零阶的第一类Bessel函数描述耦合界面径向的应力-位移关系。根据文献[23]，流体压强 {p_{\rm f}} 可以通过管壁的径向位移 {u_{\rm f}} 来表示：
{p_{\rm f}} = \frac{{{\omega ^2}{\rho _{\rm f}}}}{{k_{\rm fs}^r}}\frac{{{J_0}(k_{\rm fs}^rR)}}{{{J_0}(k_{\rm fs}^rR)}}{u_{\rm f}} (5) 式中： {J_0} 为0阶的第一类Bessel函数； \omega 为圆频率； {(k_{\rm fs}^r)^2} = k_{\rm f}^2 - k_{\rm s}^2 ， {k_{\rm s}} 为所要求解的模态的波数，流体波数{k_{\rm f}} = \omega /{c_{\rm f}}，流体场波速 {c_{\rm f}} = \sqrt {{B_{\rm f}}/{\rho _{\rm f}}} ， {B_{\rm f}} 和 {\rho _{\rm f}} 分别为流体的体积弹性模量和密度；管壁的半径为 R ； {u_{\rm f}} 为流体-管壁耦合界面处的流体径向位移。
土体被简化为各向同性的弹性介质，根据经典的弹性波理论，仅能传播横波和纵波两种类型的波。假设土体与管壁间不存在周向耦合，因此不需要讨论周向的应力-位移关系。文献[24]给出了土体的轴向应力 {\tilde \sigma _{rx}} 和径向应力 {\tilde \sigma _{rr}} 的表达式。它们可以表达为土体的轴向位移 {u_x} 和径向位移 {u_r} 的形式：
\left[ {\begin{matrix} {{{\tilde \sigma }_{rx}}} \\ {{{\tilde \sigma }_{rr}}} \end{matrix}} \right] = {\boldsymbol{T}}\left[ {\begin{matrix} {{u_x}} \\ {{u_r}} \end{matrix}} \right] (6) 式中，矩阵 {\boldsymbol{T}} 的每一项元素如下所示：
{T_{11}} = \frac{{{\mu _{\rm m}}}}{R}\frac{{k_{\rm ds}^rR{{(k_{\rm rs}^r)}^2}{R^2}}}{{k_{\rm rs}^rRk_{\rm ds}^rR\left[ {{H_0}(k_{\rm rs}^rR)/H_0'(k_{\rm rs}^rR)} \right] + k_{\rm s}^2{R^2}\left[ {{H_0}(k_{\rm ds}^rR)/H_0'(k_{\rm ds}^rR)} \right]}} (7) {T_{12}} = - {\text{i}}{\mu _{\rm m}}{k_{\rm s}}\left\{ {2 - \frac{{{{(k_{\rm rs}^r)}^2}{R^2}{H_0}(k_{\rm ds}^rR)/H_0'(k_{\rm ds}^rR)}}{{k_{\rm rs}^rRk_{\rm ds}^rR\left[ {{H_0}(k_{\rm rs}^rR)/H_0'(k_{\rm rs}^rR)} \right] + k_{\rm s}^2{R^2}\left[ {{H_0}(k_{\rm ds}^rR)/H_0'(k_{\rm ds}^rR)} \right]}}} \right\} (8) {T_{21}} = {\text{i}}{\mu _{\rm m}}{k_{\rm s}}\left\{ {2 - \frac{{{{(k_{\rm rs}^r)}^2}{R^2}{H_0}(k_{\rm ds}^rR)/H_0'(k_{\rm ds}^rR)}}{{k_{\rm rs}^rRk_{\rm ds}^rR\left[ {{H_0}(k_{\rm rs}^rR)/H_0'(k_{\rm rs}^rR)} \right] + k_{\rm s}^2{R^2}\left[ {{H_0}(k_{\rm ds}^rR)/H_0'(k_{\rm ds}^rR)} \right]}}} \right\} (9) {T_{22}} = - \frac{{{\mu _{\rm m}}}}{R}\left\{ {2 + \frac{{k_{\rm rs}^rR{{(k_{\rm rs}^r)}^2}{R^2}\left[ {{H_0}(k_{\rm rs}^rR)/H_0'(k_{\rm rs}^rR)} \right]\left[ {{H_0}(k_{\rm ds}^rR)/H_0'(k_{\rm ds}^rR)} \right]}}{{k_{\rm rs}^rRk_{\rm ds}^rR\left[ {{H_0}(k_{\rm rs}^rR)/H_0'(k_{\rm rs}^rR)} \right] + k_{\rm s}^2{R^2}\left[ {{H_0}(k_{\rm ds}^rR)/H_0'(k_{\rm ds}^rR)} \right]}}} \right\} (10) 式中： \rho_{\rm m} 为土体的密度； \lambda_{\rm m} 和 {\mu}_{\mathrm{m}} 为土体的Lamé常数； t 为时间； {H_0} 为0阶的第二类Bessel函数； H_0' = \partial {H_0}/\partial r ； H_0'' = {\partial ^2}{H_0}/{\partial ^2}r ； {(k_{\rm ds}^r)^2} = k_{\rm d}^2 - k_s^2 ； {(k_{\rm rs}^r)^2} = k_{\rm r}^2 - k_{\rm s}^2 ； {k_{\rm d}} = \omega /{c_{\rm d}} 和 {k_{\rm r}} = \omega /{c_{\rm r}} ； c_{\rm d} 和 c_{\rm r} 分别为土体中的纵波波速和横波波速，可以通过下式进行计算：
{c_{\rm d}} = \sqrt {({\lambda _{\rm m}} + 2{\mu _{\rm m}}){\rho _{\rm m}}} (11) {c_{\rm r}} = \sqrt {{\mu _{\rm m}}/{\rho _{\rm m}}} (12) 1.3 低频简化频散方程的构建
本节采用低频简化的流体荷载和土体荷载对基础方程进行了修正。通过在管道内、外表面分别建立体现流体影响和土体影响的两个耦合应力场(分别命名为 {\boldsymbol \sigma _{{\rm all}\_{\rm int}}} 和 {\boldsymbol \sigma _{{\rm all}\_{\rm ext}}} )，构成了用于求解管壁波速的流体-管壁-土体耦合系统的频散方程。本节主要阐述了最为复杂的埋地充流管道的方程修正，对于无充流或者非埋地的工况，应当相应地去除流体荷载修正环节或土体荷载修正环节。
如图4所示，基础方程修正的流程为：首先需要计算管道内壁和外壁上的位移场，进而通过流体和土体与管壁间的位移耦合，获得流体和土体的位移；然后，利用1.2节提出的流体和土体的位移-应力转换公式，得到流体和土体由于位移耦合而产生的应力，即低频简化荷载；最后，将流体和土体产生的应力通过应力耦合施加到管壁上，便可实现基础方程的修正，即得到频散方程。
1.3.1 管道位移场求解
对于F(1, 1)模态，管壁内表面位移场{{{\boldsymbol{Y}}_{\rm int}}}可以表示为：
{{{\boldsymbol{Y}}_{\rm int}}} {\text{ = }}\left[ {\begin{matrix} {{u_{\rm int}}} \\ {{v_{\rm int}}} \\ {{w_{\rm int}}} \end{matrix}} \right] = {{{\boldsymbol{D}}_{\rm int\_Y}}} \left[ {\begin{matrix} {{{\rm{Amp}}_{\rm L + }}} \\ {{{\rm{Amp}}_{\rm L - }}} \\ {{{\rm{Amp}}_{\rm SV + }}} \\ {{{\rm{Amp}}_{\rm SV - }}} \\ {{{\rm{Amp}}_{\rm SH + }}} \\ {{{\rm{Amp}}_{\rm SH - }}} \end{matrix}} \right] (13) 式中，{\boldsymbol{D_{ \rm int\_Y}}} 为一个 3 \times 6 的矩阵。将内表面参数代入式(1)，可以得到 6 \times 6 的矩阵{{{\boldsymbol{D}}_{\rm int}}} ，该矩阵的第1行~第3行即为{{{\boldsymbol{D}}_{\rm int\_Y}}} 。
外表面位移场 {{{\boldsymbol{Y}}_{\rm ext}}} 可以表示为：
{{{\boldsymbol{Y}}_{\rm ext}}} {\text{ = }}\left[ {\begin{matrix} {{u_{\rm ext}}} \\ {{v_{\rm ext}}} \\ {{w_{\rm ext}}} \end{matrix}} \right] = {{{\boldsymbol{D}}_{\rm ext\_Y}}} \left[ {\begin{matrix} {{{\rm{Amp}} _{\rm L + }}} \\ {{{\rm{Amp}} _{\rm L - }}} \\ {{{\rm{Amp}} _{\rm SV + }}} \\ {{{\rm{Amp}} _{\rm SV - }}} \\ {{{\rm{Amp}} _{\rm SH + }}} \\ {{{\rm{Amp}} _{\rm SH - }}} \end{matrix}} \right] (14) 式中，{{{\boldsymbol{D}}_{\rm ext\_Y}}}为一个 3 \times 6 的矩阵。将外表面参数代入式(1)，可以得到 6 \times 6 的矩阵{{{\boldsymbol{D}}_{\rm ext}}}，该矩阵的第1行~3行即为{{{\boldsymbol{D}}_{\rm ext\_Y}}}。
L(0, 1)模态位移场的计算与F(1, 1)模态是类似的，只需要去除周向位移量，并去除SH+和SH−两种波。管壁内表面位移场 {{{\boldsymbol{Y}}_{{\rm{int}}}}} 可以表示为：
{{{\boldsymbol{Y}}_{\rm int}}}{\text{ = }}\left[ {\begin{matrix} {{u_{\rm int}}} \\ {{v_{\rm int}}} \end{matrix}} \right] = {{{\boldsymbol{D}}_{\rm int\_Y}}} \left[ {\begin{matrix} {{{\rm{Amp}} _{\rm L + }}} \\ {{{\rm{Amp}} _{\rm L - }}} \\ {{{\rm{Amp}} _{\rm SV + }}} \\ {{{\rm{Amp}} _{\rm SV - }}} \end{matrix}} \right] (15) 式中，{{{\boldsymbol{D}}_{\rm int\_Y}}}为一个 2 \times 4 的矩阵。
L(0, 1)模态的外表面位移场{{{\boldsymbol{Y}}_{\rm ext}}}可以表示为：
{{{\boldsymbol{Y}}_{\rm ext}}}{\text{ = }}\left[ {\begin{matrix} {{u_{\rm ext}}} \\ {{v_{\rm ext}}} \end{matrix}} \right] ={{{\boldsymbol{D}}_{\rm ext\_Y}}}\left[ {\begin{matrix} {{{\rm{Amp}} _{\rm L + }}} \\ {{{\rm{Amp}} _{\rm L - }}} \\ {{{\rm{Amp}} _{\rm SV + }}} \\ {{{\rm{Amp}} _{\rm SV - }}} \end{matrix}} \right] (16) 式中，{{{\boldsymbol{D}}_{\rm ext\_Y}}}为一个 2 \times 4 的矩阵。
1.3.2 流体和土体的应力场求解
流体-管壁-土体间采用了忽略周向耦合的二维耦合，因而位移耦合只考虑径向和轴向两个方向。在径向，流体-管壁-土体的位移耦合可以表示为：
{u_{\rm f}}{\text{ = }}{u_{\rm int}} (17) {u_{\rm r}}{\text{ = }}{u_{\rm ext}} (18) 式中： {u_{\rm f}} 为流体的径向位移； {u_{\rm int}} 为管壁内表面的径向位移； {u_{\rm ext}} 为管壁外表面的径向位移； {u_{\rm r}} 为土体的径向位移。
在轴向，管壁与土体的位移耦合可以表示为：
{u_x}{\text{ = }}{v_{\rm ext}} (19) 式中： {v_{\rm ext}} 为管壁外表面的轴向位移； {u_x} 为土体的轴向位移。
可以发现，求解T(0, 1)扭转模态的基础方程式(4)中，不存在径向和轴向的分量。这说明，该模态与流体和土体是解耦的。因此，无需通过流体荷载和土体荷载对基础方程进行修正。对于F(1, 1)模态，由式(13)和式(17)可得到流体的径向位移：
{u_{\rm{f}}} = {{\boldsymbol D_{\rm int\_row1\_F}}} {\bf {Amp}} (20) 式中，{{\boldsymbol D_{\rm int\_row1\_F}}} 为式(13)中矩阵的第1行。
将式(20)代入式(5)，则流体的径向应力 {p_{\rm{f}}} 可以表示为：
{p_{\rm{f}}} = \frac{{{\omega ^2}{\rho _{\rm{f}}}}}{{k_{\rm{fs}}^r}}\frac{{{J_0}(k_{\rm{fs}}^rR)}}{{{J_0}(k_{\rm{fs}}^rR)}} {{\boldsymbol D_{\rm int\_row1\_F}}} {\bf {Amp}} (21) 由式(14)、式(18)和式(19)，可以得到F(1, 1)模态的土体位移场：
{u_r} = {{\boldsymbol D_{\rm ext\_row1\_F}}} {\bf {Amp}} (22) {u_x} = {{\boldsymbol D_{\rm ext\_row2\_F}}} {\bf {Amp}} (23) 将式(22)和式(23)代入式(6)，可以得到F(1, 1)模态的土体应力场：
\left[ {\begin{matrix} {{{\tilde \sigma }_{rx}}} \\ {{{\tilde \sigma }_{rr}}} \end{matrix}} \right] = \boldsymbol T \left[ {\begin{matrix} { {{\boldsymbol D_{\rm ext\_row2\_F}}} {\bf {Amp}} } \\ { {{\boldsymbol D_{\rm ext\_row2\_F}}} {\bf {Amp}} } \end{matrix}} \right] (24) 类似地，将式(15)代入式(17)即可得到针对L(0, 1)模态的流体径向位移：
{u_{\rm{f}}} = {{\boldsymbol D_{\rm int\_row1\_L}}} {\bf {Amp}} (25) 式中， {{\boldsymbol D_{\rm int\_row1\_L}}} 为式(15)中矩阵的第1行。
将式(25)代入式(5)，则L(0, 1)模态的流体径向应力 {p_{\rm{f}}} 可以表示为：
{p_{\rm{f}}} = \frac{{{\omega ^2}{\rho _{\rm{f}}}}}{{k_{\rm{fs}}^r}}\frac{{{J_0}(k_{\rm{fs}}^rR)}}{{{J_0}(k_{\rm{fs}}^rR)}}{{\boldsymbol D_{\rm int\_row1\_L}}} {\bf {Amp}} (26) 由式(16)、式(18)和式(19)可以得到L(0, 1)模态的土体位移场：
{u_r} = {{\boldsymbol D_{\rm ext\_row1\_L}}} {\bf {Amp}} (27) {u_x} ={{\boldsymbol D_{\rm ext\_row2\_L}}} {\bf {Amp}} (28) 将式(27)和式(28)代入式(6)，可以得到L(0, 1)模态的土体应力场：
\left[ {\begin{matrix} {{{\tilde \sigma }_{rx}}} \\ {{{\tilde \sigma }_{rr}}} \end{matrix}} \right] = \boldsymbol T \left[ {\begin{matrix} {{{\boldsymbol D_{\rm ext\_row2\_L}}} {\bf {Amp}} } \\ {{{\boldsymbol D_{\rm ext\_row1\_L}}} {\bf {Amp}} } \end{matrix}} \right] (29) 1.3.3 基于应力耦合的基础方程修正
基础方程的修正实质上是，将流体和土体产生的应力场通过应力耦合条件施加在管壁应力场的过程。由于采用了忽略周向耦合的二维耦合，应力耦合只考虑径向、轴向两个方向。根据应力耦合条件，对于F(1, 1)模态，两个耦合面的应力场 {\boldsymbol \sigma _{\rm all\_ext}} 和 {\boldsymbol \sigma _{\rm all\_int}} 可以表示为：
{\boldsymbol \sigma _{\rm all}} {\text{ = }}\left[ {\begin{matrix} { {\boldsymbol \sigma _{\rm all\_ext}}} \\ { {\boldsymbol \sigma _{\rm all\_int}}} \end{matrix}} \right] = \left[ {\begin{matrix} {{\sigma _{\rm all\_ext\_rr}}} \\ {{\sigma _{\rm all\_ext\_rx}}} \\ {{\sigma _{\rm all\_ext\_r\theta }}} \\ {{\sigma _{\rm all\_int\_rr}}} \\ {{\sigma _{\rm all\_int\_rx}}} \\ {{\sigma _{\rm all\_int\_r\theta }}} \end{matrix}} \right] = \left[ {\begin{matrix} {{\sigma _{\rm ext\_rr}}{\text{ + }}{{\tilde \sigma }_{\rm rr}}} \\ {{\sigma _{\rm ext\_rx}}{\text{ + }}{{\tilde \sigma }_{\rm rx}}} \\ {{\sigma _{\rm ext\_r\theta }}} \\ {{\sigma _{\rm int\_rr}}{\text{ + }}{p_\rm f}} \\ {{\sigma _{\rm int\_rx}}} \\ {{\sigma _{\rm int\_r\theta }}} \end{matrix}} \right] (30) 式中： {{\boldsymbol \sigma _{\rm int}}} 与 {{\boldsymbol \sigma _{\rm ext}}} 为管壁的内表面应力场和外表面应力场； {p_{\rm{f}}} 为流体的径向应力； {\tilde \sigma _{rr}} 和 {\tilde \sigma _{rx}} 分别为土体的径向应力和轴向应力。
将式(2)、式(21)和式(24)代入式(30)，并进行整理，即可得到针对F(1, 1)模态的修正后的基础方程，即F(1, 1)模态频散方程，可以表示为：
\left[ {\begin{matrix} {{\sigma _{\rm all\_ext\_rr}}} \\ {{\sigma _{\rm all\_ext\_rx}}} \\ {{\sigma _{\rm all\_ext\_r\theta }}} \\ {{\sigma _{\rm all\_int\_rr}}} \\ {{\sigma _{\rm all\_int\_rx}}} \\ {{\sigma _{\rm all\_int\_r\theta }}} \end{matrix}} \right] = {{\boldsymbol D_{\rm F\_s}}} \left[ {\begin{matrix} {{\rm{Amp}} _{\rm{L + }}} \\ {{\rm{Amp}} _{\rm{L - }}} \\ {{\rm{Amp}} _{\rm{SV + }}} \\ {{\rm{Amp}} _{\rm{SV - }}} \\ {{\rm{Amp}} _{\rm{SH + }}} \\ {{\rm{Amp}} _{\rm{SH - }}} \end{matrix}} \right] (31) 式中， {{\boldsymbol D_{\rm F\_s}}} 为 6 \times 6 的矩阵。
L(0, 1)模态的频散方程的生成流程与F(1, 1)是相似的，由应力耦合条件而生成耦合面的应力场 {\boldsymbol \sigma _{\rm all\_ext}} 和 {\boldsymbol \sigma _{\rm all\_int}} 可以表示为：
{{\boldsymbol \sigma _{\rm all}}}{\text{ = }}\left[ {\begin{matrix} {{\boldsymbol \sigma _{\rm all\_ext}}} \\ {{\boldsymbol \sigma _{\rm all\_int}}} \end{matrix}} \right] = \left[ {\begin{matrix} {{\sigma _{\rm all\_ext\_rr}}} \\ {{\sigma _{\rm all\_ext\_rx}}} \\ {{\sigma _{\rm all\_int\_rr}}} \\ {{\sigma _{\rm all\_int\_rx}}} \end{matrix}} \right] = \left[ {\begin{matrix} {{\sigma _{\rm ext\_rr}}{\text{ + }}{{\tilde \sigma }_{rr}}} \\ {{\sigma _{\rm ext\_rx}}{\text{ + }}{{\tilde \sigma }_{rx}}} \\ {{\sigma _{\rm int\_rr}}{\text{ + }}{p_\rm f}} \\ {{\sigma _{\rm int\_rx}}} \end{matrix}} \right] (32) 将式(3)、式(26)和式(29)代入式(32)，即可得到针对L(0, 1)模态的频散方程：
\left[ {\begin{matrix} {{\sigma _{\rm all\_ext\_rr}}} \\ {{\sigma _{\rm all\_ext\_rx}}} \\ {{\sigma _{\rm all\_int\_rr}}} \\ {{\sigma _{\rm all\_int\_rx}}} \end{matrix}} \right] = {{\boldsymbol D_{\rm L\_s}}} \left[ {\begin{matrix} {{\rm{Amp}} _{\rm{L + }}} \\ {{\rm{Amp}} _{\rm{L - }}} \\ {{\rm{Amp}} _{\rm{SV + }}} \\ {{\rm{Amp}} _{\rm{SV - }}} \end{matrix}} \right] (33) 式中， {{\boldsymbol D_{\rm L\_s}}}为 4 \times 4 的矩阵。
由于T(0, 1)模态无需进行修正，针对T(0, 1)模态的修正后的基础方程，即频散方程，可以直接由式(4)得到：
\left[ {\begin{matrix} {{\sigma _{\rm all\_ext\_r\theta }}} \\ {{\sigma _{\rm all\_int\_r\theta }}} \end{matrix}} \right] = {{\boldsymbol D_{\rm T\_s}}} \left[ {\begin{matrix} {{\rm{Amp}} _{\rm{SH + }}} \\ {{\rm{Amp}} _{\rm{SH - }}} \end{matrix}} \right] (34) 式中，{{\boldsymbol D_{\rm T\_s}}}为 2 \times 2 的矩阵。
2 波速求解方法
对于L(0, 1)、T(0, 1)和F(1, 1)三种不同的模态，流体-管壁-土体耦合系统中耦合面的应力场 {{{\boldsymbol{\sigma}} _{\rm all}}} 可以表示为 {{{\boldsymbol{D}}_{\rm s}}} 与波强度 {{{\bf{Amp}}}} 的乘积。将流体、管壁和土体三种介质的几何参数和材料参数代入方程，即可得到需要求解的频散方程。基于数值方法，如所示，频散方程的求解主要可以分为频率范围确定、波数确定条件和网格化搜索三个步骤。
1) 频率范围确定
本文的频散曲线求解仅涉及T(0, 1)、L(0, 1)和F(1, 1)三种模态，因而仅需考虑低于环频率的频率范围，从而减少数值计算的计算量。环频率 {f_{\rm r}} 可通过式(35)[30]计算：
{f_{\rm r}} = \frac{1}{{2{\rm \pi} {D_{\rm p}}}}\sqrt {\frac{{{E_{\rm p}}}}{{{\rho _{\rm p}}(1 - \nu _{\rm p}^2)}}} (35) 式中： {E_{\rm p}} 为管壁材料的弹性模量； {\rho _{\rm p}} 为管壁材料的密度； {D_{\rm p}} 为管道直径； {\nu _{\rm p}} 为管道材料的泊松比。可以发现，相同管材的管道的环频率会随着管径的增大而减少，这将导致频散曲线在频率维度上呈现“压缩”的趋势。
2) 波数确定条件
在代入了流体、管壁和土体三种介质的材料参数后， {{{\boldsymbol{D}}_{\rm s}}} 矩阵内的元素便仅由波数和频率两个参数控制。在实际情况下，耦合面处各介质间的应力应当是平衡的，即{{{\boldsymbol{\sigma}} _{\rm all}}}应当为0矩阵。因此，在某一频率下，能够使得{{{\boldsymbol{\sigma}} _{\rm all}}}为0矩阵的波数即为该频率下的管壁传播波数解。以F(1, 1)模态为例，对于某一频率，管壁传播波数求解问题可以转化为：搜寻合适的波数值，该波数值能够使得下面的方程有解。
{{\boldsymbol D_{\rm F\_s}}} \left[ {\begin{matrix} {{\rm{Amp}} _{\rm{L + }}} \\ {{{\rm {Amp}}_{L - }}} \\ {{{\rm {Amp}}_{SV + }}} \\ {{{\rm {Amp}}_{SV - }}} \\ {{{\rm {Amp}}_{SH + }}} \\ {{{\rm {Amp}}_{SH - }}} \end{matrix}} \right]{\text{ = }}\left[ {\begin{matrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{matrix}} \right] (36) 由于每种波在介质中是固定存在的，波强度矩阵{\bf {Amp}}中不存在为0的元素，要使得式(36)有解，则需要：
\det \left( {{{\boldsymbol D_{\rm F\_s}}}} \right) = 0 (37) 因此，在数值求解管壁传播波数的频散曲线中，波数的确定方法为：在某一频率下，针对某一种模态，波数能够使得对应模态的矩阵 {\boldsymbol D_{\rm T\_s}} 、 {\boldsymbol D_{\rm L\_s}} 或 {\boldsymbol D_{\rm F\_s}} 的值为0。
3) 网格化搜索
在频率-波数的二维空间中进行网格划分，对网格的节点进行遍历搜索即可实现频散方程的求解。网格划分的详细参数为：在频率维度上，网格的划分范围为环频率以下，根据经验，网格间距小于20 Hz即可满足求解要求；在波数维度上，网格的划分范围为2\omega /{c_{\rm d}}~100\omega /{c_{\rm d}}，网格间距一般可取 \omega /10 。
对于某一频率下的所有波数节点，当相邻波数节点的满足式(38)时，可以认为频散方程在该频率下的波数解存在于这两个节点之间。
\det { {{{\boldsymbol{D}}_{{\rm s\_s1}}}} } \cdot {\det { {{{\boldsymbol{D}}_{{\rm s\_s2}}}} }} < 0 (38) 管壁波速可以分为群速度 {c_{\rm g}} 和相速度 {c_{\rm p}} 两类，计算公式如下：
{c_{\rm g}} = \frac{{\Delta \omega }}{{\Delta k}} = \frac{{{\rm d}\omega }}{{{\rm d}k}} (39) {c_{\rm p}} = \frac{{\bar \omega }}{{\bar k}} = \frac{\omega }{k}\;\;\;\; (40) 3 计算算例
为验证管壁波速求解方法的可行性，并初步探讨部分因素对管壁波速的影响，本节以常见类型管道为对象，开展了管壁波速传播的理论求解。算例考虑了不同流体和土体工况下DN200的钢管、HDPE管和球墨铸铁管。表2给出了流体-管壁-土体耦合系统的材料参数。
3.1 钢管
由式(35)可得，DN200钢管的环频率为4263 Hz。图5考虑了5种DN200钢管工况，对比了环频率以下模态的相速度和群速度。由于不与流体和土体耦合，T(0, 1)模态波速仅与管道材料特性有关，5种工况下的相速度和群速度均无频散性。由于与管壁的密度相差过大，对于L(0, 1)模态和F(1, 1)模态，水和天然气只能在低频段能管壁波速产生非常小的影响。密度更大的土体对频散曲线的影响较强，其对管道的约束作用在环频率以下的全频段均对管壁波速产生了非常显著的影响。
表 2 材料参数Table 2. Material parameters介质名称 密度
\(kg/m³)弹性模量
\GPa体积模量
\GPa泊松比 HDPE 950.0 0.90 − 0.38 球墨铸铁 7300.0 154.00 − 0.30 Q235钢 7850.0 205.00 − 0.30 水 1000.0 − 2.25 0.50 天然气 1.3 − 1.42×10-4 0.50 土体 2000.0 0.01 − 0.30 3.2 HDPE管
由式(35)可得，DN200 HDPE管道的环频率为837 Hz。图6给出837 Hz以下频段的5种工况的频散曲线，探讨了流体和土体对相速度和群速度的影响。同样，由于不与流体和土体耦合，T(0, 1)模态波速仅与管道材料特性有关，5种工况下的相速度和群速度均无频散性。对于L(0, 1)模态和F(1, 1)模态，水和天然气几乎不能影响管壁的相速度和群速度。由于土体的密度远大于HDPE的密度，埋地对HDPE管道的管壁影响影响非常显著。土体对波速的影响非常复杂：相对于非埋地管道的L(0, 1)模态群速度，埋地管道的波速出现了比较大的波动；在600 Hz以下的频段，土体减小了F(1, 1)模态的相速度，在更高的频段，土体则使得相速度出现了上升；对于F(1, 1)模态群速度，土体显著地增加了管壁波速
3.3 球墨铸铁管
根据式(35)，DN200球墨铸铁管道的环频率为3832 Hz。由于球墨铸铁管道一般不用于输送气体，图7考虑了非埋地无充流、非埋地充水和埋地充水共3种管道工况，对比了环频率以下模态的相速度和群速度。由于T(0, 1)模态与管道外土体和管道内流体不耦合，三种工况下T(0, 1)模态的相速度和群速度均为不随频率变化的固定值。对于L(0, 1)模态和F(1, 1)模态，水对波速的影响很小，这是由流体的密度及管道的密度相差过大决定的。由于土体的密度较大，土体对频散曲线的影响较强：对于L(0, 1)模态的相速度，土体使得2000 Hz以下的波速大幅度下降，频率越低则波速降幅越大；土体使得L(0, 1)模态群速度在的范围内呈现先快速上升、后下降的趋势；对于F(1, 1)模态的相速度，埋地管道在1000 Hz以下的波速明显增大；对于F(1, 1)模态的群速度，土体使得2000 Hz以下的波速出现了小幅度的下降。
4 埋地充水管道的波速标定实验
4.1 实验原理与实验方法
根据管壁固体声波传播理论，环频率以下频段存在着T(0, 1)、L(0, 1)和F(1, 1)三种模态。其中，只有L(0, 1)和F(1, 1)会引发管壁的径向振动，因而任意位置的管壁径向振动均可以表示为这两种模态相叠加的形式。基于图2和图3建立的坐标系，在与声源轴向距离为 x 、管道周向角度为 \theta 的位置，对应于L(0, 1)模态和F(1, 1)模态的管壁声波的表达式分别为：
{V_{\rm L}}{\text{(}}x,t{\text{) = }}{A_{\rm s}}{{\rm e}^{{\text{i}}({k_{\rm L}}x - \omega t)}} (41) {V_{\rm F}}{\text{(}}x{\text{,}}\theta {\text{,}}t{\text{) = }}{A_{\rm s}}{{\rm e}^{{\text{i}}({k_{\rm F}}x - \omega t)}}{{\rm e}^{{\text{i}}\theta }} (42) 式中： {V_{\rm L}} 和 {V_{\rm F}} 分别为两个模态在 x 位置的管壁固体声波； {A_{\rm s}} 为声源引发的固体声波强度； k 为波数； \omega 为圆频率； t 为时间； \theta 为管道的周向坐标。
通过欧拉公式对式(41)和式(42)进行整理和叠加，并忽略式中的虚部，则环频率以下的任意位置的管壁径向振动可以表示为式(43)。这实际上即为布置在管壁上的径向加速度传感器采集到的振动信号。
{V_{\rm radial}}(x,\theta ,t) = {A_{\rm s}}\left[ {\cos ({k_{\rm F}}x - \omega t)\cos \theta + \cos ({k_{\rm L}}x - \omega t)} \right] (43) 实验通过在管道上布置磁致激振器来激振管道。磁致激振器是一种能够以设定的频率对管道施加径向力的装置，其可以产生标定所需要的持续稳定、单一频率的声波。此外，标定实验需要在管壁的不同位置布置4个径向加速度传感器同时采集振动信号。传感器1和传感器2被布置在与激振位置的轴向距离为 {x_1} 的位置1。这两个传感器以轴对称的方式布置，它们的周向角分别被设为 {\theta _1} 和 {\theta _2} ，其中， {\theta _2}{\text{ = }}{\theta _1}{\text{ + }}{180^ \circ } 。传感器3和传感器4则被布置在激振器另一侧的位置2，其与激振器的间距为 {x_2} 。在周向角方面，它们的布置情况和位置1的传感器是相同的。由式(43)可得，标定实验的传感器1~传感器4采集到的振动信号分别可以表示为 {S_1}{\text{(}}t{\text{)}} 、 {S_2}{\text{(}}t{\text{)}} 、 {S_3}{\text{(}}t{\text{)}} 和 {S_4}{\text{(}}t{\text{)}} 。
为了分别求解L(0, 1)模态和F(1, 1)模态的管壁波速，首先需要解决的是，如何将叠加在一起的模态分离开。由于F(1, 1)模态的振动是反对称的，将某一位置轴对称布置传感器的信号进行相加或相减即可得到该位置的L(0, 1)模态和F(1, 1)模态的振动表达式：
{V_{{\rm L},{x_1}}}{\text{ = }}{A_{\rm s}}\cos ({k_{\rm L}}{x_1} - \omega t){\text{ = }}\frac{{{S_1}{\text{(}}t{\text{)}} + {S_2}{\text{(}}t{\text{)}}}}{2} (44) {V_{{\rm F},{x_1}}}{\text{ = }}{A_{\rm s}}\cos ({k_{\rm F}}{x_1} - \omega t)\cos {\theta _1}{\text{ = }}\frac{{{S_1}{\text{(}}t{\text{)}} - {S_2}{\text{(}}t{\text{)}}}}{2} (45) {V_{{\rm L},{x_2}}}{\text{ = }}{A_{\rm s}}\cos ({k_{\rm L}}{x_2} - \omega t){\text{ = }}\frac{{{S_3}{\text{(}}t{\text{)}} + {S_4}{\text{(}}t{\text{)}}}}{2} (46) {V_{{\rm F},{x_2}}}{\text{ = }}{A_{\rm s}}\cos ({k_{\rm F}}{x_2} - \omega t)\cos {\theta _1}{\text{ = }}\frac{{{S_3}(t) - {S_4}(t)}}{2} (47) 式中： {V_{{\rm L},{x_1}}} 和 {V_{{\rm F},{x_1}}} 分别为位置1的L(0, 1)模态和F(1, 1)模态的振动表达式； {V_{{\rm L},{x_2}}} 和 {V_{{\rm F},{x_2}}} 分别为位置2的L(0, 1)模态和F(1, 1)模态的振动表达式。
对于与声源间距已知的位置1和位置2，通过确定声波到达两个位置的时间差，即时延，便可实现管壁波速的计算。互相关法是当前确定时延最为常用的方法之一。该方法的基础是，不同位置的模态振动信号的波形是基本一致的，仅在时间轴上发生平移。因此，将两个位置的固体声信号进行互相关分析，便可以进行时延计算。两种模态的互相关函数峰值所对应的时间即为两种模态的时延。
将时延分别代入式(48)和式(49)，便可以获得所求解的L(0, 1)模态和F(1, 1)模态的波速。
{c_{\rm L}} = \frac{{\left| {{x_1} - {x_2}} \right|}}{{\Delta {t_{\rm L}}}} (48) {c_{\rm F}} = \frac{{\left| {{x_1} - {x_2}} \right|}}{{\Delta {t_{\rm F}}}} (49) 式中： {c_{\rm L}} 和 {c_{\rm F}} 为固体声波沿管壁传播的速度； \Delta {t_{\rm L}} 为位置1与位置2的L(0, 1)模态时延； \Delta {t_{\rm F}} 为位置1与位置2的F(1, 1)模态时延。
基于上述实验原理，在环频率以下的k种频率的激振下，实验可以采集两个位置的轴对称布置的加速度传感器信号。依据传感器编号和激振频率的不同，这四个传感器的信号集合被命名为 \left\{{D}_{1,{f}_{1}}\text{(}n\text{)}，\cdots,{D}_{1,{f}_{k}}\text{(}n\text{)}\right\} 、 \left\{{D}_{2,{f}_{1}}\text{(}n\text{)}，\cdots,{D}_{2,{f}_{k}}\text{(}n\text{)}\right\} 、 \left\{{D}_{3,{f}_{1}}\text{(}n\text{)}，\cdots,{D}_{3,{f}_{k}}\text{(}n\text{)}\right\} 和 \left\{{D}_{4,{f}_{1}}\text{(}n\text{)}，\cdots,{D}_{4,{f}_{k}}\text{(}n\text{)}\right\} 。下面阐述了信号处理方法的5个步骤。
1) 信号降噪
该步骤采用以激振频率为中心、100 Hz带宽的带通滤波器对四个传感器的信号集合进行带通滤波。该处理可以避免环境振动对后续波速计算的干扰。
2) 样本分割
经过信号降噪处理的4个传感器信号集合中的数据将被分割为时长1s的样本。样本分割可以避免部分异常信号对波速计算结果的影响。经过分割后的样本被命名为 Sampl{e_{i,j,{f_k}}}(n) ，其中， i 为传感器编号， j 为样本编号， {f_k} 为激振频率， n 为样本中的时序数据编号。
3) 模态分离
对于每一个激振频率，将同一时间的由4个传感器采集的样本代入式(44)~式(47)进行模态分离，可以得到每个频率下的L(0, 1)模态和F(1, 1)模态的振动信号。 {L_{1,j,{f_k}}}(n) 和 {F_{1,j,{f_k}}}(n) 分别为在 {f_k} 频率的激振下的位置1处的L(0, 1)模态和F(1, 1)模态的第 j 个样本； {L_{2,j,{f_k}}}(n) 和 {F_{2,j,{f_k}}}(n) 分别为在 {f_k} 频率的激振下位置2处的L(0, 1)模态和F(1, 1)模态的第 j 个样本。
4) 基于互相关的时延估计
对不同模态的信号分别计算互相关函数，出现峰值的时间即为所求解的时延。对于每一个激振频率和每一种模态，不同样本的时延会被对比，并剔除显著偏离的异常结果。
5) 频散曲线
将L(0, 1)模态和F(1, 1)模态的样本时延结果分别代入式(48)和式(49)，即可得到每一个样本的管壁波速结果。将同一个激振频率下的波速结果取均值，可以得到该频率的管壁波速。将不同频率的管壁波速进行连线，即可获得实验标定的管壁波速频散曲线。
4.2 实验平台与实验工况
受限于实验条件，本文只进行了埋地充水球墨铸铁管道的标定实验来验证所提出的理论求解方法。实验平台的主体为一根总长度约40 m的埋地充水管道。管道为DN200的球墨铸铁管，埋地深度为0.8 m，内压保持在0.3 MPa。如图8所示，该管道共设有4个阀井。其中，阀井A、B和D可布置加速度传感器，而阀井C可布置磁致激振器。
磁致激振器(图9(a))被放置在阀井C中来在管壁上产生不同频率的声波。在振动信号采集设备方面，实验采用了NI9234数采设备，采样率为10 240 Hz。加速度传感器的型号为Lance LC0115T，其灵敏度为4949 mV/g，频响范围为0.5 Hz~2000 Hz。在传感器布置方面，由于实验仅需在两个阀井布置传感器，本文选择了间距较大且位于激振位置两侧的阀井A和阀井D来布置传感器，阀井B并未用于实验。如图9(b)所示，实验在两个阀井中分别以轴对称的形式布置了2个加速度传感器。在实验工况方面，以100 Hz为间隔采集了200 Hz~2000 Hz间的19个标定实验工况。
4.3 信号处理与标定结果
对传感器采集到的信号进行了时域和频域的初步分析。以距离激振位置1.7 m的传感器为例，图10(a)~图10(d)给出了500 Hz和1000 Hz激振频率下的原始信号的时域图和频域图。可以发现，由于实验场地的干扰噪声较多，磁致激振器产生的500 Hz和1000 Hz的管壁振动被环境振动淹没。因此，对传感器采集到的原始信号进行滤波降噪是十分必要的。图10(e)~图10(h)中的滤波信号的时域波形接近正弦函数，这与磁致激振器的正弦激振方式是一致的。滤波信号的频域峰值分别位于500 Hz和1000 Hz，这与激振频率也是完全一致的。因此，可以认为滤波信号即为由激振产生的管壁振动信号。
通过4.1节提出的信号处理方法，可以得到不同工况下的所有样本的波速标定结果。本文首先进行了离散性分析，以说明实验结果的有效性。为了直观地展示实验结果的离散性，对于L(0, 1)模态和F(1, 1)模态的每一个激振频率，本文统计了与均值的差值小于均值的30%的样本个数占比。如图11所示，对于两种模态的所有频率，绝大多数样本计算得到的波速与均值的差值小于均值的30%。这表明实验结果的离散性较小，从而支撑了结果的可信度。
![]()
图 11 差值小于均值的30%的样本个数占比Figure 11. Proportion of samples with difference less than 30% of the mean如图12所示，每一个频率下的波速平均值被连接从而构成了实验标定的波速频散曲线。这些实验标定结果与第3.3节的理论计算结果进行了对比，频散曲线的整体趋势较为一致，这表明所提出的低频简化频散方程法是较为可靠的。
5 结论
对于基于声信号的管道泄漏定位方法，沿管壁传播的声波波速是最为重要的参数之一。由于环频率以下的低频段信号更适合进行泄漏定位，本文提出了一种该频段T(0, 1)模态、L(0, 1)模态和F(1, 1)模态的泄漏固体声波波速计算方法，重要结论如下：
(1) T(0, 1)模态所对应的周向不与管道内流体和管道外土体耦合，在不同工况下的T(0, 1)模态仅由管道材料决定。
(2) 由于管道内流体与管壁的密度相差较大，流体不会对L(0, 1)模态和F(1, 1)模态的波速产生显著影响。
(3) 对于L(0, 1)模态和F(1, 1)模态的波速，管道外土体会产生非常明显且复杂的影响。
(4) 本文理论方法计算的埋地充水管道频散曲线与现场测试的结果基本吻合，表明所提出的方法具有工程应用价值
-
![]()
图 1 低频简化频散方程的建立流程
Figure 1. Process for establishing the low frequency simplified dispersion equation
![]()
图 11 差值小于均值的30%的样本个数占比
Figure 11. Proportion of samples with difference less than 30% of the mean
表 1 理论求解方法对比
Table 1 Theoretical solution method comparison
方法类别 计算工况 计算模态 计算
准确性多层
管壁管外
土体管内
流体T(0, 1) L(0, 1) F(1, 1) S1 圆柱壳波法 是 否 否 是 是 是 否 好 数值模拟法 是 否 否 是 是 是 否 差 局部体积波叠加法 是 是 否 是 是 是 否 中 低频近似求解法 否 是 是 是 是 否 是 好 
下载: 导出CSV
表 2 材料参数
Table 2 Material parameters
介质名称 密度
\(kg/m³)弹性模量
\GPa体积模量
\GPa泊松比 HDPE 950.0 0.90 − 0.38 球墨铸铁 7300.0 154.00 − 0.30 Q235钢 7850.0 205.00 − 0.30 水 1000.0 − 2.25 0.50 天然气 1.3 − 1.42×10-4 0.50 土体 2000.0 0.01 − 0.30
下载: 导出CSV
-
[1] 中华人民共和国住房和城乡建设部. 2020年城乡建设统计年鉴 [R]. 北京: 中华人民共和国住房和城乡建设部, 2021. Ministry of Housing and Urban-Rural Development of the People's Republic of China. 2020 urban and rural construction statistical yearbook [R]. Beijing: Ministry of Housing and Urban-Rural Development of the People's Republic of China, 2021.
[2] 北京地下管线综合管理研究中心. 2021年全国地下管线事故统计分析报告 [R]. 北京: 中国城市规划协会地下管线专业委员会, 2021. Beijing Underground Pipeline Comprehensive Management Research Center. 2021 national underground pipeline accident statistical analysis report [R]. Beijing: Underground Pipeline Professional Committee of China Urban Planning Association, 2021. (in Chinese)
[3] LI S Z, GONG C Y, LIU Z L. Field testing on a gas pipeline in service for leak localization using acoustic techniques [J]. Measurement, 2021, 182: 109791. doi: 10.1016/j.measurement.2021.109791
[4] ABDULSHAHEED A, MUSTAPHA F, GHAVAMIAN A. A pressure-based method for monitoring leaks in a pipe distribution system: A review [J]. Renewable and Sustainable Energy Reviews, 2017, 69: 902 − 911. doi: 10.1016/j.rser.2016.08.024
[5] BOAZ L, KAIJAGE S, SINDE R. An overview of pipeline leak detection and location systems [C]// Proceedings of the 2nd Pan African International Conference on Science, Computing and Telecommunications. Arusha: IEEE, 2014: 133-137.
[6] 李俊花, 孙昭晨, 崔莉. 基于BP神经网络原理的长输管道泄漏点定位及其实验研究[J]. 工程力学, 2010, 27(8): 169 − 173. LI Junhua, SUN Zhaochen, CUI Li. Leakage localization for pipelines based on BP neural network and experimental verification [J]. Engineering Mechanics, 2010, 27(8): 169 − 173. (in Chinese)
[7] 范重, 张康伟, 张郁山, 等. 视波速确定方法与行波效应研究[J]. 工程力学, 2021, 38(6): 47 − 61. doi: 10.6052/j.issn.1000-4750.2020.05.S010 FAN Zhong, ZHANG Kangwei, ZHANG Yushan, et al. Study on apparent wave velocity calculation method and on travelling wave effect [J]. Engineering Mechanics, 2021, 38(6): 47 − 61. (in Chinese) doi: 10.6052/j.issn.1000-4750.2020.05.S010
[8] 巩晨阳. 连续焊接燃气管道泄漏声发射源定位技术研究[D]. 上海: 同济大学, 2020. GONG Chenyang. Research on location technology of acoustic emission source of continuous welded gas pipeline leakage [D]. Shanghai: Tongji University, 2020. (in Chinese)
[9] NIUTTA C B, TRIDELLO A, PAOLINO D S, et al. Residual properties in damaged laminated composites through nondestructive testing: a review [J]. Materials, 2021, 14(16): 4513. doi: 10.3390/ma14164513
[10] ZHENG M F, MA H W, LYU Y, et al. On the dispersion of cylinder guided waves propagating in a multilayer composite hollow cylinder made of anisotropic materials [J]. Aerospace Science and Technology, 2019, 95: 105432. doi: 10.1016/j.ast.2019.105432
[11] HAYWOOD-ALEXANDER M, DERVILIS N, WORDEN K, et al. Structured machine learning tools for modelling characteristics of guided waves [J]. Mechanical Systems and Signal Processing, 2021, 156: 107628. doi: 10.1016/j.ymssp.2021.107628
[12] SZLASZYNSKI F, LOWE M J S, HUTHWAITE P. Short range pipe guided wave testing using SH0 plane wave imaging for improved quantification accuracy [J]. Sensors, 2022, 22(8): 2973. doi: 10.3390/s22082973
[13] YANG Z Y, WU Z J, ZHANG J Q, et al. Acoustoelastic guided wave propagation in axial stressed arbitrary cross-section [J]. Smart Materials and Structures, 2019, 28(4): 045013. doi: 10.1088/1361-665X/aadb6e
[14] JOSEPH R, LI L F, HAIDER M, et al. Hybrid SAFE-GMM approach for predictive modeling of guided wave propagation in layered media [J]. Engineering Structures, 2019, 193: 194 − 206. doi: 10.1016/j.engstruct.2019.04.082
[15] SARAVANAN T J. Convergence study on ultrasonic guided wave propagation modes in an axisymmetric cylindrical waveguide [J]. Mechanics of Advanced Materials and Structures, 2022, 29(13): 1856 − 1873. doi: 10.1080/15376494.2020.1842949
[16] 徐涛龙, 邵常宁, 兰旭彬, 等. 粒子法和离散元法在管土耦合分析中的应用[J]. 工程力学, 2022, 39(S1): 239 − 249. XU Taolong, SHAO Changning, LAN Xubin, et al. Application of particle method and discrete element method in pipe-soil coupling analysis [J]. Engineering Mechanics, 2022, 39(S1): 239 − 249. (in Chinese)
[17] GHAVAMIAN A, MUSTAPHA F, BAHARUDIN B T H T, et al. Detection, localisation and assessment of defects in pipes using guided wave techniques: a review [J]. Sensors, 2018, 18(12): 4470. doi: 10.3390/s18124470
[18] PARRINELLO A, KESOUR K, GHIRINGHELLI G L, et al. Diffuse field transmission through multilayered cylinders using a Transfer Matrix Method [J]. Mechanical Systems and Signal Processing, 2020, 136: 106514. doi: 10.1016/j.ymssp.2019.106514
[19] LIU C C, YU J G, ZHANG B, et al. Analysis of Lamb wave propagation in a functionally graded piezoelectric small-scale plate based on the modified couple stress theory [J]. Composite Structures, 2021, 265: 113733. doi: 10.1016/j.compstruct.2021.113733
[20] GAO J, LYU Y, ZHENG M F, et al. Application of state vector formalism and Legendre polynomial hybrid method in the longitudinal guided wave propagation analysis of composite multi-layered pipes [J]. Wave Motion, 2021, 100: 102670. doi: 10.1016/j.wavemoti.2020.102670
[21] GEORGIADES E, LOWE M J S, CRASTER R V. Leaky wave characterisation using spectral methods [J]. The Journal of the Acoustical Society of America, 2022, 152(3): 1487 − 1497. doi: 10.1121/10.0013897
[22] 林亨, 吴冬雁, 赵俊亮. 考虑淤砂层的理想流体层中悬浮隧道管体动力水荷载研究—SV波[J]. 工程力学, 2022, 39(8): 114 − 121. doi: 10.6052/j.issn.1000-4750.2021.04.0295 LIN Heng, WU Dongyan, ZHAO Junliang. Effect of sediment model on dynamic pressures of submerged floating tunnel due to SV wave incidence [J]. Engineering Mechanics, 2022, 39(8): 114 − 121. (in Chinese) doi: 10.6052/j.issn.1000-4750.2021.04.0295
[23] MUGGLETON J M, KALKOWSKI M, GAO Y, et al. A theoretical study of the fundamental torsional wave in buried pipes for pipeline condition assessment and monitoring [J]. Journal of Sound and Vibration, 2016, 374: 155 − 171. doi: 10.1016/j.jsv.2016.03.035
[24] GAO Y, LIU Y Y, MUGGLETON J M. Axisymmetric fluid-dominated wave in fluid-filled plastic pipes: Loading effects of surrounding elastic medium [J]. Applied Acoustics, 2017, 116: 43 − 49. doi: 10.1016/j.apacoust.2016.09.016
[25] LIU Y, HABIBI D, CHAI D, et al. A numerical study of axisymmetric wave propagation in buried fluid-filled pipes for optimizing the vibro-acoustic technique when locating gas pipelines [J]. Energies, 2019, 12(19): 3707. doi: 10.3390/en12193707
[26] XIAO R, JOSEPH P F, MUGGLETON J M, et al. Limits for leak noise detection in gas pipes using cross correlation [J]. Journal of Sound and Vibration, 2022, 520: 116639. doi: 10.1016/j.jsv.2021.116639
[27] LU P, SHENG X Z, GAO Y, et al. Influence of shear effects on the characteristics of axisymmetric wave propagation in a buried fluid-filled pipe [J]. Chinese Journal of Mechanical Engineering, 2022, 35(1): 74. doi: 10.1186/s10033-022-00710-7
[28] LING E H, RAHIM R H A. A review on ultrasonic guided wave technology [J]. Australian Journal of Mechanical Engineering, 2020, 18(1): 32 − 44. doi: 10.1080/14484846.2017.1373385
[29] PAVLAKOVIC B, LOWE M. Disperse user’s manual [R]. London: Imperial College, 2013.
[30] MOORE S. A review of noise and vibration in fluid-filled pipe systems [C]// Proceedings of the Acoustics. Brisbane, Australia, 2016.
计量
- 文章访问数: 214
- HTML全文浏览量: 69
- PDF下载量: 44
