管道明-满流及混合流广泛存在于输水工程^[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]实现数值求解，有效地解决了数值计算困难的问题，实现了不同流态变化管道流水动力过程求解。

1 模型与方法 1.1 可压缩流体扩散波动模型

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

$ \frac{\partial\left(\rho_\alpha A\right)}{\partial t}+\frac{\partial\left(\rho_\alpha Q\right)}{\partial x}=0, $

(1)

$ \frac{\partial\left(\rho_\alpha Q\right)}{\partial t}+\frac{\partial\left(\rho_\alpha \frac{Q^2}{A}+A \bar{p}\right)}{\partial x}=\rho_\alpha\left(s_0-s_f\right) g A. $

(2)

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

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

$ \frac{\partial h}{\partial x}=s_0-s_{\mathrm{f}} . $

(3)

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

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

$ Q=-K \frac{\partial H}{\partial x}. $

(4)

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

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

$ \frac{\partial(\rho A)}{\partial t}-\frac{\partial\left(\rho K \frac{\partial H}{\partial x}\right)}{\partial x}=0. $

(5)

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

$ \frac{\partial A}{\partial t}-\frac{\partial\left(K \frac{\partial H}{\partial x}\right)}{\partial x}=0. $

(6)

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

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

$ \kappa_{\mathrm{w}}=\frac{\mathrm{d} \rho}{\rho \mathrm{d} p}. $

(7)

其中ρ为压强, Pa。

则ρ与p的关系可表达为

$ \rho=\rho_0+\rho \kappa_{\mathrm{w}} \Delta p \text {. } $

(8)

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

$ \Delta p=\rho g \Delta h=\rho g \Delta H. $

(9)

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

$ \rho=\frac{\rho_0}{1-\rho_0 g \kappa_{\mathrm{w}} \Delta H} . $

(10)

则ρ的全微分表达式为

$ \mathrm{d} \rho=\frac{\rho_0^2 g \kappa_{\mathrm{w}}}{\left(1-\rho_0 g \kappa_{\mathrm{w}} \Delta H\right)^2} \frac{\partial H}{\partial t} \mathrm{~d} t+\\ \frac{\rho_0^2 g \kappa_{\mathrm{w}}}{\left(1-\rho_0 g \kappa_{\mathrm{w}} \Delta H\right)^2} \frac{\partial H}{\partial x} \mathrm{~d} x . $

(11)

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

$ \begin{gathered} \left(\frac{\rho_0 \kappa_{\mathrm{w}} g A}{1-\rho_0 \kappa_{\mathrm{w}} g \Delta H}+W\right) \frac{\partial H}{\partial t}- \\ \frac{\rho_0 \kappa_{\mathrm{w}} g K}{1-\rho_0 \kappa_{\mathrm{w}} g \Delta H}\left(\frac{\partial H}{\partial x}\right)^2-\frac{\partial\left(K \frac{\partial H}{\partial x}\right)}{\partial x}=0. \end{gathered} $

(12)

其中：$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)可以简化为

$ \begin{aligned} & \left(\rho_0 \kappa_{\mathrm{w}} g A+W\right) \frac{\partial H}{\partial t}-\rho_0 \kappa_{\mathrm{w}} g Q \frac{\partial H}{\partial x}- \\ &\;\;\;\;\;\;\;\;\; \frac{\partial\left(K \frac{\partial H}{\partial x}\right)}{\partial x}=0 . \\ & \end{aligned} $

(13)

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

1.2 多流态模型

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

$ R e=\frac{4 Q}{\nu \pi D}. $

(14)

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

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

$ \begin{aligned} & Q=D^2 \sqrt{g D h_{\mathrm{f}} / L}\left\{\left(\frac{128 \nu}{\pi D \sqrt{g D h_{\mathrm{f}} / L}}\right)^4+\right. \\ & 1.153\left[\left(\frac{415 \nu}{D \sqrt{g D h_{\mathrm{f}} / L}}\right)^8-\right. \\ & \left.\left.\ln \left(\frac{\varepsilon}{3.7 D}+\frac{1.775 \nu}{D \sqrt{g D h_{\mathrm{f}} / L}}\right)\right]^{-4}\right\}^{-0.25} . \\ & \end{aligned} $

(15)

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

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

$ \begin{aligned} & K= D^2 \sqrt{\frac{g D L}{\left|h_{\mathrm{f}}\right|}}\left\{\left(\frac{128 \nu}{\pi D \sqrt{g D\left|h_{\mathrm{f}} / L\right|}}\right)^4+\right. \\ & 1.153\left[\left(\frac{415 \nu}{D \sqrt{g D\left|h_{\mathrm{f}} / L\right|}}\right)^8-\right. \\ &\left.\left.\ln \left(\frac{\varepsilon}{3.7 D}+\frac{1.775 \nu}{D \sqrt{g D\left|h_{\mathrm{f}} / L\right|}}\right)\right]^{-4}\right\}^{-0.25}. \end{aligned} $

(16)

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

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

1.3 管道流模型对比

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

$ C=\gamma\left(\frac{1}{E_{\mathrm{w}}}+\frac{1}{E_{\mathrm{c}}}\right) \pi r_{\mathrm{c}}^2 . $

(17)

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

$ C=\frac{1}{E_{\mathrm{w}}} \rho g \pi r_{\mathrm{c}}^2=\rho g \kappa_{\mathrm{w}} A_{\mathrm{f}}. $

(18)

其中A_f表示满流管道截面积，m²。

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

$ \begin{gathered} C \frac{\partial H}{\partial t}-\frac{\partial\left(K \frac{\partial H}{\partial x}\right)}{\partial x}=0 ; \\ C= \begin{cases}\rho_0 \kappa_{\mathrm{w}} g A_{\mathrm{f}}, & \text { 满流；} \\ W, & \text { 明流. }\end{cases} \end{gathered} $

(19)

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

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

$ \begin{aligned} & A= \begin{cases}\frac{1}{4} d^2 \arccos \left(1-\frac{2 h}{d}\right)+ & \\ h\left(h-\frac{d}{2}\right) \sqrt{\frac{d}{h}-1}, & h <d ; \\ \frac{1}{4} \pi d^2, & h \geqslant d ;\end{cases} \\ & D=\left\{\begin{array}{l} \frac{d^2 \arccos \left(1-\frac{2 h}{d}\right)+h\left(h-\frac{d}{2}\right) \sqrt{\frac{d}{h}-1}}{d \arccos \left(1-\frac{2 h}{d}\right)}, \\ h<d ; \\ d, \quad h \geqslant d . \end{array}\right. \\ & \end{aligned} $

(20)

其中h < d时为明流，h≥d时为满流。

此外，Rob模型仅考虑明流和满流时管道储水系数的差异，未考虑满流时密度变化所引起的质量对流，质量方程不完全守恒，而多流态模型通过增加式(13)等号左侧第2项来修正该误差。因此，在采用Rob模型计算时，会出现奇异值和较大的偏差，而多流态模型取得更好的计算结果。

1.4 数值算法

在由空间域Ω和时间域[t, t_end]共同构建的计算域上，求解式(13)和(16)的数值解H的时空间分布，其初值和边界条件已知，可表示为

$ \begin{cases}H=H_0, & \text { on } \quad {\mathit{\Omega}}, \quad t=0 ; \\ H=H_{\mathrm{D}}, & \text { on } {\mathit{\Gamma}}_{\mathrm{D}}, \quad t \in\left[0, t_{\mathrm{end}}\right] ; \\ -K \frac{\partial H}{\partial x} \cdot \vec{n}=q_{\mathrm{N}}, & \text { on } {\mathit{\Gamma}}_{\mathrm{N}}, \quad t \in\left[0, t_{\mathrm{end}}\right] .\end{cases} $

(21)

其中：H₀为初始时刻的水头值，Γ_D为Dirichlet边界(水头边界)，H_D为Γ_D上的水头值，Γ_N为Neumann(流量)边界，q_N为Γ_N上的流量值，$\vec{n}$为Γ_N的外法向量，t_end表示计算终止时刻。

假设方程真解为H(x, t)，采用Galerkin有限元法对Ω进行空间离散形成单元集合，则该问题的数值解可表达为

$ H_{\mathit{\Omega}}(x, t)=\sum\limits_{j=1}^n H_j(t) \varphi_j(x) . $

(22)

其中: H_j表示Ω上节点j的值，φ_j表示Ω上节点j的基函数，n表示Ω上的节点数目。将式(21)—(22)代入式(13)，利用Einstein求和公式表达，可得到数值弱解表达形式为

$ \begin{array}{l} \int_\mathit{\Omega } {\left( {{\rho _0}{\kappa _{\rm{w}}}gA + W} \right)} \frac{{\partial {H_j}}}{{\partial t}}{\varphi _j}{\varphi _k}{\rm{d}}\mathit{\Omega } - \\ \;\;\;\;\;\int_\mathit{\Omega } {{\rho _0}} {\kappa _{\rm{w}}}gQ\frac{{\partial {H_j}}}{{\partial x}}{\varphi _j}{\varphi _k}{\rm{d}}\mathit{\Omega } - \\ \int_\mathit{\Omega } {\frac{{\partial \left( {K\frac{{\partial {H_j}}}{{\partial x}}{\varphi _j}} \right)}}{{\partial x}}} {\varphi _k}{\rm{d}}\mathit{\Omega } + \int_\mathit{\Omega } {{q_m}} {\varphi _k}{\rm{d}}\mathit{\Omega } = 0. \end{array} $

(23)

利用Gauss公式，式(23)可整理得到

$ \begin{gathered} \int_{\mathit{\Omega }}\left(\rho_0 \kappa_{\mathrm{w}} g A+W\right) \frac{\partial H_j}{\partial t} \varphi_j \varphi_k \mathrm{~d} \mathit{\Omega }- \\ \int_{\mathit{\Omega }} \rho_0 \kappa_{\mathrm{w}} g Q \frac{\partial H_j}{\partial x} \varphi_j \varphi_k \mathrm{~d} \mathit{\Omega }+ \\ \int_{\mathit{\Omega }} K H_j \frac{\partial \varphi_j}{\partial x} \frac{\partial \varphi_k}{\partial x} \mathrm{~d} \mathit{\Omega }+\int_{\mathit{\Omega }} q_m \varphi_k \mathrm{~d} \mathit{\Omega }+\int_{\partial \mathit{\Omega }} q_{\mathrm{N}} \varphi_k \mathrm{~d} \mathit{\Omega }=0 . \end{gathered} $

(24)

采用一阶隐式Euler法进行时间离散，则式(24)可表示为

$ \begin{gathered} \int_{\mathit{\Omega }}\left(\rho_0 \kappa_{\mathrm{w}} g A+W\right) \varphi_j \varphi_k \frac{H_j^{t+\Delta t}-H_j^t}{\Delta t} \mathrm{~d} \mathit{\Omega }- \\ \int_{\mathit{\Omega }} \rho_0 \kappa_{\mathrm{w}} g Q \frac{\partial \varphi_j}{\partial x} \varphi_k H_j^{t+\Delta t} \mathrm{~d} \mathit{\Omega }+ \\ \int_{\mathit{\Omega }} K \frac{\partial \varphi_j}{\partial x} \frac{\partial \varphi_k}{\partial x} H_j^{t+\Delta t} \mathrm{~d} \mathit{\Omega }+ \\ \int_{\mathit{\Omega }} q_m \varphi_k \mathrm{~d} \mathit{\Omega }+\int_{\partial \mathit{\Omega }} q_{\mathrm{N}} \varphi_k \mathrm{~d} \mathit{\Omega }=0 . \end{gathered} $

(25)

其中：t和Δt分别表示时间和时间步长，可知式(25) 等号左侧各系数项都与要求解变量H^t+Δt有关，是非线性对流-扩散方程，一般采用非线性迭代方法求解，如Newton迭代法等。前人研究表明该类方程存在强烈的非线性特征，Newton迭代法不能保证计算的稳定性和收敛性，且管道混合流流动过程中，明流和满流流动机理迥异，数值突变显著，采用Newton迭代法不容易数值收敛，因此学者们提出了Picard方法、修正Picard方法以及L方法等改善算法收敛性^[35-36]。经过测试，本文采用L方法，在满足该方法线性化参数L_n≥d条件下，保证了算法的无条件收敛。

采用L方法求解多流态模型，首先定义：

$ \mu=\int_0^H C \mathrm{~d} H=\int_0^H\left(\rho_0 \kappa_{\mathrm{w}} g A+W\right) \mathrm{d} H . $

其中μ表示单位长度管道储水量，由流体密度和过水断面面积共同决定，m²。对于管径为d的圆形管道，μ可以表达为

$ \mu=\left\{\begin{array}{lc} {\left[r^2+\rho_0 g \kappa_{\mathrm{w}} r^2(h-r)\right] \arccos \left(1-\frac{h}{r}\right)+} \\ {\left[\rho_0 g \kappa_{\mathrm{w}} r^2+\frac{1}{3} \rho_0 g \kappa_{\mathrm{w}} h(2 r-h)+(h-r)\right] \cdot} \\ \sqrt{h(h-r)}, & h \geqslant d ; \\ \pi r^2\left(1+\rho_0 g \kappa_{\mathrm{w}}(h-r)\right), & h <d . \end{array}\right. $

(26)

其中r为管道半径，m。

图 3展示了混合流动过程中，Rob模型和多流态模型μ随管道水深比h/d的变化趋势，显然，μ的函数在明流和满流的临界处(图 3圆圈位置)连续可导。

因此，经过Galerkin有限元法离散和L方法求解，方程数值解可以表示为

$ \begin{aligned} & \left(C_{t+\Delta t}^{n-1}-\Delta t \times V_{t+\Delta t}^{n-1}+\Delta t \times K_{t+\Delta t}^{n-1}\right) H_{t+\Delta t}^n= \\ & C_{t+\Delta t}^{n-1} \times H_{t+\Delta t}^{n-1}-Q_m-Q_N+M_t-M_{t+\Delta t}^{n-1}. \end{aligned} $

(27)

其中：

$ \begin{gathered} C_{t+\Delta t}^{n-1}=\int_{\mathit{\Omega }}\left(\rho_0 \kappa_{\mathrm{w}} g A_{t+\Delta t}^{n-1}+W_{t+\Delta t}^{n-1}\right) \varphi_j \varphi_k \mathrm{~d} \mathit{\Omega }, \\ V_{t+\Delta t}^{n-1}=\int_{\mathit{\Omega }} \rho_0 \kappa_{\mathrm{w}} g Q_{t+\Delta t}^{n-1} \frac{\partial \varphi_j}{\partial x} \varphi_k \mathrm{~d} \mathit{\Omega }, \\ K_{t+\Delta t}^{n-1}=\int_{\mathit{\Omega }} K_{t+\Delta t}^{n-1} \frac{\partial \varphi_j}{\partial x} \frac{\partial \varphi_k}{\partial x} \mathrm{~d} \mathit{\Omega }, \\ Q_m=\Delta t \int_{\mathit{\Omega }} q_m \varphi_k \mathrm{~d} \mathit{\Omega } , \\ Q_{\mathrm{N}}=\Delta t \int_{\mathit{\Omega }} q_{\mathrm{N}} \varphi_k \mathrm{~d} \mathit{\Omega }, \\ M_{t+\Delta t}^{n-1}=\int_{\mathit{\Omega }} \mu_{t+\Delta t}^n \varphi_k \mathrm{~d} \mathit{\Omega }, \\ M_t=\int_{\mathit{\Omega }} \mu_t \varphi_k \mathrm{~d} \mathit{\Omega } . \end{gathered} $

其中：n和(n-1)表示(t+Δt)时刻的不同的迭代水平；Q_N表示计算域边界Δt时段内进水量或出水量，该项只存在于Neumann边界单元。

2 结果 2.1 Van Rijn案例

基于DWM，Van Rijn等^[37]针对宽浅渠道给出了上游水位突变条件下稳定状态下沿程水深理论解：

$ h=\left(\frac{13}{3}(L-x) \frac{Q^2}{B^2} n^2+h_L^{\frac{13}{3}}\right)^{\frac{3}{13}} . $

其中：x表示距离渠道上游的距离，m；B表示渠道宽度，m；h_L表示x=L位置处压力水头。该理论解假设流量和水力梯度满足Manning-Strickler公式。

为验证代码的准确性，设计宽浅渠道计算案例: s₀=0，L=500 m，B=20 m，Q=0.1 m³/s，下游恒定水深h_L=0.02 m，初始时刻渠道各处水深为0.02 m, n=0.02 s·m^-1/3。数值求解过程中，设置x步长，Δx=1 m，Δt=10 s，渠道水深随时间变化过程如图 4所示。显然上游持续通水，渠道沿程水深逐渐增加并在约6 h时达到稳定状态，数值计算结果和理论值完全重合，证明了计算代码的准确性。

2.2 Rossman混合流案例

该案例是Rossman为测试SWMM软件极限功能所构建的计算工况^[38]，由5节变管径管道组成，各节长度都为1 000 ft (1 ft=0.304 8 m)，第1、3、5节管道管径为12 ft，第2、4节管道管径为3 ft，s₀为0.05%，n为0.02 m^-1/3 s。管道初始时刻均为干燥状态，在上游段注水，从开始到15 min，注水量线性增加至50 ft³/s，此后注水量保持不变直到3 h，然后在15 min内线性减少为0，计算总时长为6 h。在此过程中，管道会经历明流-满流-明流的演化过程，且管径的变化也会导致同时刻管道沿途出现明满流交替现象。为研究本文所提出多流态模型的适用性及准确性，对该案例进行数值模拟。数值计算过程中，设置Δx=100 ft，Δt=30 s，并比较了多流态模型、Rob模型和SWMM的计算结果, 结果如图 5所示。

SWMM采用FDM进行求解，理论上计算结果更准确，因此将SWMM结果作为基准。图 5表明，相比Rob模型，多流态模型与SWMM计算结果更加接近，尤其是水深计算，相对误差更小。Rob模型虽然整体上和SWMM结果也相近，但在满流流动状态下局部位置出现奇异值，水深和流量发生突变。因为在明流逐渐发展成为满流过程中，流体密度发生改变且其作用不可忽略，会产生密度波，其波速一般比明流流动产生的重力波波速高2个数量级，是满流产生水击现象或水击压力的本质原因。因此临近满流后，必须考虑流体密度变化所带来的质量守恒问题，但Rob模型并未考虑，故在计算临界流动和满流流动时会发生水深或流量的突变，而多流态模型修正了密度变化所带来质量对流作用，即式(13)左侧第2项，计算结果更准确合理。

SWMM的计算结果在明-满流交替处容易发生数值振荡，而Rob模型和多流态模型的计算结果在此处基本无振荡，这是SWMM采用显式时间差分方法而导致的。此外，SWMM基于FDM模型，计算过程必须满足CFL条件^[39]，即|CFL=vΔt/Δx|≤1，故文[38]中Δx=100 ft时设置Δt=5 s，如果Δt再增大，计算会不收敛。本文所采用多流态模型经过对动量方程的简化，对时间步长限制更少(CFL条件无需满足)，经测试当Δx=100 ft时，设置Δt=150 s，求解仍然收敛并得到稳定的数值解，这表明在动量可简化的条件下，多流态模型非常适合长期水文过程预测问题。

图 5还表明，多流态模型的计算结果虽然与SWMM结果非常接近，但仍然存在一定的差异，主要是由于SWMM采用FDM求解，并未简化动量方程；此外SWMM默认流动过程一直处于湍流流态，采用Manning-Strickler模型，所以在来水和退水过程(实际存在层流-湍流流态演化过程)水深较小时流量计算差异较大，同理Rob也采用Manning-Strickler模型，忽略流态演化过程，因此Rob模型尽管会出现局部位置奇异值，但来水和退水过程和SWMM计算结果更接近。

2.3 Capart管道流实验

Capart等^[39]开展了单相管道流实验，d为0.145 m，总长12.74 m，由3段节管道组成，其管道长度依次为3.48、5.75和3.51 m，坡度依次为1.954%、1.704%和1.225%。实验中控制上游和下游水箱水位，并在管道沿程测定断面测压水头，监测管道内水深变化。本文采用多流态模型对该实验进行数值模拟，并与Capart实验结果和模拟结果进行了比对，结果见图 6，其中C3、C4和C6分别表示距离上游入口3.06、5.50和7.64 m的管道断面位置。

图 6表明，在C3、C4、C6不同断面位置，多流态模型模拟结果与实验结果都非常吻合，相比采用QMD的Capart模拟结果，多流态模型模拟结果与实验结果相对误差更小，在涨水阶段表现尤为明显。上述结果表明在动量变化不剧烈，重力、水压力和摩擦力近似平衡的条件下，采用多流态模型可以更准确刻画管道流的水动力学过程。

3 结论

本文考虑流体可压缩性，假设管道明-满流及混合流流动过程中流体密度连续变化，统一了管道明-满流流动过程，利用扩散波动近似简化了动量方程并考虑管道多流态特征，引入Swamee-Swamee公式表征了管道层流-过渡流-湍流演化过程，通过理论推导分析构建了变密度管道混合流多流态扩散波动模型，并开展了相应的数值实验，得到以下结论：

1) 无论明流或满流，随着水深、流量的变化，管道流动过程中存在层流-过渡流-湍流多流态转化交替过程，层流、过渡流态也不可忽略。

2) 多流态模型的计算结果与SWMM计算结果相近，满足工程需求，且数值稳定性更高，对时间步长要求低，更适合长期水文动力过程预测；计算结果与管道流实验的对比，进一步验证了模型的适用性和准确性。在动量简化条件满足的条件下，多流态模型可以准确描述管道混合流动过程。

3) 与Rob模型计算结果相比，多流态模型计算精度更高。多流态模型考虑了明-满流过程中流体密度变化特征，修正了密度变化所带来质量对流作用，保证了满流流动过程中质量守恒，进一步减少了计算误差，消除了Rob模型计算结果出现的局部奇异值。

4) 与Rob模型计算结果相比，多流态模型收敛性更好。多流态模型假设明-满流流动过程中流体密度连续变化，保证了管道储水系数的连续性，避免了数值突变；多流态模型采用L方法进行数值求解，定义单位长度管道储水量，保证其连续可导性，当L_n≥d时，算法无条件收敛。

在研究过程中也发现，多流态模型仍存在缺陷，利用扩散波动近似对动量方程简化，则对于突变型管道流动(如阀门关闭)等动量突变状况，模型难以适用。为解决该问题，需要考虑管道流动中完整动量变化过程，将连续密度假设引入FDM，建立管道混合流更精确的求解模型，这也是后续需要开展的工作。