1.一种基于瞬变流频率响应分析的管道泄漏检测方法，其特征在于：包括以下步骤：步骤S1：采集各项管道系统参数；

步骤S2：将步骤S1获得的系统参数代入广义多单元K-V模型进行数值模拟，从而计算得到管道系统的频率响应函数；

步骤S3：利用频率响应函数获取管道系统的FRD；

步骤S4：推导FRD与管道泄漏位置的关系，获得各奇次谐波的幅值|hb|随谐波次数周期性变化规律；

步骤S5：根据各奇次谐波的幅值|hb|随谐波次数周期性变化规律，进行管道泄漏诊断和定位。

2.如权利要求1所述的基于瞬变流频率响应分析的管道泄漏检测方法，其特征在于：步骤S2中包括步骤S2.1建立粘弹性管道的时域控制方程：时域水锤方程中描述加压管道一维瞬态流动的动量方程是：式中，g为重力加速度(m/s2)，A为管道的横截面积(m2)，Q为流量(m3/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)中的蠕变柔度函数利用K-V模型描述：该模型包括一个弹性元件和N个串联的粘弹性元件，弹性元件由弹性模量E0的单个弹簧表示，粘弹性元件由粘度ηk的缓冲器和弹性模量Ek并联的弹簧组成；

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)。

3.如权利要求2所述的基于瞬变流频率响应分析的管道泄漏检测方法，其特征在于：步骤S2中包括步骤S2.2建立粘弹性管道的频域控制方程：根据稳定振荡流的概念，并将压头H、流量Q和应变ε用平均状态下的振荡表示，方程(1)(3)(8)分别改写为：式中，q*、h*和 代表平均状态下流量、压头和延迟应变的振荡分量；R’代表单位长度的阻力； 代表第K个K-V元件的振荡延迟应变；

R’可以用稳定摩擦部分Rs和非稳定摩擦部分Rμs的总和来描述，即：R’＝Rs+Rus      (12)

式中，Rs＝fQ0/(gDA2)为光滑管道湍流；f为Darcy-Weisbach摩擦系数；Rus的表达式取决于所选择的动态摩擦模型：式中，ν是流体的运动粘度；Cus是剪切衰减系数；Cus＝7.41/Rk，K＝log10(14.3/R0.05)；R是雷诺数；

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

式中，q和h为频域中的流量(m3/s)和水头振荡(m)；

将时域的傅立叶变换应用于式(11)并对所得方程实施数学处理：联立方程(7)(15)(16)，将方程改写为：

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

q＝c1 sinh(μL)+c2 cosh(μL)       (18)h＝-Z[c1cosh(μx)+c2sinh(μx)]           (19)式中，c1和c2是根据边界条件确定的常数系数，μ是传播算子，Z是特性阻抗；

n

对于长度为L的均匀管道段，利用上游的边界条件，即x＝0，常数系数c1＝-(1/Z)h ，c2＝qn，其中上标是“n”表示管段的上游端；下游端的解，即x＝L，可表示为：其中，上标“n+1”表示管段的下游端。

4.如权利要求3所述的基于瞬变流频率响应分析的管道泄漏检测方法，其特征在于：步骤S2中还包括步骤S2.3粘弹性管道的频率响应函数：时域内，在管壁粘弹性的瞬态压力波的计算中引入了额外的耗散以及色散，采用多单元K-V模型研究粘弹性对完整粘弹性管道频率响应的影响；

忽略摩擦问题，将阻力系数R的值线性化，将式(20)和(21)中的摩阻设为零，方程式简化为：式中，T是水锤波传播的周期，对于粘弹性管道，T的值总是固定的；

对于由恒定水头储存器和振荡式在线阀门限定的均匀粘弹性管道，振荡阀的传递矩阵为：式中，ΔHvo是直列阀的稳态压头损失； 是稳态下的无量纲阀开度； 是无量纲阀开度振荡的振幅；Zv0是直列阀在稳态下的阻抗，等于2Hv0/Qv0；其中Qv0是通过直列阀的稳态排放；

利用传递矩阵法和管道传递矩阵(22)，求得振荡阀上游侧的复合头振荡为：式中，hd是管道上游端的头部振荡；

对于无摩擦粘弹性管道，将式(23)(24)代入(27)，使用三角函数简化双曲函数，并将hd端除以主动输入项 使其正常化，式(27)变为：式中，hb是管道下游端的标准化头部振荡。

5.如权利要求4所述的基于瞬变流频率响应分析的管道泄漏检测方法，其特征在于：步骤3具体包括：首先，以管道末端电磁阀开度的单频率变化作为输入信号；然后，通过不断改变电磁阀的开度获得管道系统在该频率输入下对应的频率响应；接着，将不同阀门开度下的频率响应参数代入频率响应函数中，得到与之对应的压力振荡幅值；最后，根据频率响应函数，以ω/ωth为横坐标，以断面b处的压力振荡幅值|hb|为纵坐标，得到管道系统的FRD。

6.如权利要求5所述的基于瞬变流频率响应分析的管道泄漏检测方法，其特征在于：步骤S4具体包括以下步骤：S4.1:泄漏状态下管道瞬变流的频域理论模型分析：管道起始断面与末端断面流动参数之间的关系为：

Zc＝PvZb＝Pv(F2PLF1)Za＝UZa     (29)式中：Za为起始断面处的状态向量， Zc为末端断面处的状态向量；F1为上游端至泄漏孔之间管段的场传递矩阵；F2为泄漏孔至下游端之间管段的场传递矩阵；PL、Pv分别为泄漏孔口和末端电磁阀的点传递矩阵；U为总传递矩阵；

进出口断面状态向量之间的场传递矩阵为：

式中，l为管道的长度；Z为管道的特征阻抗；

2 2

Z＝aμ/(gAiω)； R＝fQ0/gDA ；

对于管道中的阀门，阀门进出口断面状态向量之间的关系为：写成扩展状态向量和扩展点矩阵的形式为：

式中，上标L0和R0分别表示左侧和右侧；下标v表示阀门断面； 和 分别称为阀门进出口断面的扩展状态向量；P′v称为阀门的扩展点传递矩阵；ΔHv0和Qv0分别表示恒定流状态下的过阀水力损失和流量；τ0为恒定流状态下的无量纲开度；k为阀门开度变化的幅值；

将管线中的泄漏孔视作固定开度的旁通阀，则主管上泄漏孔上、下游断面状态向量之间的点传递矩阵为：式中，ΔHL0和QL0分别表示恒定流状态下流经泄漏孔的水力损失和泄漏流量，两者之间满足孔口方程 CdL和AL分别为泄漏孔的流量系数和过流面积；

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

对于理想流体，系统断面a和断面c的流动参数之间的关系为：Z′c＝P′vZ′b＝P′vF′2P′vF′1Z′a     (35)将式(30)(33)(34)代入(35)可得：

系统固有角频率为ωth＝πa/2L，对于奇次谐波：ω＝(2m+1)*ωth＝(2m+1)*πa/2L m＝0，1，2，....(37)将式(37)代入式(36)，令 可得：

实际管道系统中，泄漏孔面积占主管道面积的比例很小，可近似认为则式(38)简化为：

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

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

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

7.如权利要求6所述的基于瞬变流频率响应分析的管道泄漏检测方法，其特征在于：步骤S5具体包括：首先，将管道系统FRD中的|hb|取倒数，绘制1/|hb|～m函数；然后，对1/|hb|～m函数进行快速傅里叶变换(FFT)得到管道系统的频谱图；接着，观察频谱图中是否存在主频分量并确定主频所对应的幅值和相位；最后，根据频谱图中的主频和对应的相位确定 即确定泄漏点的位置参数。