管道明-满流及混合流广泛存在于输水工程[1-2]、输气工程[3]、岩溶管道及暗河[4-6]中。一般情况下，管道过水断面远小于输运长度，水流横向运动在工程尺度可以忽略，因此管道流通常被视作一维问题，并基于此建立了一维管道流控制方程——Saint-Venant方程[7-8]。Saint-Venant方程考虑了管道流动过程中完整动量过程，因此被称为FDM(full dynamic model)[9]。FDM数值精度高，但求解复杂，对空间离散和时间步长要求严格，面对长期水文过程预测问题，计算成本很大，而实际管道流动过程中很多工况下动量变化稳定，所以学者们基于特定问题对Saint-Venant方程的动量方程进行简化，分别提出了QDM (quasi-steady dynamic wave model)、DWM (diffusion wave model)、KWM(kinematic wave model)[9-10]。以上模型都是FDM的简化结果，虽然牺牲了部分数值精度，但大幅降低了对于网格划分和时间步长的要求，在满足工程精度要求的前提下大大提升了计算效率，故简化模型在工程问题中应用前景广泛，如DWM和KWM应用于泥石流运移计算[11-13]，DWM应用于岩溶区管道流动模拟等[14-16]。

Saint-Venant方程是描述管道或渠道明流流动的方程，DWM作为该模型的近似，也难以直接应用于管道满流流动模拟。因此，采用DWM的研究均无法准确描述岩溶管道满流过程。Basha等[17-18]将管道与岩溶区地层渗流耦合，分别建立了忽略流体密度变化的明流模型和满流模型，无法模拟明-满流交替及其混合流动过程；Murad等[19]考虑岩溶区基质-裂隙-管道结构，构建了忽略管道满流形式的岩溶多尺度流动模型，假设岩溶管道都处于部分充水的明流状态。但实际自然水文过程和工程中，管道明、满流交替的混合流动形式是非常普遍的现象。为将DWM应用于管道混合流，针对岩溶区岩溶管道流动，Cornaton等[20]考虑满流时流体压缩性以及管壁变形，提出了满流流动管道储水系数模型；De Rooij等[21]分别考虑管道明流和满流不同储水能力，提出了混合流对流扩散模型(以下简称为Rob模型)，用于预测岩溶区管道流动过程; Malenica等[22]采用Rob模型构建岩溶区管道-基质水动力过程，并与地质模型试验结果对比，验证了模型的适用性。虽然Rob模型考虑了管道明-满流变化过程，但研究表明，该模型仍存在数值计算困难[21]、明流-满流界定不明确[21-22]等问题。此外，前人的研究大多忽略流态变化，假设管道恒为湍流，极少考虑层流、过渡流、湍流的交替演化过程。

本文假设明-满流及混合流流动过程中流体密度连续变化，统一明-满流动过程，采用扩散波动近似简化动量方程，并考虑管道流动的多流态过程，构建了变密度多流态管道流混合扩散对流模型(以下简称为多流态模型)。模型采用Galerkin有限元法和L方法[23]实现数值求解，有效地解决了数值计算困难的问题，实现了不同流态变化管道流水动力过程求解。

采用Leon所构建的单相管道流流动模型，Saint-Venant方程和水击方程可统一表达为[24-25]：

其中：ρα为流体密度(明流时密度为常值，满流时密度改变)，kg/m3；下标α表示不同的流动形式(明流或满流)；A表示过流横截面积，m2；Q表示流量，m3/s；p为过流横截面的流体压力, Pa；g为重力加速度，m/s2；s0为管道坡度，sf为水力梯度, 均为无因次数。

扩散波动近似假设局部对流加速度和时变加速度的影响在式(2)可以忽略，静水压力、重力以及沿程摩擦阻力达到平衡，则动量方程可简化为[21]

其中：h为压力水头或水深，m；按照静水压力计算满足p = ραgh。

此外，基于管道流试验，学者们建立了不同流态条件下关联Q、流速u和sf的经验与半经验模型，从而构建Q与水头的表达式：

其中：H表示总水头(H=h+z，z表示位置水头)，m；K表示体积传导系数，m3/s；可表达为K=φ(h, sf, ge, f)，是h、sf、几何形状ge(迂曲度、粗糙度等)以及流体特性f(运动黏度等)的函数。目前应用广泛的K的计算模型包括基于渠道湍流的Manning-Strickler模型[26-27]和圆管湍流的Darcy-Weibach模型[28-29]，以及层流条件下Poiseuille理论模型[30]。

将式(4)代入式(1)，则可以得到管道混合流扩散波动模型：

其中：ρ为流体密度，会随着流动状态状改变，kg/m3。式(5)考虑流体可压缩性，但实际数值计算或工程应用中，一般认为是水为不可压缩流体，则式(5)可简化为

式(6)即为DWM的最初形式。管道内发生满流流动时，流体密度发生改变且不可忽略，同时横截面积不变，因此式(6)只适用于管道明流流动，无法描述满流流动和混合流流动。

为将DWM推广用于满流及混合流流动问题，必须考虑流体可压缩性。本文认为在明-满流及其混合流流动过程中，流体密度连续可变，故引入流体压缩系数κw(单位为Pa-1) [31-32]：

其中ρ为压强, Pa。

则ρ与p的关系可表达为

其中：ρ0为流体参考密度，kg/m3；p由流体的压力水头h产生，且位置水头即管道高程不发生改变，则：

结合式(8)和(9)，进一步得到ρ与H的关系：

则ρ的全微分表达式为

将式(11)代入式(5)，展开微分项，可得到考虑密度变化的管道流控制方程：

其中：$W=\frac{\partial A}{\partial h}$表示液面顶部宽度，m。对于弱不可压缩流体，其流体压缩系数很小(如水在5 atm(1 atm=101.325 kPa)压强作用下的压缩系数仅为0.538×10-9 Pa-1)，则当管道水深变化较小时，$\rho_0 \kappa_{\mathrm{w}} g \Delta h \ll 1$，则式(12)可以简化为

其中：W代表横截面几何形状的储水能力，ρ0κwgA代表由流体密度和几何形状相互作用的储水能力。可定义C=ρ0κwgA+W，表示管道混合流动过程中储水系数，m。显然明流时管道断面变化是储水能力的主导，而满流时管道断面恒定，密度改变是储水能力的主导。式(13)即本文所构建的基于管道混合流扩散对流模型，通过构建考虑多流态K的计算模型来表征流动过程中的流态变化特征具体见1.2节，以下简称为多流态模型。

基于管道流实验以及河流、输水工程等工程测量数据，学者们建立了管道流流态、流速以及水头损失的经验与半经验模型，表 1列出了管道流常用的经验公式及其适用范围，其中：d为管径，m；ε为粗糙度，m；Rh为水力半径，m；n为Manning系数，m-1/3s。但实际管道流体运移过程中，源汇分布以及管径等结构的变化，会导致流态交替转化，存在层流、过渡流以及湍流演化过程。一般选择Re作为流态判断指标：

式中：D表示水力直径，m；ν表示流体运动黏度，m2/s。工程上认为Re < 2 000为层流，2 000≤Re≤4 000为过渡流，Re>4 000为湍流。

表 1中模型只针对层流或湍流状态，缺乏对过渡流的描述，也无法涵盖流态演化过程，导致计算精度降低。为准确描述管道流多流态特征，本文采用Swamee-Swamee模型[33]，该模型是Swamee基于管道流全流态实验数据构建的全流态经验公式：

其中：hf为管道水头损失，m；L为管道长度，m；D为管道水力直径，m；

结合式(4)，可得到K表达式为

需要注意的是，Swamee-Swamee公式是基于管道满流所构建的模型，但研究表明当管道明流(部分充满)时，只要D仍表示该状态下的水力直径，仍可保证方程的适用性[34]。图 1比较了相同条件下不同经验模型计算结果，其中d=0.2 m，ε=0.000 1 m，n=0.012 s·m-1/3, ν=10-6m2/s，并根据Re确定了流态。Re下、上限分别采用Poiseuille公式和Manning-Strickler公式计算。显然，无论满流或明流下，水力梯度的增大都会导致流动逐渐向湍流发展，且明流时管道内水深增加(sf恒为0.000 1) 也会使流态逐渐演化为湍流。当水深较小时，流动仍处于层流状态，因此管道混合流过程中必然存在层流-过渡流-湍流的变化过程。

Swamee-Swamee公式可以准确地刻画管道流动全流态演化过程：在Re < 2 000时，与Poiseuille层流模型一致；在Re>4 000时，与Darcy-Weibach湍流模型结果一致；而在2 000≤Re≤4 000时，实现了层流到湍流的光滑过渡，描述了过渡流的水力特征。本文采用该模型来表征管道流动的流态演化过程，即式(13)中的K通过式(16)计算得到。

Cornaton等[20]在构建岩溶管道模型时，同时考虑流体压缩和管壁弹性形变提出了满流状况下管道储水系数表达式：

其中：γ表示流体容重，N/m3；Ew表示流体体积模量，GPa；Ec表示管壁体积模量，GPa；rc表示管道半径，m。不考虑管壁变形时，式(17)等价于

其中Af表示满流管道截面积，m2。

在此基础上，De Rooij[21]修正扩散波动模型的管道储水系数，分别给出满流和明流时的储水系数，构建管道混合流扩散波动模型即Rob模型：

Rob模型忽略明流和满流的统一变化过程，在明流和满流过渡状态时，C值发生突变，导致数值求解困难[21]。为保证数值稳定性，Rob模型将接近满流的明流作为满流处理，导致明流和满流的界定模糊不明确。图 2比较了Rob模型和多流态模型C值随h/d的变化规律，显而易见：Rob模型的C值是不连续的，在明流和满流交界处(图 2中圆圈位置)发生间断跳跃，导致数值计算困难，且在明流阶段并未考虑密度变化影响；多流态模型的C值是连续的，在明流阶段，同时考虑了截面几何与流体密度变化对储水系数的影响。

相比Rob模型对明流和满流的界限模糊处理，多流态模型明确了满流和明流的界定：