给水管网非恒定流动数值计算方法

D01:10.13374.isml00103x.2009.04.019 第31卷第4期 北京科技大学学报 Vol.31 No.4 2009年4月 Journal of University of Science and Technology Beijing Apr.2009 给水管网非恒定流动数值计算方法 马飞1 曲世琳D吴民2》 1)北京科技大学土木与环境工程学院。北京1000832)北京联合大学机电学院，北京100020 摘要为计算给水管网中非恒定流动运行工况使系统可以快速平稳地完成工况转换，以非恒定流动理论为基础分别利 用重分阻尼系数法划分时间步长的特征线法和引入摩阻附加项的玻尔兹曼网格法计算给水管网非恒定流动的数学模型，并 应用于某小区给水管网的非恒定流动工况分析中.结果表明.在管网非恒定流工况下，当R7000时利用玻尔兹曼网格法 计算否则利用特征线法计算将提高计算的准确性和计算效率. 关键词给水管网：非恒定流：特征线法：玻尔兹曼网格法 分类号TV131202418 Numerical calculation methods of water distribution systems in hydraulic transient conditions MA Fe.QU Shi-lin WU Yi-min2) 1)School of Civil and Envimnmental Engineering,University of Science and Technology Beijing.Beijng 100083.China 2)Mechanical and Electrical Engineering School Beijing Union University.Beijng 100020.China ABSTRACT Transient phenomena often happen in water distribution systems(WDS).In order to accelerate system transformation fiom unsteady state to steady state two mathematical models were proposed based on the theory of transient flow.which are the method of characeristics (MOC)with re-divide damping coefficient and the lattice Boltzmann equation (LBE).These methods were used to calculate an example in WDS.It is show n that LBE is more accurate when Rey nolds number is less than 7000,otherwise MOC is more precise KEY WORDS water distribution system;unsteadly flow:method of characteristics lattice Boltzmann equation 由于给水管网是动态的非线性系统.管网的运 1给水管网非恒定流动物理模型 行状态往往是从一个稳定的运行状态过渡到另一个 稳定的运行状态刂.建立在恒定流动理论基础之上 非恒定流动中液体质点的运动要素不仅随空间 的管网水力模型是非常成熟的，在管网的平稳运行 位置变化，而且随时间过程变化？.因此分析这类 工况建模中，这种方法能够精确计算管网各控制点 运动，要考虑由于运动要素随时间变化所引起的惯 的压力、流量等运行参数.然而，在过渡阶段这种方 性力的作用，即其运动要素为时间和空间的函数. 法很难准确地确定管网的运行工况.因此，本文从 由质量和能量守恒原理得管道非恒定流动控制方 非恒定流动的角度出发，分别利用重分阻尼系数法 程内： 划分时间步长的特征线法(method of characteris- 盟+a2Y=0 (1) dt'g dx is,MOC)和引入摩阻附加项的玻尔兹曼网格法 (lattice Boltzmann equation,LBE)求解城市给水管 影+g要+川川=0 (2) 网系统非恒定流动数学模型，进一步完善给水管网 式中，V为管内某断面的平均流速，H为某点测压 动态分析方法， 管水头，D为管径，x为沿流线的位置坐标，【为时 收稿日期：200805-22 基金项目：国家“十一五”科技支撑计划重大资助项目(No.2006BA03A06) 作者简介：马飞(1968一).男.教授，博土，E-mail:ustb巧eke@.ina.com

给水管网非恒定流动数值计算方法 马 飞1) 曲世琳1) 吴一民2) 1)北京科技大学土木与环境工程学院, 北京 100083 2)北京联合大学机电学院, 北京 100020 摘 要 为计算给水管网中非恒定流动运行工况, 使系统可以快速平稳地完成工况转换, 以非恒定流动理论为基础, 分别利 用重分阻尼系数法划分时间步长的特征线法和引入摩阻附加项的玻尔兹曼网格法计算给水管网非恒定流动的数学模型, 并 应用于某小区给水管网的非恒定流动工况分析中.结果表明, 在管网非恒定流工况下, 当 Re≤7 000 时利用玻尔兹曼网格法 计算, 否则利用特征线法计算将提高计算的准确性和计算效率. 关键词 给水管网;非恒定流;特征线法;玻尔兹曼网格法 分类号 TV 131.2;O 241.8 Numerical calculation methods of water distribution systems in hydraulic transient conditions MA Fei 1), QU Shi-lin 1), WU Yi-min 2) 1)School of Civil and Environment al Engineering , Universit y of Science and Technology Beijing , Beijing 100083 , China 2)Mechanical and Electrical Engineering S chool, Beijing Union University , Beijing 100020 , China ABSTRACT Transient phenomena often happen in w ater distributio n systems (WDS).In order to accelerate sy stem transformation f rom unsteady state to steady state , tw o mathematical models w ere proposed based on the theory of transient flow , which are the method of characteristics (MOC)with re-divide damping coefficient and the lattice Boltzmann equation (LBE).These methods were used to calculate an ex ample in WDS .It is show n that LBE is more accurate when Rey nolds number is less than 7 000 , o therwise MOC is mo re precise. KEY WORDS w ater distributio n sy stem ;unsteady flow ;method o f characteristics;lattice Boltzmann equation 收稿日期:2008-05-22 基金项目:国家“十一五”科技支撑计划重大资助项目(No .2006BAJ03A06) 作者简介:马 飞(1968—), 男, 教授, 博士, E-mail:ustbyeke@sina.com 由于给水管网是动态的非线性系统 .管网的运 行状态往往是从一个稳定的运行状态过渡到另一个 稳定的运行状态[ 1] .建立在恒定流动理论基础之上 的管网水力模型是非常成熟的 , 在管网的平稳运行 工况建模中 ,这种方法能够精确计算管网各控制点 的压力 、流量等运行参数 .然而, 在过渡阶段这种方 法很难准确地确定管网的运行工况.因此 , 本文从 非恒定流动的角度出发, 分别利用重分阻尼系数法 划分时间步长的特征线法(method of characteris￾tics , MOC)和引入摩阻附加项的玻尔兹曼网格法 (lattice Boltzmann equation , LBE)求解城市给水管 网系统非恒定流动数学模型, 进一步完善给水管网 动态分析方法. 1 给水管网非恒定流动物理模型 非恒定流动中液体质点的运动要素不仅随空间 位置变化 ,而且随时间过程变化[ 2] .因此分析这类 运动 ,要考虑由于运动要素随时间变化所引起的惯 性力的作用 , 即其运动要素为时间和空间的函数. 由质量和能量守恒原理得管道非恒定流动控制方 程 [ 3-5] : H t +a 2 g V x =0 (1) V t +g H x + f 2D V V =0 (2) 式中 , V 为管内某断面的平均流速 , H 为某点测压 管水头, D 为管径, x 为沿流线的位置坐标 , t 为时 第 31 卷 第 4 期 2009 年 4 月 北 京 科 技 大 学 学 报 Journal of University of Science and Technology Beijing Vol .31 No.4 Apr.2009 DOI :10.13374/j .issn1001 -053x.2009.04.019

。424 北京科技大学学报 第31卷 间坐标，g为重力加速度，∫为管内水流的沿程阻力 重分阻尼损失系数为： 系数，α为水锤波在管段内的行波速度. fLy f= (7) 2采用重分阻尼系数法划分时间步长的特 aAtintar△+0.5 征线法 对于复杂的给水管网，此公式通过对管道进行 特征线法是将方程(1)和(2)在计算空间和时间 分段误差控制，在满足时步一致性和不调整管道水 上进行离散，其特征线方程为： 击波速的前提下，满足伯努利水头损失方程在重分 c:(Qs)+()+ 前后的一致性. 此法在时步处理上有如下改进： f0DlOP△x=0 (1)保证调整后的管道和原管道理论水击波速 (3) 2gDA2 相同，即没有改变管道固有水击特性； Cg号(0-Q）-(Hp-He+ (2)重分阻尼系数法使调整后的管道和原来管 f0mlQm△x=0 道总阻尼损失相同，严格满足伯努利方程； 2gDA (4) (3)此方法在计算机编程中很容易实现稳定 性很容易得到满足， 式中，C为t一△t到t方向特征线C为t到t十 △t方向特征线，Qp、Hp为计算点在t时刻的流量 3引入摩阻附加项的玻尔滋曼网格法 和压力，QA、HA为计算点在前一网格1一△1时刻 (LBE) 的流量与压力，QB、HB为计算点后一网格t十△1 玻尔兹曼网格模型(LBE)基于流体宏观运动速 时刻的流量与压力，A为计算管过流断面面积，△x 度为粒子运动的平均速度的原理，将微观与宏观模 为距离步长，QP为A点与P点间的平均流速，Q即 型相结合，提供了间接求解宏观模型方程的方 为B点与P点间的平均流速. 由于给水管网是由多根管道组成的复杂串联管 法.不同于传统的直接求解宏观偏微分方程数值 解的计算方法，LBE是先获得玻尔兹曼变量的一般 道系统，只有当时间步长取很小值时，才能满足步长 解，宏观的状态变量是通过对玻尔兹曼变量积分获 的一致性.这样不仅延长了计算时间，而且给计算 的精度和稳定性带来了很大影响.因此，本文采用 得.其计算过程如图1所示. 一致性 重分阻尼系数法选取步长.即根据计算精度的需 波尔兹曼微分方程 宏观守恒率 要，不断调整时间步长，根据时间步长将管道分为整 波尔兹曼 差分方程 数段，同时将原管道的阻尼再重新分配到调整后加 波尔兹曼方程的解 反算 一宏观量的解 长或缩短的计算管道上.这样使整个管道的阻尼损 图1LBE法求解过程 失和原管道的等价，调整后管道依旧严格满足伯努 Fig.I Comput at ional process of LBE 利水头损失方程，而波速可用求得的真实水击波速. 初选△1，将管网中计算管段分别进行分段，如 粒子分布方程满足玻尔兹曼网格方程(LBE): 多余长度△L≤05a△t,则忽略不计.否则N= fa(x+ea△i,t十△t)-fa(x,t)=n+n(8) N+1,N为： 式中，f。是粒子分布函数：e。是粒子的运动速度向 Nj=int +0.5 量，即相邻网格点连线形成的向量，一维运动，粒子 (5) 运动速度为e={0,1,一1：a表示粒子运动方向， 式中，N为每根管道的分段数，L,为重分阻尼之前 粒子共有b十1个允许的运动方向：2为碰撞函数， 每根管道的长度，a,为重分阻尼前每根管道的波 0=-【fa(x)-f公(x,:8为粒子在(x 速. t)处的平衡分布函数：2。为与摩阻相关的附加项： 当波速不变时，等价变换后管道的计算长度为： τ为松弛时间系数. LJ=ay△int ar,+0 (6) LBE方程9： 式中，L]为调整后管道的计算长度，称之为虚拟长 度 +, a=0 亚(u) (9) 三0

间坐标, g 为重力加速度 , f 为管内水流的沿程阻力 系数 , a 为水锤波在管段内的行波速度 . 2 采用重分阻尼系数法划分时间步长的特 征线法 特征线法是将方程(1)和(2)在计算空间和时间 上进行离散,其特征线方程为[ 6-7] : C +: a gA(QP -QA)+(HP -HA)+ fQAP QAP Δx 2gDA 2 =0 (3) C - : g gA (QP -QB)-(HP -HB)+ fQBP QBP Δx 2gDA 2 =0 (4) 式中 ,C +为 t -Δt 到 t 方向特征线, C -为 t 到 t + Δt 方向特征线 , QP 、HP 为计算点在 t 时刻的流量 和压力 , QA 、HA 为计算点在前一网格 t -Δt 时刻 的流量与压力, QB 、HB 为计算点后一网格 t +Δt 时刻的流量与压力 , A 为计算管过流断面面积, Δx 为距离步长 , QAP为 A 点与 P 点间的平均流速, QBP 为 B 点与 P 点间的平均流速. 由于给水管网是由多根管道组成的复杂串联管 道系统,只有当时间步长取很小值时,才能满足步长 的一致性.这样不仅延长了计算时间, 而且给计算 的精度和稳定性带来了很大影响.因此 , 本文采用 重分阻尼系数法选取步长 .即根据计算精度的需 要,不断调整时间步长,根据时间步长将管道分为整 数段, 同时将原管道的阻尼再重新分配到调整后加 长或缩短的计算管道上.这样使整个管道的阻尼损 失和原管道的等价 ,调整后管道依旧严格满足伯努 利水头损失方程 ,而波速可用求得的真实水击波速 . 初选 Δt ,将管网中计算管段分别进行分段 , 如 多余长度 ΔL ≤0.5a JΔt , 则忽略不计 .否则 N = NJ +1 , NJ 为: NJ =int LJ aJΔt +0.5 (5) 式中, NJ 为每根管道的分段数, LJ 为重分阻尼之前 每根管道的长度, a J 为重分阻尼前每根管道的波 速. 当波速不变时, 等价变换后管道的计算长度为 : L′J =aJΔtint LJ aJΔt +0.5 (6) 式中, L′J 为调整后管道的计算长度, 称之为虚拟长 度. 重分阻尼损失系数为: f′= fLJ aJΔtint LJ aJΔt +0.5 (7) 对于复杂的给水管网 ,此公式通过对管道进行 分段误差控制 ,在满足时步一致性和不调整管道水 击波速的前提下, 满足伯努利水头损失方程在重分 前后的一致性. 此法在时步处理上有如下改进: (1)保证调整后的管道和原管道理论水击波速 相同,即没有改变管道固有水击特性 ; (2)重分阻尼系数法使调整后的管道和原来管 道总阻尼损失相同 ,严格满足伯努利方程 ; (3)此方法在计算机编程中很容易实现, 稳定 性很容易得到满足 . 3 引 入 摩 阻 附 加 项 的 玻 尔 兹 曼 网 格 法 (LBE) 玻尔兹曼网格模型(LBE)基于流体宏观运动速 度为粒子运动的平均速度的原理, 将微观与宏观模 型相结合 , 提供了间接求解宏观模型方程的方 法 [ 8] .不同于传统的直接求解宏观偏微分方程数值 解的计算方法 , LBE 是先获得玻尔兹曼变量的一般 解 ,宏观的状态变量是通过对玻尔兹曼变量积分获 得 .其计算过程如图 1 所示. 图 1 LBE 法求解过程 Fig.1 Comput ational process of LBE 粒子分布方程满足玻尔兹曼网格方程(LBE): fα(x +eαΔt , t +Δt)-fα(x , t)=Ψα+Ψ′α (8) 式中, fα是粒子分布函数;eα是粒子的运动速度向 量 ,即相邻网格点连线形成的向量 , 一维运动 ,粒子 运动速度为 eα={0 , 1 , -1};α表示粒子运动方向, 粒子共有 b +1 个允许的运动方向 ;Ψα为碰撞函数, Ψα=- 1 τ [ fα(x , t)-f 0 α(x , t)] ;f 0 α为粒子在(x , t)处的平衡分布函数;Ψ′a 为与摩阻相关的附加项; τ为松弛时间系数. LBE 方程[ 9] : ∑ b α=0 f 0 α t + ∑ b α=0 eαf 0 α xj = ∑ b α=0 Χ(u) (9) · 424 · 北 京 科 技 大 学 学 报 第 31 卷

第4期 马飞等：给水管网非恒定流动数值计算方法 425。 旁+-0 (12) =0 一=0 (10) 如果假设h与v分别为： 当设约束条件为人=会=会= h=gH (13) v=a2V 则式(11)和(12)即可与式(1)和(2)达成统一. 之=v时LBE方程的宏观模型形式为： 4实例分析 图2为某小区管网简图.系统总流量为 ++办川=0 (11) 0.29m3s1. 0.021 ⑥ 0.031 5 ③ 91 50-6 ② 2501160 200-760 0.034 0.021 ① ① 300-610 1 ⑦ 200-975 ③ 250-1067 6/0.028 ④ o 0.028 200-915 0.023 ① ⑧ 200 ⑧ 250-975 70.024 ⑤ 1四/ 6.021 ⑤ 200-1067 300-610 0.018 250-1220 200-1160 ② ⑨ 1 ⑨ 200-975 0.024 0.021 © ☑ 图2某小区管网（各管段的管径和管长显示单位分别为mm和m) Fig 2 Exampe of pipe network (the units of pipe diameter and length are mm and m) 在初始情况下，水库HP(0,t)=150m时，系统 库水位的高度 平均流速为0.85ms1,最大流速为2.06ms1,最 小流速为0.21ms1.平均压头为48.52m,最大水 头为56.44m,最小水头为44.27m.节点9的流量 70 60 随时间变化如表1所示. 50 表1节点9流量变化 40 Table I Demand of Node 9 with time 时间/s 0 36912 20 10 流量/(m3·s-1)0034004200500570071 2 10 6 节点9的流量在12s内增量为0.037m3·s1, 时间 由于管网较小，节点9的瞬时变化对与其相连的节 点5、6和10影响相对较大.分别利用M0C法与 图3H=150m时M0C法计算各点压力值 LBE法计算节点9、10、5、6、7和1压力随时间变化 Fig 3 Node pressure when H=150m by MOC 如图3与图4所示.由图3和图4分析可知，当节 点9流量连续变化时将引起其他各节点压力和管段 在本管网的节点1、5、6、7、9和10的消火栓处 流量不同程度的变化两种算法得到管段上各计算 设置压力传感器，进行压力时时监测记录.将节点9 点压力随时间变化趋势相同.节点9的动作对与其 与节点10实测压力值与MOC法和LBE法计算值 相邻节点6和节点10压力影响最大，而对节点1的 比较如图5和图6所示.由图5可知，两种算法与 影响最小.此时，节点1的压力变化主要决定于水 实测值在节点9处极值发生的时刻相同，最大发生

∑ b α=0 eαf 0 α t + ∑ b α=0 eαeαf 0 α x j =0 (10) 当设约束条件为 ∑ b α=0 fα= ∑ b α=0 f 0 α=v , ∑ b α=0 eαfα= ∑ b α=0 eαf 0 α=h , ∑ b α=0 Χ(u)=- f 2D v v , ∑ b α=0 eαeαf 0 α= ∑ b α=0 f 0 α=v 时 , LBE 方程的宏观模型形式为: v t + h x + f 2D v v =0 (11) h t + v x =0 (12) 如果假设 h 与 v 分别为 : h =gH v =a 2 V (13) 则式(11)和(12)即可与式(1)和(2)达成统一. 4 实例分析 图 2 为某 小 区管 网 简 图.系统 总 流量 为 0.29 m 3·s -1 . 图 2 某小区管网(各管段的管径和管长显示单位分别为 mm 和 m) Fig.2 Example of pipe network (the units of pipe diameter and length are mm and m) 在初始情况下, 水库 HP(0 , t)=150 m 时 ,系统 平均流速为0.85m·s -1 , 最大流速为2.06m·s -1 , 最 小流速为 0.21m·s -1.平均压头为 48.52m , 最大水 头为 56.44 m ,最小水头为 44.27 m .节点 9 的流量 随时间变化如表 1 所示. 表 1 节点 9 流量变化 Table 1 Demand of Node 9 with time 时间/ s 0 3 6 9 12 流量/(m 3·s -1) 0.034 0.042 0.05 0.057 0.071 节点 9 的流量在 12 s 内增量为 0.037 m 3·s -1 , 由于管网较小, 节点 9 的瞬时变化对与其相连的节 点5 、6 和 10 影响相对较大 .分别利用 MOC 法与 LBE 法计算节点 9 、10 、5 、6 、7 和 1 压力随时间变化 如图 3 与图 4 所示.由图 3 和图 4 分析可知 , 当节 点 9 流量连续变化时将引起其他各节点压力和管段 流量不同程度的变化, 两种算法得到管段上各计算 点压力随时间变化趋势相同.节点 9 的动作对与其 相邻节点 6 和节点 10 压力影响最大,而对节点 1 的 影响最小.此时, 节点 1 的压力变化主要决定于水 库水位的高度. 图 3 H =150 m 时 MOC 法计算各点压力值 Fig.3 Node pressure w hen H =150 m by MOC 在本管网的节点 1 、5 、6 、7 、9 和 10 的消火栓处 设置压力传感器 ,进行压力时时监测记录 .将节点9 与节点 10 实测压力值与 MOC 法和 LBE 法计算值 比较如图 5 和图 6 所示 .由图 5 可知 ,两种算法与 实测值在节点 9 处极值发生的时刻相同 , 最大发生 第 4 期 马 飞等:给水管网非恒定流动数值计算方法 · 425 ·

。426 北京科技大学学报 第31卷 如图7和图8所示.由图7和图8可知，管网中各 点利用LBE法计算结果与实测结果相对值HBE/ 70 H。在L.07~0.91范围内，绕HE/H。=1的曲线 上下脉动，其中H。为实测值.利用MOC法计算结 50 果与实测结果相对值HMoc/H。在1.28~084范 40 围内，存在着某些节点实测值与计算值相对误差较 30 20 大情况，最大误差为28%.由此分析，此时LBE法 要优于MOC法. 12 13 时间S 图4H=150m时LBE法计算各点压力值 Fig 4 Node pressure when H=150m by LBE 09 在4s时刻，最小值在2s时刻.LBE法计算的相对 误差为1.08%，M0C法计算的相对误差为41%， 6 1012 时间s LBE法优于MOC法.由图6可知，在节点10处 LBE法计算相对误差为2.8%，MOC法计算的相对 图7MOC法计算各点压力值与实测值比较 Fig.7 误差为3.6%两种算法基本相同.为了比较两种算 Comparison of measured and simulated pressure by MOC 法的准确性，将各点两种算法压力值与实测值比较， 3 。96 ◆-10 4-7 70 -TEST 00 1.0 0.9 0.8 30 10 12 时间s 10 时间/s 图8BE法计算各点压力值与实测值比较 Fig.8 Comparison of measured and simulated pressure by LBE 图5节点9计算值与实测值比较 Fig.5 Comparison of calculation and experimental pressure valwes at 当水库Hp(0,t)=270m时，节点5、6、7、10和 Node 9 9的流量压力变化情况，分别如图9和图10所示. 180 50 160 40 140 MOC 120 LBE 100 20 TEST 0 1012 时间s 时间 图6节点10计算值与实测值比较 Fig.6 Comparison of calculation and experimental presaure values at 图9H=270m时M0C法计算各点压力值 Node 10 Fig 9 Node pressure when H=270m by MOC

图 4 H =150 m 时 LBE 法计算各点压力值 Fig.4 Node pressure w hen H =150 m by LBE 在4 s 时刻 ,最小值在 2 s 时刻.LBE 法计算的相对 误差为 1.08 %, MOC 法计算的相对误差为 4.1 %, LBE法优于 MOC 法 .由图 6 可知 , 在节点 10 处 LBE 法计算相对误差为 2.8 %,M OC 法计算的相对 误差为3.6 %, 两种算法基本相同 .为了比较两种算 法的准确性,将各点两种算法压力值与实测值比较 , 图 5 节点 9 计算值与实测值比较 Fig.5 Comparison of calculation and experimental pressu re values at Node 9 图 6 节点 10 计算值与实测值比较 Fig.6 Comparison of calculation and experimental pressu re values at Node 10 如图 7 和图 8 所示 .由图 7 和图 8 可知, 管网中各 点利用 LBE 法计算结果与实测结果相对值 HLBE/ Ho 在 1.07 ～ 0.91 范围内 ,绕 HLBE/ Ho =1 的曲线 上下脉动 ,其中 Ho 为实测值 .利用 MOC 法计算结 果与实测结果相对值 H MOC/ Ho 在 1.28 ～ 0.84 范 围内 ,存在着某些节点实测值与计算值相对误差较 大情况, 最大误差为 28 %.由此分析 , 此时 LBE 法 要优于 MOC 法 . 图 7 MOC 法计算各点压力值与实测值比较 Fig.7 Comparison of measured and simulated pressure by MOC 图 8 LBE 法计算各点压力值与实测值比较 Fig.8 Comparison of measured and simulated pressure by LBE 当水库 HP(0 , t)=270 m 时, 节点 5 、6 、7 、10 和 9 的流量压力变化情况, 分别如图 9 和图 10 所示 . 图 9 H =270 m 时 MOC 法计算各点压力值 Fig.9 Node pressure w hen H =270 m by MOC · 426 · 北 京 科 技 大 学 学 报 第 31 卷

第4期 马飞等：给水管网非恒定流动数值计算方法 427。 计算给水管网非恒定流动，将提高计算的准确性，同 时提高计算效率. 200 150 100 参考文献 50 0 [I]Tao J K.Studies on dynamic hydraulic model of shanghai munici- -50 雷 -100 pal water distribu tion netw ork.China Water Wastewater,1999. -150 15(4):11 200 -250 (陶建科.建立上海市计算机给水管网动态水力模型研究.中国 给水排水，1999,15(4)：11) D A202 时间 Wylie E B.Fluid Transients in Systems.Engewood.Prentice B Hal.1993 3 Mohamed S G.Bryan W K.Duncan A M.Energy estim ates for 图10H=270m时LBE法计算各点压力值 discretization errors in w ater hammer problems.J Hydraul Eng, Fig 10 Node pressure when H=270m by LBE 1998124(4):384 [4 Sericola B.Transient analysis of stochastic fluid models.Perform 比较图9与图10，发现M0C法计算结果与 Eval.199832(4):245 LBE法计算结果相差较大.此时，MOC法计算结果 [5 Uspuras E.Kaliatka A.Dunduis G.Analysis of potential water hammer at Ignalina NPP using thermal-hydraulic and structumal 与水库水位为150m时压头变化趋势相同，并符合 analysis codes Nucl Eng Des.2001,203(1):1 流动变化趋势：而LBE法计算结果不稳定.当水库 [6 Cheng X Y.Computational F luid Dynamics:Numerical Solu- 水位为270m时，此时初始计算时Re=15654.这 tion of Partial Differential Eqation.Beijing:Science Press 与Hou等10利用256×256网格计算流体流动得出 2000 (程心一.计算流体动力学：偏微分方程的数值解法.北京：科 的Re从1.0到10000，压力与其他方法结果误差小 学出版社，2000) 于5%，当Re7500时LBE容易失稳的结论相符 I7 Fox JA.Hydrau lic Ana lysis for Unsteady Flow of Water Sup- 合.因此，在利用LBE计算管网非恒定流动过程 ply Pipeline.Beijing:Petrokum Industry Press 1996 中，应时刻根据计算结果计算不同时刻的Re. (福克斯」A.管网中不稳定流动的水力分析.陈祖泽译。北 京：石油工业出版社.1996) 5结论 8 QuS L Yuan Y X.Wu Y B.et al.Analysis on transient flow of long-distance water transmission pipelne.China Water Wasew- (1)提出了利用重分阻尼系数法进行时间步长 aler,2005.21(12):59 的特征线法(M0C),此方法在不改变管道固有水击 (曲世琳袁一星，伍悦滨，等.长距离输水管线非恒定流动分 特性的前提下，保证调整后管道与原管道水击波速 析.中国给水排水，2005,21(12)：59) 相同，并且严格满足伯努利水头损失方程. [9 QuS L.Zhao H B Wu Y B.Lattice Boltzmann equation for (2)引入摩阻附加项的玻尔兹曼网格法(LBE) transient flow of long-distance water supply pipeline.J Nanjing Univ Sci Technol,2004.28(6):663 求解简单，几何边界易处理，流场计算时不须网格转 (曲世琳赵洪滨，伍悦滨.玻尔兹曼网格方程在长输管道瞬 换因此很方便进行编程. 变流动中的应用.南京理工大学学报，2004,28(6)：663) (3)实例分析表明，Re≤7000时采用玻尔兹曼 [10 Hou S Zou Q.ChenS et al.Simulation of cavity flow by the 网格法(LBE),Re>7O00时采用特征线法(MOC) lattice Boltzmann method.J Comput Phys.1995.118(2):329

图10 H =270 m 时 LBE 法计算各点压力值 Fig.10 Node pressure w hen H =270 m by LBE 比较图 9 与图 10 , 发现 M OC 法计算结果与 LBE 法计算结果相差较大.此时 ,MOC 法计算结果 与水库水位为 150 m 时压头变化趋势相同, 并符合 流动变化趋势;而 LBE 法计算结果不稳定.当水库 水位为 270 m 时, 此时初始计算时 Re =15 654 .这 与Hou 等 [ 10] 利用 256 ×256 网格计算流体流动得出 的 Re 从 1.0 到 10000 ,压力与其他方法结果误差小 于 5 %, 当 Re >7 500 时 LBE 容易失稳的结论相符 合.因此 , 在利用 LBE 计算管网非恒定流动过程 中,应时刻根据计算结果计算不同时刻的 Re . 5 结论 (1)提出了利用重分阻尼系数法进行时间步长 的特征线法(MOC),此方法在不改变管道固有水击 特性的前提下, 保证调整后管道与原管道水击波速 相同 ,并且严格满足伯努利水头损失方程. (2)引入摩阻附加项的玻尔兹曼网格法(LBE) 求解简单 ,几何边界易处理,流场计算时不须网格转 换,因此很方便进行编程 . (3)实例分析表明 , Re ≤7 000 时采用玻尔兹曼 网格法(LBE), Re >7 000 时采用特征线法(MOC) 计算给水管网非恒定流动, 将提高计算的准确性 ,同 时提高计算效率. 参 考 文 献 [ 1] Tao J K.Studies on dynami c hydraulic model of shanghai munici￾pal w ater distribu tion netw ork.China Wa ter Wastewater , 1999 , 15(4):11 (陶建科.建立上海市计算机给水管网动态水力模型研究.中国 给水排水, 1999 , 15(4):11) [ 2] Wylie E B .Flu id Transients in Systems .Englewood:Prentice Hall, 1993 [ 3] Mohamed S G , Bryan W K , Duncan A M .Energy estim ates for discretization errors in w at er hammer problems.J Hyd raul Eng , 1998 , 124(4):384 [ 4] Sericola B.Transient analysis of stochastic fluid models.Perform Eval , 1998 , 32(4):245 [ 5] Uspuras E , Kaliatka A , Dunduis G.Analysis of potential water hammer at Ignalina NPP using thermal-hydrauli c and structu ral analysis codes.Nucl Eng Des, 2001 , 203(1):1 [ 6] Cheng X Y.Compu tationa l F luid Dynamics :Nu merical Solu￾tion of Partial Diff erential E quation.Beijing :Science Press, 2000 (程心一.计算流体动力学:偏微分方程的数值解法.北京:科 学出版社, 2000) [ 7] Fox J A .Hydrau lic Ana lysis for Unsteady Flow of Water S up￾p ly Pipeline .Beijing :Petroleum Industry Press, 1996 (福克斯 J A .管网中不稳定流动的水力分析.陈祖泽, 译.北 京:石油工业出版社, 1996) [ 8] Qu S L , Yuan Y X, Wu Y B , et al.Analysis on transient flow of long-distance w ater transmission pipeline .China Wat er Wastew￾ater , 2005 , 21(12):59 (曲世琳, 袁一星, 伍悦滨, 等.长距离输水管线非恒定流动分 析.中国给水排水, 2005 , 21(12):59) [ 9] Qu S L , Zhao H B, Wu Y B .Lattice Boltzmann equation for transient flow of long-distance w at er supply pipeline .J Nanjing U niv Sci Technol , 2004 , 28(6):663 (曲世琳, 赵洪滨, 伍悦滨.玻尔兹曼网格方程在长输管道瞬 变流动中的应用.南京理工大学学报, 2004 , 28(6):663) [ 10] Hou S , Zou Q , Chen S , et al.Simulation of cavit y flow by the lattice Boltzmann method.J Comp ut Phys, 1995 , 118(2):329 第 4 期 马 飞等:给水管网非恒定流动数值计算方法 · 427 ·