4. 管道瞬变流数学模型

泄漏检测与压力瞬变分析的核心，是把管道内“压力—流量—密度”随时间和位置的变化用一组偏微分方程描述清楚。工程上通常采用一维瞬变流模型：沿管道轴向取坐标 x，以稳态流动方向为正方向，取时间 t，未知量为压力 p(x,t)、流速 v(x,t) 和密度 ρ(x,t)。三个基本方程分别来自质量守恒、动量守恒和流体状态关系，缺一不可。

连续性方程：质量守恒的微分表达

取长度为 dx 的管段微元，管内质量变化率等于流入与流出质量流量之差。对圆管截面面积 A，连续性方程可写为：

∂ρ/∂t + ρ·∂v/∂x + v·∂ρ/∂x = 0
等价形式： (1/ρ)·(∂ρ/∂t + v·∂ρ/∂x) + ∂v/∂x = 0

式中第一项反映密度随时间的变化，后两项反映对流与流速梯度引起的质量积聚。对于液体管道，密度变化很小，常用体积弹性模量 K 把密度变化换算为压力变化，从而把连续性方程改写为以压力和水头为变量的形式：

∂p/∂t + v·∂p/∂x + ρ·a²·∂v/∂x = 0
a = √(K/ρ)   （a 为压力波速，m/s）

这里的波速 a 是瞬变分析最关键的参数之一。它并非声速的简单套用，而要计入管壁弹性与流体压缩性，常用公式为 a = √[(K/ρ) / (1 + K·D/(E·e)·C)]，其中 D 为管径、e 为壁厚、E 为管材弹性模量、C 为管道约束系数。以常见成品油管道为例，K 约 1.3×10⁹ Pa、ρ 约 830 kg/m³，若忽略管壁弹性，波速约 1250 m/s；计入钢管弹性后通常降到 1000～1100 m/s 量级。这个量级直接决定了泄漏压力波到达两端传感器的时间差，是定位计算的基础。

动量方程：力与加速度的平衡

动量方程描述管段内流体所受压力梯度、重力分量、摩阻与惯性之间的平衡。以测压管水头 H 和流量 Q 为变量，常见形式为：

∂H/∂x + (1/(gA))·∂Q/∂t + (f·Q·|Q|)/(2g·D·A²) = 0

三项依次代表压力梯度项、惯性项和摩阻项。摩阻项中的 f 为达西摩阻系数，用 |Q|·Q 保证流向反转时摩阻方向随之改变。若管道存在高程变化，还需在压力项中叠加地形坡度影响。动量方程是泄漏瞬变“波形衰减”和“相位滞后”的来源：摩阻越大，压力波沿程衰减越快，远端传感器看到的泄漏特征越弱，这也是小泄漏难以被远端识别的物理原因。

状态方程与简化形式

状态方程建立压力与密度（或水头与压力）之间的关系。对液体，常用线性化形式：

p = p₀ + K·(ρ - ρ₀)/ρ₀
或写成水头形式： H = p/(ρg) + z

把状态方程代入连续性方程，就得到工程上最常用的“水锤方程”组合——以 H、Q 为未知量的两个一阶偏微分方程。进一步做两类常见简化：一是刚性水柱模型，忽略流体与管壁弹性（令 a→∞），只保留惯性项，适用于慢关阀、长周期工况的粗略估算；二是稳态摩阻近似，把非恒定摩阻用瞬时流量代入稳态摩阻公式计算，适用于瞬变较弱、频率较低的场合。对泄漏检测而言，若泄漏引起的压力波动幅度小、频率低，稳态摩阻近似通常够用；但若要分析快速关阀或大口径爆破，必须保留非恒定摩阻与波速的准确取值，否则波形幅值和到达时间都会失真。

一维瞬变流三方程与特征线离散示意 连续性方程（质量守恒） ∂p/∂t + v·∂p/∂x + ρa²·∂v/∂x = 0 波速 a 由 K、ρ、管壁弹性决定 动量方程（力平衡） ∂H/∂x + (1/gA)·∂Q/∂t + fQ|Q|/(2gDA²)=0 摩阻项决定波形衰减与相位滞后 状态方程（压缩性） p = p₀ + K·(ρ−ρ₀)/ρ₀ 把密度变化换算为压力/水头 代入状态方程 → 水锤方程（以 H、Q 为未知量） 简化：刚性水柱（a→∞，忽略弹性）｜稳态摩阻近似（瞬变弱、频率低时适用） 特征线法（MOC）离散：沿 C⁺、C⁻ 特征线求交点 x t C⁺：dx/dt = +a C⁻：dx/dt = −a P 点（t+Δt, x） 上游边界 下游边界 离散求解（Δt = Δx/a） C⁺：H_P = H_A − B(Q_P − Q_A) − R·Q_A|Q_A| C⁻：H_P = H_B + B(Q_P − Q_B) + R·Q_B|Q_B| B = a/(gA)，R = f·Δx/(2gDA²)

特征线法（MOC）的离散思路

上述偏微分方程组无法直接解析求解，工程上最成熟的方法是特征线法。其核心思想是：把两个偏微分方程沿两条特征线方向组合成两个常微分方程，再在 x–t 平面上按网格离散求解。取时间步长 Δt 与空间步长 Δx 满足 Δx = a·Δt（库朗条件），则从上一时刻相邻节点 A、B 出发的两条特征线恰好交于下一时刻的 P 点，形成规则的矩形网格。

沿 C⁺ 和 C⁻ 特征线分别得到两个代数方程，联立即可解出 P 点的 HP 与 QP：

C⁺:  H_P = H_A − B·(Q_P − Q_A) − R·Q_A·|Q_A|
C⁻:  H_P = H_B + B·(Q_P − Q_B) + R·Q_B·|Q_B|
其中 B = a/(gA)，R = f·Δx/(2gDA²)

联立得：Q_P = [H_A − H_B + B(Q_A + Q_B) − R(Q_A|Q_A| − Q_B|Q_B|)] / (2B)
再回代任一式求 H_P。

实际编程时，先按管段长度和波速确定 Δx，再由库朗条件反推 Δt；若管段长度不能被 Δx 整除，需做波速调整或插值处理。每个内节点按上式逐时步推进，计算量随节点数与时间步数线性增长，这也是 MOC 适合长输管道实时泄漏监测的原因。

边界条件处理

内节点有两个方程可解两个未知量，但管道两端各只有一个特征线可用，必须补充边界条件才能闭合。常见处理方式如下：

  • 上游为定压储罐/泵出口：给定 H 或泵扬程曲线，用 C⁺ 方程反求 QP。
  • 下游为定压罐或阀门：给定 H 用 C⁻ 方程求 QP；若为阀门，则联立阀门流量特性 Q = Cv·√(ΔH) 迭代求解。
  • 泄漏点：作为内部边界，把泄漏流量 QL 从上游侧流量中扣除，两侧压力相等，联立 C⁺、C⁻ 与孔口出流公式 QL = CdAL√(2gHL) 求解。
  • 关阀/启泵：用随时间变化的开度或转速曲线作为边界输入，逐时步更新。
建模时务必先统一符号约定：本文取稳态流向为 x 正方向，v、Q 以该方向为正，摩阻项用 |Q|·Q 保证反向流时符号正确。若改用“由下游向上游为正”，则所有流速、流量及特征线方向需同步翻转，否则泄漏定位会出现系统性偏差。

以一段 20 km、内径 0.4 m 的成品油管道为例（示例参数）：取 a = 1050 m/s，Δx = 500 m，则 Δt ≈ 0.476 s，共 40 个计算段。稳态流量 300 m³/h、上游水头 120 m。若在距上游 8 km 处发生当量孔径 20 mm 的泄漏，按上述 MOC 推进可观察到：泄漏点压力先出现一个负压突降，随后以 a 的速度向两端传播，约 7.6 s 到达上游、11.4 s 到达下游。两端压力波到达时间差 Δtd 与泄漏位置满足 xL = (L + a·Δtd)/2，这正是负压波定位的基本公式。需要说明，上述数值为教学示例的合理假设，实际工程中波速、摩阻系数和泄漏当量孔径均需通过现场标定与历史数据反演确定。

掌握三方程与 MOC 离散后，泄漏检测的下一步工作就有了明确的物理基础：波速决定定位精度，摩阻决定可检测的最小泄漏量，边界条件决定模型能否真实复现泵站与阀室的调节行为。后续的传感器布点、采样率选择与定位算法，都应回到这三个方程和特征线网格上来校核。

← 上一章 📖 返回目录 下一章 →