管道各个截面的液流平均速度和压力不随时间变化或变化很小的流动称为稳定流动，反之称为非稳定流动。管道中工作元件工作状态的变化或外界干扰，都可能使液体管道的流动变为不稳定，例如开关阀门、启停泵机组和发生泄漏事故等。管道内的不稳定流动虽属于从一个稳定工况到另一个稳定工况的过渡工况，但其引发的问题不容忽视。由于不稳定流动产生的过高或过低的压力以及管道、设备的振动会危及管道安全，因此，研究传输、传动和控制工程中管道的不稳定流动已逐渐形成一门独立的学科分支———流体传输管道动力学。

一.   管道内非稳定流动及其分析方法

研究液体管道不稳定流动的方法主要有频域分析法和时域分析法^{〔1, 2〕}。频域分析法是指在频率域内对管内流体的频率特性进行分析。流体管道属于分布参数系统，可用传输线理论进行分析。由于频率特性分析中主要针对管道和系统元件的阻抗特性进行分析，因此又将频域分析称为阻抗特性分析。

对管内非稳定流动的时域特性分析，一般称为管道瞬态特性分析，它主要是在不同的边界条件下，对管道的瞬态响应进行时域仿真计算，求得系统的动态特性。描述液体管道内一维非稳定流动的基本方程为运动方程和连续性方程^〔3〕，即

式中  H———管内液体压头，m；

g———重力加速度，m/s²；

A———管道截面积，m²；

Q———流量，m³/s；

V———管内液体流速，m/s；

h_u———非稳定流动时因管道内液体与管壁的摩擦阻力引起的压头损失，m；

a———波速，m/s。

式（1）中压头损失项h_u直接影响管道中各点压力和流量的计算，如何计算其值是管道非稳定流动分析领域中的重要问题。

二.   传统的瞬态摩阻计算方法

稳定流动时不可压缩单相流体沿程阻力一般利用达西公式计算：

式中  h_s———稳定流动时的压头损失，m；

λ———达西摩阻系数；

l———管道长度，m；

d———管道内直径，m。

达西摩阻系数是液体性质、流动状态和管子粗糙程度的函数。对于层流流动：

式中  Re———雷诺数；

v———液体的运动粘度。

对于紊流流动下λ的计算，各国学者提出的公式不下百个。通常的作法是将紊流流动区域再划分为水力光滑区、混合摩擦区和阻力平方区，分别给出各流动区域的计算公式。文献[1]指出，这种分区是人为的分区，根据Moody图，λ值从水力光滑区到混合摩擦区是一个渐变的过程，由分区引起的λ值的突变在理论和实践上都无法解释^〔4〕。柯尔布鲁克公式能够很好地符合Moody图，用其计算整个紊流区的λ值是切实可行的。柯尔布鲁克公式为：

式中  ε———相对粗糙度；

k———等值绝对粗糙度，m。

对于非稳定流动下的摩擦阻力及相应的压头损失的计算，通常的作法是进行所谓的“拟稳态”假设，即在管道某一截面上平均流速为V的瞬时摩阻，等于该截面在同一流速V下稳定流动时的摩阻。这种假设也用来模拟阀、水力机械、孔口等非管元素。

不少学者通过外加的可改变频率的振动来验证这种“拟稳态”假设^〔2〕。Safwat和Polder（1973）作了层流流动下的外加振动试验，发现管壁的剪切应力不仅是雷诺数的函数，而且与振动的频率有关。Latelier和Leutheusser（1976）也得出结论，层流流动的起始区域和振动区域的摩擦阻力不符合“拟稳态”假设。Ramaprian和Tu（1980）作了层流区、过渡区和紊流光滑区的外加振动试验，得出了紊流状态下的平均能量耗散大于按“拟稳态”假设得到的值的结论。

如果只考虑非稳定流动引起的瞬时压力的极值而接受“拟稳态”假设，以达西公式计算管道非稳定流动下的摩阻损失已经足够了。然而，对于长时间的动态模拟以达到新的稳定状态而言，考虑非稳定流动的影响，更准确地计算摩阻损失则是十分必要的。

三.   非稳定流动下摩阻的计算方法

如前所述，粘性流体在管道内非稳定流动时，管内流体与管壁间的摩阻力损失同流体脉动频率大小有关，具有频率相关的特性，称为频率相关摩阻损失。为了更准确地计算非稳定流动下的摩阻，必须引入频率相关摩阻模型。

在频域分析领域中，Zielke（1968）从N-S方程出发，推导出层流流动下的单位长度管壁压头损失的计算公式：

式中  ρ———流体密度，kg/m³；

r———管子半径，m；

W（t-t₁）———权函数。

式（6）中等式右边第一项对应于层流流动的稳态压头损失，第二项对应于频率相关压头损失。Zielke公式被公认为频率相关摩擦损失的精确模型。

引入无量纲时间因子和权函数W（t）=W（τ），并取：

将式（7）代入式（6）即可得到精度很高的频率相关摩阻的近似。

Zielke公式的重要性主要还在于它发展了频率相关摩擦损失在紊流流动区域及在时域的计算方法。Trikha（1975）提出一种近似的权函数W_app（τ）去逼近W（τ），并假设3个条件，即

（1）；

（2）W_app（r）=W（r），（当r=0.000 1，0.001，0.01时）；

（3），（当r>0.01时）。

根据上述条件可得到近似的权函数为：

紊流在流动时仍取频率相关压头损失项为

将式（8）代入式（6），由此得到的近似模型以一阶惯性量组成，因此，可以方便地通过拉普拉斯（Laplace）逆变换到时域，式（6）可变换为：

由此可推出：

式（10）即为时域中频率相关压头损失的差分方程，其中m_i=40，8.1，1；n_i=8 000，200，26.4（i=1，2，3）。

四.   包含频率相关摩阻的非稳定流动微分方程求解

式（1）、式（2）合称波动方程，常见的求解方法有特征线法、直接差分法、有限单元法等。由于特征线法物理意义明确，并可非常方便地用计算机进行快速运算，其优越性已大大超过从前用于管道非稳定流动分析的解析法和图解法，因此得到越来越多的工程应用，由于式（1）、式（2）可构成一阶拟线性偏微分方程组，因此存在两条特征线，沿两条特征线可将原方程组化为常微分方程组。

式中  C_w———惯性水击常数。

必须注意，式（11）只有沿特征线式（12）才成立。

对于液体管道，有V≪a，则。上述方程组可以通过构造差分格式，用有限差分的方法求解。

如图 1所示，将长为l的管子等分为n段，每段长度为Δx=l/n，纵轴为时间轴，且令每个单元的长度为Δt=Δx/a，称为时间步长。在如此构成的差分网格中，对角线AP满足Δx=aΔt，而对角线BP则满足Δx=-aΔt，称为特征网格。

图  1  特征网格

下载: 全尺寸图片 幻灯片

显然，若已知t₀时刻A、B的流量和压头Q_A、H_A、Q_B和H_B，则将式（11）分别沿特征线C⁺：Δx=aΔt和C^-：Δx=-aΔt积分，可得：

由于h_u是流量的函数，式（13）和式（14）中的项不能简单求积，可用梯形求积公式计算：

将式（15）和式（16）分别代入式（13）和式（14）并作整理，可得：

C_L、C_R称为左右特征项，则：

在式（22）中，左右特征项C_L、C_R均为Q_P的函数，可表示为Q_P=Φ（Q_P）的形式，用埃特金迭代公式求解Q_P。取Q_P，0=（Q_A+Q_B）/2作为迭代初值，迭代格式为：

式（23）中的k=0，1，2，…，为迭代次数。显然，当管内液体紊流流动时，因为柯尔布鲁克公式为隐式，所以达西摩阻系数λ也必须进行迭代计算。