一种基于瞬变流频率响应分析的管道泄漏检测方法与流程

本发明涉及粘弹性管道泄漏检测技术领域，特别是涉及一种基于瞬变流频率响应分析的管道泄漏检测方法。

背景技术：

在我国，管道泄漏检测及定位技术的研究已有几十年的历史，但由于管道输送介质、所处环境的多样性及泄漏成因的复杂性，目前尚未形成一种通用、快捷、精确的管道泄漏检测方法。基于瞬变流频域响应分析的管道泄漏检测属于实时模型法，是一种正在发展的新型检测技术，也是目前管道泄漏检测中准确性、可靠性较高的方法。

mpesha于2001年首次提出应用频率响应法来进行泄漏检测，认为系统频响图(frequencyresponsediagram，frd)中存在泄漏导致的次生压力振幅峰值；lee认为frd中不存在泄漏导致的压力振幅峰值，提出通过快关末端旁通阀激发瞬变流，以旁通阀流量变化过程为输入信号来获取系统frd，通过分析frd中奇次谐波幅值随谐波次数的变化规律来进行泄漏检测，但流量变化过程往往难以精确获取；chaudhry以末端振荡阀的开度变化为输入信号来获取系统frd，通过分析frd中偶次谐波幅值随谐波次数的变化规律来进行泄漏检测，但该方法难以准确判断泄漏面积参数；gong采用传递矩阵法直接计算管道系统的frd，通过数值试验对奇次谐波法和偶次谐波法进行了比较，认为奇次谐波法具有更好的鲁棒性，但实际工程中由于泄漏位置未知，无法通过传递矩阵法来直接获取管道系统frd。

针对以上问题，本发明以管道末端电磁阀的开度变化为输入信号，通过计算管道系统的频率响应函数获取管道系统的frd，并基于管道瞬变流的频域理论，推导出管道系统frd中的奇次谐波幅值随谐波次数的变化规律与泄漏位置之间的关系，进而实现管道泄漏的检测与定位。

技术实现要素：

本发明所要解决的技术问题是：为了克服现有技术中的不足，本发明提供一种基于瞬变流频率响应分析的管道泄漏检测方法。

本发明解决其技术问题所要采用的技术方案是：一种基于瞬变流频率响应分析的管道泄漏检测方法，包括以下步骤：

s1：采集管道系统参数

采用实时监测系统和实地测量方法获取试验管道系统的各项参数，包括：流量q(l/s)，测压头h(m)，管道长度l(m)，流体瞬态波的弹性波速ae(m/s)，α为管道约束系数，管道内径d(mm)，流体密度ρ(g/cm³)，管道壁厚e(mm)。

s2：将步骤s1获得的系统参数代入广义多单元k-v模型进行数值模拟，从而得到管道系统的频率响应函数。

s2.1：粘弹性管道的时域控制方程

时域水锤方程中描述加压管道一维瞬态流动的动量方程是：

式中，g为重力加速度(m/s²)，a是管道的横截面积(m²)，q为流量(m³/s)，h是压头(m)，t为时间(s)，x为水锤波沿管道传输的距离(m)，hf是由于摩擦导致的每单位长度的水头损失可以看作是稳态分量hfs和非稳态分量hfu的总和：

hf＝hfs+hfu(2)

在考虑管壁线性粘弹性时，压力管道一维瞬态流动的连续性方程为：

其中，ae为流体瞬态波的弹性波速(m/s)，εr是延迟应变。对于具有线性粘弹性的材料，对某一瞬时应力的应变响应可以写成：

ε(t)＝εe+εr(t)(4)

其中ε是总应变，εr是瞬时弹性应变。对于粘弹性管道，总轴向应变可以描述为：

式中，α为管道约束系数，d为管道内径，ρ为流体密度，e为管道壁厚，下标“0”表示初始条件，j0为瞬时弹性蠕变柔度，j(t)为蠕变柔度函数。式(5)中的第二项对应于延迟应变εr。

式(5)中的蠕变柔度函数利用kevin-voigt(k-v)模型描述：

该模型包括一个弹性元件和n个串联的粘弹性元件。弹性元件由弹性模量e0的单个弹簧表示，粘弹性元件由粘度ηk的缓冲器和弹性模量ek并联的弹簧组成。

利用kevinvoigt(k-v)模型描述蠕变柔度函数，广义kevinvoigt(k-v)模型是一种现象学模型，已被普遍用于描述粘弹性固体的力学行为。

式中，j0等于1/eo，它是第一弹簧的蠕变柔度，jk等于1/ek，它是第k个k-v元素弹簧的蠕变柔度，τk等于ηk/ek，它是第k个k-v元素缓冲器的延迟时间。延迟应变可以写为：

其中，εrk是第k个k-v元素中的延迟应变。将k-v模型应用于粘弹性管道，忽略方程(5)中a、d和e的时间变化，每个延迟应变的时间导数可写成：

其中，c＝αdρg/(2e)。

s2.2：粘弹性管道的频域控制方程

使用稳定振荡流的概念，并假设压头(h)、流量(q)和应变(ε)可以用平均状态下的振荡来表示，方程为(1)、(3)和(8)分别改写为：

式中，q*、h*和代表平均状态下流量、压头和延迟应变的振荡分量，r’代表单位长度的阻力，代表第k个k-v元件的振荡延迟应变。

r’可以用稳定摩擦部分rs和非稳定摩擦部分rμs的总和来描述，即：

r’＝rs+rus(12)

其中rs＝fq0/(gda²)为光滑管道湍流，f为darcy-weisbach摩擦系数。rus的表达式取决于所选择的动态摩擦模型:

其中ν是流体的运动粘度，cus是剪切衰减系数，cus＝7.41/r^k，k＝log10(14.3/r^0.05)，r是雷诺数。

将时间的傅立叶变换应用于方程(9)和(10)得出：

其中q和h为频域中的流量(m³/s)和水头振荡(m)。

将傅立叶变换应用于式(11)并对所得方程实施数学处理，得出

将式(16)代入式(7)，然后将结果代入式(15)，将方程改写为

递推方程(14)和(17)的通解可导出为

q＝c1sinh(μl)+c2cosh(μl)(18)

h＝-z[c1cosh(μx)+c2sinh(μx)](19)

其中，c1和c2是根据边界条件确定的常数系数，μ是传播算子，z是特性阻抗。

对于长度为l的均匀管道段，利用上游的边界条件(即x＝0)，常数系数c1＝-(1/z)hⁿ，c2＝qⁿ，其中上标“n”表示所考虑管段的上游端。该段下游端(即x＝l)的解可表示为：

其中上标“n+1”表示管段的下游端。

s2.3：粘弹性管道的频率响应函数

在时域内，管壁粘弹性对瞬态压力波引入了额外的耗散(阻尼)和色散(相移)。采用多单元k-v模型研究粘弹性对完整粘弹性管道频率响应的影响。

管道粘弹性是本研究的重点，因此摩擦问题将被忽略。将阻力系数r在方程中的值线性化。将式(20)和(21)中的摩阻设为零，方程式简化为

其中，

在模型中给出了粘弹性效应，并且与频率有关。对于弹性管道(没有管壁粘弹性)，t的值总是一致的。

对于由恒定水头储存器和振荡式在线阀门限定的均匀粘弹性管道，振荡阀的传递矩阵可写为

式中，δhvo是直列阀的稳态压头损失，是稳态下的无量纲阀开度，是无量纲阀开度振荡的振幅，zv0是直列阀在稳态下的阻抗，等于2hv0/qv0，其中qv0是通过直列阀的稳态排放。

使用传递矩阵法和管道传递矩阵(22)，上游端的标准化头部振荡为：

式中，hd是管道上游端的标准化头部振荡。

对于无摩擦粘弹性管道，将方程(23)和(24)代入等式(27)，使用三角函数简化双曲函数，并通过将hd端除以主动输入项使其正常化，等式(27)变为

式中，hb是管道下游端的标准化头部振荡。

s3：利用频率响应函数获取管道的频率响应图

首先，以管道末端电磁阀开度的单频率变化作为输入信号；然后，不断改变电磁阀的开度获得到管道系统在该频率输入下对应的频率响应；接着，将不同阀门开度下的频率响应参数代入频率响应函数中，即可得到与之对应的压力振荡幅值；最后，以ω/ωth为横坐标，根据式(28)下游断面处的压力振荡幅值|hb|为纵坐标，即可得到管道系统的frd。

s4：推导frd与管道泄漏位置的关系

s4.1：泄漏状态下管道瞬变流的频域理论模型分析

管道起始断面和末端断面流动参数之间的关系可表示为：

zc＝pvzb＝pv(f2plf1)za＝uza(29)

式中：za为起始断面处的状态向量，zc为末端断面处的状态向量；f1为上游端至泄漏孔之间管段的场传递矩阵；f2为泄漏孔至下游端之间管段的场传递矩阵；pl、pv分别为泄漏孔口和末端电磁阀的点传递矩阵；u为总传递矩阵。

任意管段i的进出口断面状态向量之间的场传递矩阵为：

式中：l为管段长度；z为管道特征阻抗；

对于管线中的阀门，阀门进出口断面状态向量之间的关系可表示为：

写成扩展状态向量和扩展点矩阵的形式为：

式中：上标l0和r0分别表示左侧和右侧；下标v表示阀门断面；和分别称为阀门进出口断面的扩展状态向量；p′v称为阀门的扩展点传递矩阵；δhv0和qv0分别表示恒定流状态下的过阀水力损失和流量；τ0为恒定流状态下的无量纲开度；k为阀门开度变化的幅值。

对于管线中的泄漏孔，将其视作固定开度的旁通阀，则主管上泄漏孔上下游断面状态向量之间的点传递矩阵可表示为：

式中：δhl0和ql0分别表示恒定流状态下流经泄漏孔的水力损失和泄漏流量，两者之间满足孔口方程cdl、al分别为泄漏孔的流量系数和过流面积。

s4.2：理想流体的情况：

对于理想流体，系统断面a和断面c的流动参数之间的关系用扩展状态向量和扩展传递矩阵表示为：

z′c＝p′vz′b＝p′vf2′p′lf1′z′a(35)

将式(30)(33)(34)代入式(35)可得：

系统的水锤波传播的周期为t＝4l/a，其固有角频率为ωth＝πa/2l。对于奇次谐波：

ω＝(2m+1)*ωth＝(2m+1)*πa/2lm＝0，1，2，….(37)

将式(37)代入式(36)，令可得：

实际管道系统中发生不易直观检测到的泄漏时，一般泄漏面积占管道面积的比例很小，因此，可近似则式(38)简化为：

由式(39)可以看出：在管道无泄漏的条件下，各奇次谐波的幅值|hb|与m值无关，均为2δhv0k/τ0；当管道发生泄漏时，各奇次谐波的幅值|hb|不再保持为固定值，而是随着m值的变化而变化。

将式(39)两边取倒数，可得：

式(40)表明：管道发生泄漏条件下，1/|hb|随m值的变化呈现余弦曲线形式的周期性变化，且其变化的频率为初相为幅值为因此，分析管道系统frd中各奇次谐波幅值的倒数1/|hb|随m值变化的规律即可确定泄漏孔的位置参数和面积参数(cdlal)，达到泄漏检测的目的。

s5：基于瞬变流频率响应分析进行管道泄漏判断，具体包括以下步骤：

(1)将管道系统frd中的|hb|取倒数，绘制1/|hb|～m函数；

(2)对1/|hb|～m图像数据进行fft变换，得到管道系统的频谱图；

(3)观察频谱图中是否存在主频分量并确定主频所对应的幅值和相位；

(4)根据频谱图中的主频和对应的相位，确定即泄漏点位置参数。

与常规的采用频率扫描法获取管道frd不同，本发明利用广义多单元kevin-voigt(k-v)模型，建立了频率响应函数的解析表达式，进而通过数值计算获取管道的频率响应图。

本发明的有益效果是：本发明提供的一种基于瞬变流频率响应分析法的管道泄漏检测方法，该方法通过建立管道系统频率响应函数，获取管网的频率响应图，进而通过管道瞬变流的频域理论分析，推导出管道瞬变frd中的奇次谐波幅值随谐波次数的变化规律与泄漏参数之间的关系，之后通过对各奇次谐波的幅值|hb|随谐波次数周期性变化规律的分析来进行管道泄漏的诊断。

附图说明

下面结合附图和实施例对本发明作进一步说明。

图1是本发明的流程图。

图2是试验管道及传感器仪表布置图。

图3是阀门开度与阀门开度系数关系图。

图4是管道系统的frd。

图5是奇次谐波幅值倒数1/|hb|随m值的变化函数图像。

图6是管道系统的频谱图。

图中：1、上游涡街流量计，2、上游压力变送器，3、上游温度传感器，4、泄漏阀，5、支管球阀，6、主管球阀，7、下游温度传感器，8、下游压力变送器，9、下游涡街流量计，10、电磁阀。

具体实施方式

现在结合附图对本发明作详细的说明。此图为简化的示意图，仅以示意方式说明本发明的基本结构，因此其仅显示与本发明有关的构成。

如图1所示，本发明的一种基于瞬变流频率响应分析的管道泄漏检测方法，包括以下步骤：

s1：采集管道系统参数

现场试验管道及传感器仪表布置图如图2所示，其中：上游涡街流量计1和下游涡街流量计9分别用于测量管道上下游的稳态流量和介质流速；上游压力变送器2和下游压力变送器8分别用于测量管道上下游的稳态压头；上游温度传感器3和下游温度传感器7分别用于测量管道系统上下游的实时温度，用于泄漏声波波速的计算；泄漏阀4用于模拟管道的泄漏，位于上游端30米处；支管球阀5和主管球阀6用于试验管道长度的调节，设定支管球阀5处于关闭状态，主管球阀6处于开放状态，即管道长度为100米；末端电磁阀10可实现管道瞬间的开闭，用于瞬态弹性波速的计算；图中箭头表示介质流向。

采用实时监测和实地测量方法获取管道系统的各项参数，如表1所示。

表1管道系统参数

s2：获取管道系统的频率响应函数

将步骤s1获得的管道系统参数代广义多单元k-v模型，得到粘弹性管道下游端的标准化头部振荡函数：

式中，hb为管道下游端压力振荡幅值；i为虚数表示的虚部单位；z为管道稳态下的特性阻抗；g为重力加速度；a为管道截面面积；t为水锤波传播周期；ae为瞬态弹性波速；ω为角频率；l为北侧管道长度。

将式(1)中各参数的求解公式代入并展开，得到管道系统的频率响应函数：

式中，hb为管道下游端压力振荡幅值；δhv0为过阀水力损失；τ0为阀门相对开度；k为阀门开度系数；qv0为过阀平均流量；ql0为泄漏流量损失；δhl0为泄漏水力损失；m为齐次谐波的表示量，取值为0,1,2,...；为泄漏点位置参数。

s3：利用频率响应函数获取管道的频率响应图；

以管道末端电磁阀10开度的单频率变化作为输入信号，通过不断改变电磁阀10的开度获得到管道系统在该频率输入下对应的频率响应。其中，电磁阀10开度分为0.1～1.0共10种开度(1.0代表电磁阀10全开)，阀门开度与阀门开度系数的关系如图3所示。

篇幅所限，只给出阀门开度为0.1时的管道系统参数具体值，如表2所示。

表2阀门开度为0.1时的管道系统参数汇总

将不同阀门开度下的各参数代入频率响应函数中，求得对应的压力振荡幅值hb，以ω/ωth为横坐标，|hb|为纵坐标，绘制管道系统的frd，如图4所示。

s4：推导frd与管道泄漏位置的关系；

由管道系统的frd可知，管道泄漏状态时，frd中的各奇次谐波的幅值|hb|随着m值的变化而变化，为了更清晰地分析|hb|与m值之间的关系，将频率响应函数两边取绝对值的倒数，可得：

式(3)表明：管道泄漏条件下，1/|hb|随m值的变化呈周期性变化，且其变化的频率为初相为幅值为因此，确定管道系统frd中各奇次谐波幅值的倒数1/|hb|随m值的变化规律即可确定泄漏孔的位置参数从而达到泄漏检测的目的。

s5：根据管道系统的frd绘制1/|hb|～m的函数图像，如图5所示。

再对1/|hb|～m函数进行傅里叶变换，得到管道系统的频谱图，如图6所示。对应的主频频率、幅值和相位及由此确定的泄漏点位置参数如表3所示。

表3泄漏检测结果汇总

由表3可知，频谱峰值对应的主频为0.664，因此由于对应的相位为0.873rad/s，位于之间，所以实际泄漏位置即泄漏点距管道上游的距离为33.6米。

以上述依据本发明的理想实施例为启示，通过上述的说明内容，相关的工作人员完全可以在不偏离本发明的范围内，进行多样的变更以及修改。本项发明的技术范围并不局限于说明书上的内容，必须要根据权利要求范围来确定其技术性范围。