基于涡黏度模型的冰封河道纵向时均流速垂向分布解析解

２０２３年１２月
水 利 学 报
ＳＨＵＩＬＩ ＸＵＥＢＡＯ第５４卷　第１２期
文章编号：０５５９－９３５０（２０２３）１２－１４７３－１２
收稿日期：２０２３－０４－１０；网络首发日期：２０２３－１１－０７
网络首发地址：ｈｔｔｐｓ：??ｋｎｓ．ｃｎｋｉ．ｎｅｔ?ｋｃｍｓ?ｄｅｔａｉｌ?１１．１８８２．ＴＶ．２０２３１１０６．１０３１．００３．ｈｔｍｌ
基金项目：中国博士后科学基金项目（２０２２Ｍ７１２９０４，２０２２Ｍ７２２８８１）；河南省重点研发与推广专项（科技攻关）项目（２３２１０２３２１０９９）；黄河实
验室（郑州大学）一流课题项目（
ＹＲＬ２２ＹＬ０３，ＹＲＬ２２ＩＲ１１）作者简介：王菲菲（１９９１－），博士，讲师，主要从事河流水冰沙动力学研究。

Ｅ－ｍａｉｌ：ｆｆｗａｎｇ＠ｚｚｕ．ｅｄｕ．ｃｎ
通信作者：李志伟（１９８５－），博士，副教授，主要从事河流动力过程及生态环境响应研究。

Ｅ－ｍａｉｌ：ｚｗｌｉ＠ｚｚｕ．ｅｄｕ．ｃｎ
基于涡黏度模型的冰封河道纵向时均流速垂向分布解析解
王菲菲１，２，槐文信３，陈海亮２，孙　斌１，２，李志伟１
，２
（１．郑州大学黄河实验室，河南郑州　４５０００１；２．郑州大学水利与交通学院，河南郑州　４５０００１；
３．武汉大学水资源与水电工程科学国家重点实验室，湖北武汉　４３００７２）
摘要：冰盖的存在显著改变了河道水流结构，使流速垂向分布关于零切应力平抛物线面呈型。

考虑到准确预测冰盖下水流流速分布是计算冰封河道流量、泥沙输移率和河床演变预估的基础，本文应用双层假定和冰盖流涡黏度模型，推导出一个物理意义明确且速度梯度连续的冰封河道纵向时均流速垂向分布解析解，并明确常见的冰盖流对数流速公式和双幂律流速公式的应用条件。

采用水槽实测数据和原型实验数据，对比分析所提出的解析解与现有对数公式和双幂律公式的有效性和精确性，并基于本文解析解探究了纵向时均流速垂向分布对各物理特征参数的敏感性。

结果表明，该解析解的流速计算值与实测值吻合较好且较现有对数公式和双幂律公式具有较高精度；雷诺数的较小变化并不会引起纵向时均流速垂向分布的较大改变，但流速分布对冰盖糙率与河床糙率的相对值比较敏感，且纵向时均流速的最大值偏向于较光滑边界一侧，该结论与以往基于ｋ－ε湍流模型所得到的数值结果是一致的，再次证实了本文解析解的有效性。

关键词：冰封河道；纵向时均流速；垂向分布；涡黏度
中图分类号：ＴＶ１３３文献标识码：Ａｄｏｉ：１０．１３２４３?ｊ．ｃｎｋｉ．ｓｌｘｂ．２０２３０２００
１　研究背景
高纬度地区河渠在冬季受低温影响容易结冰，形成冰盖［１－
２］。

冰盖的出现使河渠畅流期的自由水面被附加固壁所代替，改变了河渠的几何边界条件、热力及水力条件，形成冬季河渠特有的冰情现象［３］。

冰情问题是冬季高纬度地区河渠所面临的一个共性问题，例如，美国的密西西比河、加拿大的圣劳伦
斯河、俄罗斯的勒拿河和叶尼塞河等都是冰情相对严重的河流［４
］。

在我国东北、新疆和内蒙等高纬度地区，部分河渠在每年１２月至次年３月均面临结冰问题，在封河初期及开河期易发生冰凌灾害［５－
７］。

据相关资料统计，在１
９５１年—２０１０年间黄河宁蒙河段发生了将近３０年的凌汛灾害，沿岸居民不仅生活条件深受影响，而且遭受了巨大的经济损失；而黑龙江流域局部河段大约每３年形成一次较大规模
冰坝［８］。

实际工程中，冰情所引起的一系列水力条件变化还可能会诱发自然灾害的发生或带来比较严
重的威胁。

鉴于此，冰盖是冬季河渠水流特性研究中不容忽视的因素，其不仅影响河流水体温度与悬
浮物质的分布，而且影响水库坝体结构稳定和长距离输水工程的正常运行［９－
１
３］
，迫切需要全面了解冰封河渠的水流特性并合理预测以减少经济损失。

从水力学角度来看，天然河渠在每年冰封期常常伴随着水径流和河道流量的减少，从而引起河道
输沙量和物质输运率的降低［１
４］。

与明渠流相比，冰封河渠具有不同的水力特征，冰盖的存在施加了一个额外的阻力边界，该边界使断面湿周成倍增加、水深和总边界切应力增大，导致断面平均流速降
低、河道流量减少［１５］。

受冰盖底面和河床底部粗糙度的共同影响，纵向时均流速的垂向分布呈抛物线
型，并且关于零切应力所在平面通常是非对称的，最大纵向时均流速朝向较光滑的边界一侧，最大纵向时均流速位置处的剪切应力最小，并且剪切应力随着远离该位置而逐渐增大。

考虑到流速分布预测
是计算河道综合糙率、预测泥沙输移及河床形态演变的基础［１６－
２０］
，故探究河流冰封期纵向时均流速垂
向分布规律及其精确可靠的计算公式具有十分重要的理论意义和应用价值。

国内外学者针对冰盖流的纵向时均流速垂向分布规律采用不同的方法开展了较为丰富的研究。

一
方面，通过开展系列水槽实验，陈建国等［２１］探究了连续和不连续冰盖流条件下纵向时均流速的垂向分布；陈刚等［２２］对比研究了明流、完全冰封和部分冰封条件下矩形河道的纵向时均流速和紊动特性的垂向分布规律；Ｋ
ｉｍｉａｇｈａｌａｍ等［２３］
研究了不同冰盖覆盖度和边界粗糙度对梯形河道流速分布的影响。

另一方面，基于ｋ－ε湍流模型，茅泽育等［２４］和王志兴等［２５］分别建立了冰盖流的立面二维和三维数值模型以探究影响流速垂向分布规律的主导因素；Ｅ
ｓｓｅｌ等［２６］
将三维数值模拟结果与实测结果相结合，阐明了不同冰盖覆盖度下矩形河道的流速分布规律。

虽然ｋ－ε湍流模型能较准确计算冰封河渠纵向时均
流速的垂向分布［２
７］
，但其计算程序复杂且耗时长，故快速准确的解析模型仍是学者研究的重要课题。

鉴于此，国内外学者基于Ｅｉｎｓｔｅｉｎ阻力划分理论［２８］
提出了冰封河道的“双层假定”，进而将明渠流对
数流速公式和指数流速公式分别应用于冰盖流，建立冰封河道流速垂向分布的解析模型，即对数公式
和双幂律公式［２９－
３１］。

对于现有的冰封河道流速解析模型，对数流速公式在最大流速位置处具有速度梯
度不连续的缺点，而双幂律流速公式未明确涉及冰盖流的两组速度和长度尺度，其物理意义不直观。

本文旨在应用“双层假定”与冰盖流涡黏度模型，提出一个物理意义明确且处处满足连续光滑条件的冰封河道纵向时均流速垂向分布的解析解；其次，应用冰封河道水槽实测流速数据，对比分析该解析解与现有的对数公式和双幂律公式预测冰盖下纵向时均流速垂向分布的有效性和准确性；最后，基于本文所提出的解析解，探究各物理特征参数对冰盖流纵向时均流速垂向分布的影响规律。

２　纵向时均流速垂向分布公式
常用的冰封河道纵向时均流速垂向分布公式主要包括对数公式和双幂律公式。

上述公式均是基于
“双层假定”给出的［２
９］
，该假说已被证明可用来描述冰盖流的纵向时均流速垂向分布，它基于这样的假定：冰盖下水流可被划分为具有伪自由水面的两个独立流层（即河床区和冰盖区，见图１），划分平面通常被称为“零切应力平面”，各流层内的水流要素分别只受冰盖底面和河床底部粗糙度的影响，即各流层内的纵向时均流速Ｕ和断面上切应力τｚｘ的垂向分布与明渠流类似，两层交界面处具有最大纵向时均流速Ｕｍａｘ和零切应力，冰盖底面和和河床底部的流速为零，而切应力τｚｘ在冰盖底面和河床底部分别具有最大值τｉ和τｂ。

本节首先介绍冰盖流纵向时均流速垂向分布的对数公式和双幂律公式及其应用条件，进而基于“双层假定”和冰盖流涡黏度模型，推导出一个新的纵向时均流速垂向分布解析公式。

本文若无明确说明，所指流速均为纵向时均流速。

图１　冰封河道纵向时均流速Ｕ和切应力τｚｘ
的垂向分布示意图
２．１　对数公式　冰盖流纵向时均流速对数公式是基于“双层假定”将Ｐｒａｎｄｔｌ紊流对数流速公式直接
应用到河床区（即０＜ｚ≤ｈｂ）和冰盖区（即ｈｂ
≤ｚ＜Ｈ）中得到的，具体为［２９］：
＝１κｌｎｚｈｂ，０＜ｚ≤ｈｂ
＝１κｌｎＨ－ｚｈｉ，ｈｂ
≤ｚ＜Ｈ（１）
式中：Ｕｂ和Ｕｉ分别为河床区和冰盖区的纵向时均流速，ｍ?ｓ；ｕ ｂ和ｕ ｉ分别为河床区和冰盖区的摩阻流速，ｍ?ｓ；Ｕｍａｘ为纵向时均流速沿垂向分布的最大值，ｍ?ｓ；ｈｂ和ｈｉ分别为河床区和冰盖区的水深，ｍ；Ｈ为总水深且Ｈ＝ｈｂ＋ｈ
ｉ，ｍ；ｚ为垂向坐标，ｍ；κ为卡门常数且κ＝０．４。

于是，根据连续介质假定，对数流速公式（１）需满足以下条件：
（１）河床区和冰盖区交界面处的速度连续，即当ｚ＝ｈｂ时，Ｕｂ＝Ｕｉ＝Ｕｍａｘ。

（
２）河床区和冰盖区交界面处的速度梯度连续，即ｄＵｂ
ｄｚｚ＝ｈｂ＝ｄＵｉ
ｄｚｚ＝ｈｂ。

２．２　双幂律公式　虽然对数流速公式在河床区和冰盖区交界面处满足速度连续条件，但难以满足速度梯度连续条件（后续进行详细讨论）。

也就是说，基于对数公式，预测的纵向时均流速垂向分布在交界
面处是连续不光滑的。

鉴于此，Ｔｓａｉ［３０］
基于明渠流的指数流速公式，提出了冰盖流纵向时均流速垂向
分布的双幂律公式：
Ｕ＝Ｋ０
ｚ
Ｈ
()１ｍｂ１－ｚＨ
()
１ｍｉ
（２）
式中：Ｋ０为与单宽流量ｑ相关的待定参数，ｍ?ｓ；ｍｂ和ｍｉ分别为与河床粗糙度和冰盖粗糙度相关的参数；其余物理量含义与对数公式（１）中的相同，这里不再赘述。

当ｍｉ趋于无穷大时，式（２）即化为明渠流纵向时均流速的指数公式，此时Ｋ０与Ｕｍａｘ
相对应。

对式（２）进行求导，可得速度梯度为：
ｄＵｄｚ＝Ｕ
１ｍｂｚ－１
ｍｉ（Ｈ－ｚ）
[]
（３）
显然，速度梯度关于垂向坐标ｚ在区间（０，Ｈ）上是连续的，意味着双幂律公式（２）表征的纵向时均流
速Ｕ沿垂向是连续光滑的。

为了应用双幂律公式预测冰封河道纵向时均流速的垂向分布，首先需要获知粗糙度参数ｍｂ和ｍｉ，以及流量参数Ｋ０
，下面将详细介绍这三个物理量的确定过程。

令ｄＵ?ｄｚ＝０，可得最大纵向时均流速Ｕｍａｘ
的垂向位置，即河床区的水深：ｈｂ＝ｚｍａｘ
＝Ｈ
１＋ｍｂ?ｍｉ
（４）
根据Ｔ
ｓａｉ［３０］
的研究，参数ｍｂ和ｍｉ可分别表示为：ｍｂ＝κ
８
ｆｂ
()
０．５
，ｍｉ＝κ
８
ｆｉ
()
０．５
（５）
式中ｆｂ和ｆｉ分别为河床和冰盖边界的达西－威斯巴赫阻力系数。

在工程实践中，一般用曼宁糙率来表示河床和冰盖的阻力系数。

于是，将谢才公式与曼宁公式相结合，河床和冰盖的阻力系数ｆｂ和ｆｉ
可分别表示为：ｆｂ＝８ｇｎ２ｂＲ１?３ｂ，ｆｉ＝８ｇｎ２ｉＲ１?３
ｉ
（６）
式中：ｇ为重力加速度，ｍ?ｓ２
；ｎｂ和ｎｉ分别为河床底部和冰盖底面的曼宁糙率；Ｒｂ和Ｒｉ
分别为河床区和冰盖区的水力半径，ｍ。

Ｚａｒｅ等［３２］
通过将天然河道断面概化为矩形，给出了冰封河道水力半径的经验方程。

基于此，河床
区和冰盖区的水力半径可分别表示为：
Ｒｂ＝ＢｈｂＢ＋２ｈｂ，Ｒｉ
＝ＢｈｉＢ＋２ｈｉ
（７）
式中Ｂ为冰封河道宽度，ｍ。

于是，将式（５）—（７）代入式（４），可得河床区的水深计算式为：ｈｂ
＝Ｈ
１＋ｎｉｎｂｈｂＨ－ｈｂ·Ｂ＋２（Ｈ－ｈｂ
）Ｂ＋２ｈｂ
[]
１?６
（８）
从上述分析可知，当冰封河道的河床糙率ｎｂ和冰盖糙率ｎｉ
给定时，在式（８）基础上，利用迭代法可求得河床区水深ｈｂ，进而基于式（７）和式（６）计算河床区和冰盖区的水力半径和阻力系数，最终基于式（５）可获得参数ｍｂ和ｍｉ；对于流量参数Ｋ０，需根据经验进行率定。

于是，基于双幂律公式（２），便可预测冰封河道纵向时均流速Ｕ的垂向分布。

２．３　本文解析解　根据“双层假定”可知，河床区和冰盖区的切应力τｚｘ沿垂向均服从线性分布（见图１），两层交界面处切应力为零，河床底部和冰盖底面的切应力最大，则冰封河道的切应力τｚｘ
可表示为：τｚｘ＝τｂ－（τｂ＋τ
ｉ）ｚ
Ｈ
（９）
式中τｂ和τｉ
分别为河床区和冰盖区的边界切应力，Ｎ?ｍ２。

由于τｂ＝ρｕ２ ｂ且τｉ＝ρｕ２
ｉ
，则式（９）可化为：τｚｘ＝ρｕ２ ｂ
［１－（１＋λ２
）ξ］（１０）
式中：ρ为水体密度，ｋｇ?ｍ３
；λ＝ｕ ｉ?ｕ ｂ
；ξ＝ｚ?Ｈ。

令τｚｘ＝０，可得最大纵向时均流速的无量纲垂向位置ξｍａｘ＝１
?（１＋λ２
）。

于是，式（１０）可化为：τｚｘ＝ρｕ２
ｂ１
－ξξｍａｘ
()
（１１）
基于布辛涅斯克涡黏度模型，冰封河道切应力τｚｘ
又可表示为［３０］
：τｚｘ＝ρ（υ＋υｔ
）ｄＵ
ｄｚ
（１２）
式中：υ为运动黏度，ｍ２?ｓ；υｔ为水流紊动引起的涡黏度，ｍ２
?
ｓ。

考虑到天然冰封河道流动通常为紊流，黏性切应力ρυｄＵ?ｄｚ远小于紊流附加切应力ρυｔ
ｄＵ?ｄｚ。

于是，式（１２）可简化为：
τｚｘ＝ρυｔ
ｄＵ
ｄｚ
（１３）
联立式（１１）和式（１３），可得：
ｄＵｄｚ＝ｕ２
ｂ
υｔ１－ξξｍａｘ
()
（１４）
不难看出，求解式（１４）的关键在于涡黏度υｔ的确定。

鉴于河床底部和冰盖底面的粗糙度不同，冰
封河道湍流混合中存在两组速度和长度尺度，Ｇｕｏ等［３３］
将上述两组尺度扩展到一个涡黏度方程中，提
出一个适用于冰盖流涡黏度的解析模型：
υｔ＝２κＨｕ ｂ
βξ（１－ξ）１＋αξ?ξｃ－１()２
[]（１５）
式中：ξｃ为与涡黏度有关的参数，ξｃ＝１
１＋λｎ，其中ｎ＝５?６；α＝１－λλ－λ２ｎ；β＝λ－λ
２ｎ
２（１－λ２ｎ
）。

将式（
１５）代入式（１４），则有：ｄＵｄｚ＝ｕ ｂ２κＨβ·１－（ξ?ξｍａｘ）ξ（１－ξ）［１＋α（ξ?ξｃ－１
）２］（１６）
其中，ξ∈（０，１）。

由于ｄＵ?ｄｚ＝ｄＵ?（Ｈｄξ
），则式（１６）可化为：ｄＵｄξ＝ｕ ｂ
２κβ·１－ξ?ξｍａｘ()ξ（１－ξ）１＋αξ?ξｃ－１
()２[]（１７）
对方程（１７）进行积分，可得：
Ｕ＝ｕ ｂ２κβ∫
１－ξ?ξｍａｘ()ξ（１－ξ）１＋αξ?ξｃ
－１()２
[]ｄξ＝ｕ ｂ
２κβ·Ｉ（１８）
为了获得纵向时均流速Ｕ的解析解，将方程（１８）中的积分项Ｉ表示为部分分式，整理可得：
Ｉ＝∫ａ０
１ξ
－ｂ０
１１－ξ－ｃ０２α?ξｃ()ξξｃ－１()１＋αξ?ξｃ
－１()２＋ｄ０槡α?ξｃ１＋αξ?ξｃ－１()２{
}
ｄξ（１９）
式中：ａ０＝１１＋α，ｂ０＝１?ξｍａｘ－１１＋α（１?ξｃ－１）２，ｃ０＝（１＋α）?ξｍａｘ＋（α?ξｃ）（１?ξｃ－２
）２（１＋α）［１＋α（１?ξｃ－１）２
］，ｄ０
＝槡
α｛（１?ξｃ）［１＋α（１?ξｃ－１）］－（１?ξｍａｘ）［（１?ξｃ－１）（１＋α）］｝（１＋α）［１＋α（１?ξｃ－１
）２
］。

将式（１９）代入式（１８），积分可得：
Ｕ＝ｕ ｂ２κβ
｛ａ０ｌｎξ＋ｂ０ｌｎ（１－ξ）－ｃ０ｌｎ［１＋α（ξ?ξｃ－１）２
］＋ｄ０ａｒｃｔａｎ［槡α（ξ?ξｃ－１）］｝＋Ｃ（２０）
式中Ｃ为待定积分常数。

由于ξ＝ξｍａｘ时，Ｕ＝Ｕｍａｘ
，则式（２０）可转化为：Ｕ＝Ｕｍａｘ－ｕ ｂ２κβａ０ｌｎξｍａｘξ＋ｂ０ｌｎ１－ξｍａｘ１－ξ－ｃ０ｌｎ１＋α（ξｍａｘ?ξｃ－１
）２
１＋α（ξ?ξｃ－１
）２
{
＋ｄ０
ａｒｃｔａｎ（槡
α?ξｃ）（ξｍａｘ－ξ）１＋α（ξｍａｘ?ξｃ－１）（ξ?ξｃ－１
）}
（２１）
公式（２１）即为冰封河道纵向时均流速垂向分布的解析解，其中，ξ∈０，１()。

在应用式（２１）预测冰封河道纵向时均流速垂向分布之前，首先需要确定式（２１）中的物理参数，具体步骤为：最大纵向时均流速Ｕｍａｘ由纵向时均流速Ｕ的垂向分布实测数据确定；参数ξｍａｘ由式（８）及ξｍａｘ＝ｈｂ?Ｈ确定；进而，由ξｍａｘ＝１?（１＋λ２）确定参数λ，由ξｃ＝１?（１＋λｎ
）确定ξｃ，由公式α＝１－λλ－λ
２ｎ
和β＝λ－λ
２ｎ
２（１－λ
２ｎ
）分别确定参数α和β；对于河床区摩阻流速ｕ ｂ，基于式（２１）及实测流速数据，应用最小二乘法确定。

通过上述步骤，便可求得解析解（２１）中的待定物理参数，进而可计算冰封河道中纵向时均流速的垂向分布。

３　分析与讨论
为了验证本文解析解（２１）的可靠性和精确性，基于以往实测流速数据，将其流速计算结果与对数公式和双幂律公式的计算结果进行对比分析。

所采用的冰封河道水槽实验和野外原型实验的流速数据共包括１２种工况，其中水槽实验涵盖了光滑冰盖和粗糙冰盖且均是在均匀流状态下开展的，主要通过布置测针来监测水槽实验段内水深是否沿程发生变化来判定水流是否达到均匀流状态。

为方便起
见，Ｌａｕ和Ｋｒｉｓｈｎａｐｐａｎ［２７］的实验工况简记为５Ｃ、６Ｃ和７Ｃ，Ｓａｙｒｅ和Ｓｏｎｇ［２９］
的实验工况简记为ＡＳ、ＡＲ、ＢＳ和ＢＲ，魏良琰和黄继忠［３４］的实验工况简记为ＷＳ和ＷＲ，Ｗａｎｇ等［３５］的实验工况简记为ＷＭ，罗红春［３６］的原型实验工况简记为ＬＭ和ＬＲ。

下面针对上述１２种工况的实验方案进行简单描述。

Ｌａｕ和Ｋｒｉｓｈｎａｐｐａｎ［２７］
在长５７．３ｍ、宽０．７５６ｍ、深０．２９ｍ的矩形水槽中模拟冰封河道开展水流特性实验，槽底铺设中值粒径为０．１５ｍｍ的玻璃珠，底坡Ｓ０保持０．００１不变，流量变化范围为０．０１５４ｍ３
?ｓ～０．０２７４ｍ３?ｓ，冰盖由１．９ｃｍ厚的胶合板模拟，其底面黏合层压塑料以形成光滑冰盖；Ｓａｙｒｅ和Ｓｏｎｇ
［２９］在具有玻璃侧壁的循环水槽中模拟冰封河道并开展实验，测量段长２７．４ｍ、宽０．９１４ｍ、深０．４５ｍ，槽底铺设中值粒径为０．２５ｍｍ的石英砂，底坡可调节且不干扰水槽运行，冰盖由光滑冰盖和粗糙冰盖
两种类型组成，光滑冰盖由表面涂有油漆的胶合板模拟，粗糙冰盖由底面加肋的胶合板模拟；魏良琰
和黄继忠［３４］
在长３２ｍ、宽０．５ｍ、深０．５ｍ的固定玻璃水槽中开展实验，槽底为玻璃平板或铺以８ｃｍ碎石层，模拟冰盖为聚塑泡沫板或加肋聚塑泡沫板；所采用的Ｗａｎｇ等［３５］全冰封复式河道的主槽流速
实验数据是在长２０ｍ、宽１ｍ、深０．５ｍ的顺直玻璃水槽中测量的，水槽底部采用ＰＶＣ板以构成对称复式断面，其由一个主槽和两个河漫滩组成，河漫滩和主槽之间由１
∶１的斜边坡连接，主槽宽０．４ｍ，满槽高度为０
．１ｍ，河漫滩宽度为０．２ｍ，冰盖利用高密度ＥＰＳ泡沫板模拟；所采用的罗红春［３６］
原型实验流速数据是在黄河内蒙古段下游的什四份子弯道完整断面ｓｅｃ．２的Ｍ侧和Ｒ侧测量的，Ｍ侧靠近凹岸且有冰塞堆积，而Ｒ侧靠近凸岸且流速较大。

各实验工况的流动参数如表１所示。

表１　冰封河道水槽实验流动条件以及相关参数设置
来源
工况ｑ?（ｍ２
?ｓ）
Ｈ?ｍＢ?ｍＳ０ρ
?（ｋｇ?ｍ３）ＲｅｎｉｎｂＵｍａｘ?（ｍ?ｓ）Ｋ０?（ｍ?ｓ）Ｌａｕ和Ｋｒｉｓｈｎａｐｐａｎ
［２７］
５Ｃ
０．０２６６０．０９７６０．７５６０１．００×１０－３
１．００×１０３
１１６９７０．００９００．０２２５０．３２９３０．４２０６Ｃ０．０３３３０．１１１６０．７５６０１．００×１０－３１．００×１０３１４４０７０．００９００．０２２５０．３６７１０．４５５７Ｃ０．０３６２０．１１６００．７５６０１．００×１０－３１．００×１０３１５５８３０．００９００．０２２５０．３８２６０．４７５Ｓａｙｒｅ和Ｓｏｎｇ
［２９］
ＡＳ
０．０３５５０．１１８００．９１４０１．８０×１０－３１．００×１０３１４４２１０．０１１４０．０２９４０．３８１１０．５１５ＡＲ０．０３３２０．１２１００．９１４０１．８０×１０－３１．００×１０３１４７１８０．０２４９０．０２２８０．３２０６０．４８５ＢＳ０．０３９２０．１２２００．９１４０２．０４×１０－３１．００×１０３１７０３６０．０１０８０．０３０６０．４２１９０．５６０ＢＲ
０．０４２８０．１４８００．９１４０２．１１×１０－３１．００×１０３１７４５４０．０３２２０．０２３６０．３６０４０．５６５魏良琰和黄继忠［
３４］
ＷＳ０．１０１４０．２１５４０．５００００．６６×１０－３１．００×１０３４０３４００．０１０７０．００８８０．５４９６０．６４５ＷＲ０．１０１４０．２３７００．５０００２．５７×１０－３１．００×１０３３９１５８０．０２８８０．０２４７０．５８２５０．８７５Ｗａｎｇ
等［３５］
ＷＭ０．１０２００．３０００１．００００１．００×１０
－４１．００×１０３４５８３２０．０１０５０．００９００．２２７４０．２８５罗红春
［３６］
ＬＭ０．８６００２．５０００６００１．００×１０－４
１．００×１０３５４６５１８０．０２２００．０２９００．２７８７０．３６０ＬＲ
０．８６００
２．８８００
６００
１．００×１０
－４１．００×１０
３５４６５１８０．０１７５０．０２９０
０．４７０６
０．５９８
３．１　计算结果分析　基于上述１２种实验工况，图２给出了本文纵向时均流速解析解、对数公式、双幂律公式的计算值与实测值的比较，其中Ｕｍ为断面平均流速。

由图２可知，三种公式总体上均能较好预测纵向时均流速的垂向分布。

当冰盖底面糙率ｎｉ小于河床底部糙率ｎｂ时（例如，工况５Ｃ、６Ｃ、７Ｃ、ＡＳ、ＢＳ和ＬＲ），最大纵向时均流速Ｕｍａｘ偏向于冰盖区；相似的，当河床底部糙率ｎｂ小于冰盖底面糙率ｎｉ时（例如，工况ＢＲ和ＷＲ），最大流速Ｕｍａｘ偏向于河床区；当河床底部糙率ｎｂ与冰盖底面糙率ｎｉ相接近时（例如，工况ＡＲ、ＷＳ、ＷＭ和ＬＭ），最大流速Ｕｍａｘ位于水深中部，纵向时均流速Ｕ近似呈轴对称分布。

不难得出，上下边界糙率对纵向时均流速垂向分布具有不可忽视的影响。

对比图２（ａ）—（ｃ）可以发现，本文解析解相较于对数公式和双幂律公式，能较好预测河床区的纵向时均流速垂向分布；由图２
（ｄ）—（ｆ）又可发现，本文解析解对于冰盖区纵向时均流速的垂向分布具有较高预测精度。

上述结果证实了本文解析解预测冰封河道纵向时均流速垂向分布的可靠性和有效性。

为了量化本文纵向时均流速解析解与对数公式、双幂律公式的预测精度并阐明计算结果与实测结果的相关性，这里采用平均相对误差（
ＭＲＥ）和皮尔逊相关系数（ＣＯＲ）来说明，如表２所示。

平均相对误差（ＭＲＥ）定义为计算值和实测值之差与实测值比值的平均值：
ＭＲＥ＝１Ｎ∑Ｎｉ＝１Ｕｃｏｍ，ｉ－Ｕｍｅａ，ｉ
Ｕｍｅａ，ｉ
×１００％
（２２）
皮尔逊相关系数（ＣＯＲ）定义为：
ＣＯＲ＝
∑Ｎ
ｉ＝１Ｕｃｏｍ，ｉＵｍｅａ，ｉ－１Ｎ∑Ｎ
ｉ＝１Ｕｃｏｍ，ｉ∑Ｎ
ｉ
＝１Ｕｍｅａ，ｉ∑Ｎ
ｉ
＝１Ｕ２
ｃｏｍ，ｉ－１
Ｎ∑
Ｎ
ｉ＝１Ｕｃｏｍ，ｉ()２()∑Ｎ
ｉ
＝１Ｕ２
ｍｅａ，ｉ－１
Ｎ∑
Ｎ
ｉ＝１Ｕｍｅａ，ｉ()２(
)
槡（２３）
式中：Ｎ为各工况沿垂向所布置的流速测点数量；Ｕｍｅａ，ｉ为测点ｉ处纵向时均流速的实测值；Ｕｃｏｍ，ｉ表示与实测值Ｕｍｅａ，ｉ
相对应的计算值。

图２　纵向时均流速解析解与对数公式、双幂律公式计算值和实测值的比较表２　不同流速垂向分布公式的预测精度比较
工况
本文解析解对数公式双幂律公式
ＭＲＥ?％ＣＯＲＭＲＥ?％ＣＯＲＭＲＥ?％ＣＯＲ
５Ｃ１．８８０．９８６６４．５１０．９５２７４．２７０．９９３６６Ｃ３．６７０．９８２１８．８００．９２２９８．４１０．９７６４７Ｃ４．１６０．９７９５１０．２５０．９２７９１２．７７０．６９１６ＡＳ４．３７０．９８０９１０．８８０．９７４６７．２８０．９８３６ＡＲ２．７３０．９６７０６．０９０．９７６４２．５４０．９６６５ＢＳ４．２４０．９８１９６．７２０．９４６９４．２３０．９７８７ＢＲ２．０１０．９８７８６．４４０．９８１６２．０４０．９９１０ＷＳ４．２７０．９１６６６．７９０．８８１３３．７１０．９３５９ＷＲ１３．１５０．９９１９３１．３６０．９９２３４９．１１０．９８９６ＷＭ２．０１０．９９８４６．６４０．９８８５２．９８０．９９７０ＬＭ１．４９０．９９９２２．１６０．９９８３３．９６０．９９２２ＬＲ２．０６０．９９７５２．１３０．９９５８３．１８０．９８６７
由表２可知，本文解析解的计算误差ＭＲＥ均在１０％以下（工况ＷＲ除外），相关系数ＣＯＲ均在０．９１以上；然而，对数流速公式关于工况７Ｃ、ＡＳ、ＷＲ的计算误差均在１０％以上，双幂律流速公式关于工况ＷＲ的计算误差ＭＲＥ高达４９．１１％且工况７Ｃ的相关系数低于０．７０。

不难发现，本文基于冰盖流涡黏度模型提出的纵向时均流速垂向分布解析解相较于常用的对数流速公式和双幂律流速公式具有较高准确度和可靠性。

３．２　讨论　基于所提出的纵向时均流速解析解，图３以工况６Ｃ、ＡＲ和ＷＲ为例，探讨了单宽流量ｑ、水深Ｈ和边界糙率ｎ对冰盖下纵向时均流速垂向分布的影响，其中Ｒｅ的变化范围为７×１０３～１．４５×１０５。

在其他参数保持不变的情况下，图３（ａ）展示了单宽流量ｑ对纵向时均流速垂向分布的影响，可以看出流速垂向分布并不随流量的改变而变化，意味着单宽流量对流速垂向分布的影响可忽略不计；由图３（ｂ）可知，流速垂向分布对水深的变化相对于流量的变化比较敏感，随着水深Ｈ的增加，流速Ｕ?Ｕ
ｍ随之增大，但无论水深如何变化，并不会改变最大纵向时均流速的垂向位置，这主要是由于最大流速位置是由冰盖底面糙率和河床底部糙率的相对大小决定的，图３（ｃ）刚好证实了这一结论。

由图３（ｃ）可
的增大，最大纵向时均流速逐渐偏向于河床区一知，无论是工况６Ｃ、ＡＲ还是ＷＲ，随着冰盖糙率ｎ
ｉ
侧，即最大纵向时均流速靠近相对光滑边界一侧；对比图３（ｂ）与图３（ｃ）可知，当水深变化幅度较大时，其对纵向时均流速垂向分布的影响才较为明显，但冰盖糙率和河床糙率相对值的较小波动，便会
Ｒ?υ＝ｑ?［２υ（１＋Ｈ?Ｂ）］，单宽流量ｑ和水深Ｈ的变引起流速垂向分布的较大变化。

考虑到雷诺数Ｒｅ＝Ｕ
ｍ
化反映了雷诺数变化。

因此，流速垂向分布对单宽流量和水深变化不敏感意味着雷诺数对冰封河道流速垂向分布的影响并不显著，影响流速垂向分布的主导因素是河床和冰盖边界糙率的相对大小，这与王志兴等［２５］基于三维ｋ－ε湍流模型得到的研究结果是一致的。

图３　单宽流量ｑ、水深Ｈ、边界糙率ｎ对纵向时均流速垂向分布的影响
值得注意的是，虽然对数公式也能较好预测冰盖流纵向时均流速的垂向分布，但流速分布曲线在最大流速位置处是不光滑的（见图２），即在河床区与冰盖区交界面处不满足流速梯度连续，违背了水力学知识中的液体连续介质假定。

这一结果主要是将明渠流Ｐ
ｒａｎｄｔｌ对数流速公式直接拓展到冰盖流中造成的。

下面，从理论角度证实对数流速公式梯度不连续这一结论。

通过对式（１）求导，可得：
河床区（０＜ｚ≤ｈｂ
）：ｄＵｂｄｚ＝ｕ ｂ
κｚ
（２４）
冰盖区（ｈｂ
≤ｚ＜Ｈ）：ｄＵｉｄｚ＝－
ｕ ｉκ（Ｈ－ｚ）
（２５）由式（２４）（２５）易知，河床区的流速梯度ｄＵｂ?ｄｚ具有正值，而冰盖区的流速梯度ｄＵｂ?ｄｚ具有负值。

因此，在交界面ｚ＝ｈｂ
处，ｄＵｂｄｚ
ｚ＝ｈｂ
≠
ｄＵｉ
ｄｚｚ＝ｈｂ。

这一结果表明对数流速分布曲线在河床区和冰盖区交界
面处难以满足连续光滑条件。

虽然双幂律流速公式克服了对数流速公式在最大流速位置处梯度不连续的缺点，但该公式未明确涉及冰盖区与河床区的速度和长度尺度，物理意义不直观，并且公式中的Ｋ０值是根据经验进行率定的，其通用表达式尚需进一步探讨。

基于本文所率定的各工况下的Ｋ０值（见表１），图４展示了参数Ｋ０随单宽流量ｑ的变化规律，并给出了参数Ｋ０与单宽流量ｑ之间的拟合关系式。

不难看出，Ｋ０
值随单宽流量ｑ的增大而增加，即参数Ｋ０
与单宽流量ｑ呈正相关关系，该结果与Ｂａｉ和Ｄｕａｎ［３７］
所得结论是一致的。

拟合结果表明，参数Ｋ０
与单宽流量ｑ存在如下关系式：Ｋ０＝０
．２８０２ｌｎｑ＋１．２２ｑ＋１．３９３（２６）
图４　参数Ｋ０
随单宽流量ｑ的变化但由于所采用的实验数据有限，上述公式的普适性仍需基于大量数据进行验证。

本文解析解（即式（２１））是在冰盖流涡黏度模型基础上导出的，其考虑了冰盖流中因河床底部和冰盖底面粗糙度不同而存在的两组速度和长度尺度，并且较对数流速公式和双幂律流速公式具有较高精度和可靠性。

值得注意的是，在推导本文解析解过程中，由于假设黏性底层的黏滞切应力可忽略，而只考虑了脉动流速产生的紊流附加切应力，则对于雷诺数较小的冰封渠道流动型态（如层流、层流到紊流的过渡区、紊流的水力光滑区及水力粗糙过渡区），其黏滞切应力对水流的影响不可忽略，此时解析解将不能用于计算纵向时均流速的垂向分布。

但在实际工程中，天然冰封河道通常具有较高雷诺
数且流动型态为紊流，则应用本文解析解预测冰盖下纵向时均流速垂向分布时，只需要确定实测最大纵向时均流速Ｕｍａｘ和河床区摩阻流速ｕ ｂ。

具体地，首先可利用式（８）计算河床区水深ｈｂ，即得到最大纵向时均流速Ｕｍａｘ的垂向位置，进而测量冰封河道该位置处流速以获得Ｕｍａｘ
；其次，应用水文站常用的定点测流法获得特定位置处纵向时均流速的实测值［３４］，并基于最小二乘法得到摩阻流速ｕ ｂ
；最后，应用公式（２１）便可计算得到冰封河道纵向时均流速Ｕ的垂向分布。

４　结论
本文基于“双层假定”断面上切应力沿垂向服从线性分布以及布辛涅斯克涡黏度模型，推导出冰盖流纵向时均流速垂向分布的解析解。

进而，将该解析解与常用的对数流速公式和双幂律流速公式的计算结果对比分析，所得结论如下：
（
１）所提出的解析解能较好预测冰封河道纵向时均流速的垂向分布，流速分布曲线在最大纵向时流速位置处满足连续光滑条件，并且在河床和冰盖边界糙率影响下，最大纵向时均流速总是偏向于光滑边界一侧，在上下边界糙率相接近情况下，最大纵向时均流速位于水深中部，流速垂向分布呈轴对称分布，这与常用的对数公式和双幂律公式所预测的流速分布曲线相一致。

基于误差统计分析，证实了所提出的纵向流速解析解的可靠性，并且相较于对数流速公式和双幂律流速公式具有较高精度。

（２）敏感性分析表明，河床和冰盖糙率相对值的较小改变便会引起纵向时均流速垂向分布的较大变化，并且最大纵向时均流速位置由上下边界底部的相对糙率所决定，结果与先前学者基于三维ｋ－ε湍流模型所得结果是一致的。

（
３）从理论角度阐明了对数公式所预测的纵向时均流速分布曲线在河床区和冰盖区交界面处不光滑的原因，指出了双幂律公式因需率定参数Ｋ０而存在的局限性，明确了如何在工程实际中应用本文所提出的解析解对冰封河道纵向时均流速的垂向分布进行预测。

本文虽然给出了参数Ｋ０与单宽流量ｑ之间的拟合关系式，但其普适性仍需进一步验证。

参　考　文　献：
［１］　樊霖，茅泽育，吴剑疆，等．松花江白山河段冰情数值计算及分析［Ｊ］．水科学进展，２０１６，２７（６）：
８９０－８９７．
［２］　宋春山，朱新宇，韩红卫，等．河道特征对黑龙江上游冰坝生消影响［Ｊ］．水利学报，２０２０，５１（１０）：
１２５６－１２６６．
［３］　ＫＯＹＵＮＣＵＢ，ＬＥＴＢ．Ｏｎｔｈｅｉｍｐａｃｔｓｏｆｉｃｅｃｏｖｅｒｏｎｆｌｏｗｐｒｏｆｉｌｅｓｉｎａｂｅｎｄ［Ｊ］．ＷａｔｅｒＲｅｓｏｕｒｃｅｓＲｅｓｅａｒｃｈ，
２０２２，５８（９）：ｅ２０２１ＷＲ０３１７４２．
［４］　ＷＡＮＧＦＦ，ＨＵＡＩＷＸ，ＧＵＯＹＫ，ｅｔａｌ．Ｔｕｒｂｕｌｅｎｃｅｓｔｒｕｃｔｕｒｅａｎｄｍｏｍｅｎｔｕｍｅｘｃｈａｎｇｅｉｎｃｏｍｐｏｕｎｄｃｈａｎｎｅｌ
ｆｌｏｗｓｗｉｔｈｓｈｏｒｅｉｃｅｃｏｖｅｒｅｄｏｎｔｈｅｆｌｏｏｄｐｌａｉｎｓ［Ｊ］．ＷａｔｅｒＲｅｓｏｕｒｃｅｓＲｅｓｅａｒｃｈ，２０２１，５７（４）：ｅ２０２０ＷＲ０２８６２１．
［５］　郭新蕾，王涛，付辉，等．河渠冰水力学研究进展和趋势［Ｊ］．力学学报，２０２１，５３（３）：６５５－６７１．［６］　翟必篧，张宝森，刘璐，等．冰盖前缘冰坝形成及演进过程的离散元－水动力学耦合分析［Ｊ］．水动力学研
究与进展（Ａ辑），２０２１，３６（５）：６５７－６６７．
［７］　王涛，余弘婧，郭新蕾，等．黄河内蒙古河段典型凌情过程模拟［Ｊ］．水利学报，２０２１，５２（５）：５３５－５４５．［８］　王军，苏奕垒，侯智星，等．冰盖条件下桥墩局部冲刷研究进展［Ｊ］．水利学报，２０２０，５１（１０）：１２４８－１２５５．［９］　秦子鹏，宗全利，田艳，等．寒区引水明渠抽水融冰试验［Ｊ］．武汉大学学报（工学版），２０１７，５０（２）：
１６８－１７３．
［１０］　王军，苏奕垒，侯智星，等．冰盖条件下桥墩局部冲刷研究进展［Ｊ］．水利学报，２０２０，５１（１０）：１２４８－１２５５．［１１］　ＧＥＢＲＥＳ，ＢＯＩＳＳＹＴ，ＡＬＦＲＥＤＳＥＮＫ．Ｓｅｎｓｉｔｉｖｉｔｙｔｏｃｌｉｍａｔｅｃｈａｎｇｅｏｆｔｈｅｔｈｅｒｍａｌｓｔｒｕｃｔｕｒｅａｎｄｉｃｅｃｏｖｅｒｒｅｇｉｍｅ
ｏｆｔｈｒｅｅｈｙｄｒｏｐｏｗｅｒｒｅｓｅｒｖｏｉｒｓ［Ｊ］．ＪｏｕｒｎａｌｏｆＨｙｄｒｏｌｏｇｙ，２０１４，５１０：２０８－２２７．。