适用于明渠流及冰盖流的统一泥沙起动流速计算方法初探

刘晓民,杨耀天,张红武,王文娟

(1.内蒙古农业大学 水利与土木建筑工程学院,呼和浩特 010018; 2.黄河流域内蒙段水资源与水环境综合治理协同创新中心,呼和浩特 010018; 3.清华大学 黄河研究中心,北京 100084;4.内蒙古金华源环境资源工程咨询有限责任公司,呼和浩特 010020)

1 研究背景

冰凌是一种自然现象,主要发生在我国黄河流域和东北地区等北方高寒地区。黄河凌汛已成为我国冬春季节大江大河中最突出、最主要的汛情之一[1],尤其黄河流域内蒙古段是冰凌洪水经常发生的区域(见图1)。河冰的出现改变了河流的水力条件、热力条件和边界几何条件等[2],冰盖或冰塞的形成,改变了河流的运动状态,由原来的明流变成了覆盖流,过水断面湿周和阻力明显增加,断面平均流速大幅减小,导致水沙输移特征发生变化,产生了明显异于明渠流的特有泥沙输移现象。

图1 黄河内蒙古段位置Fig.1 Location of the Yellow River in Inner Mongolia reach

目前,明渠流条件下的泥沙起动流速公式已有很多研究[3],但是,关于冰盖流泥沙输移的研究成果并不多见。同明渠相比,冬季河道中冰盖的存在影响了水流的水力特性,比如流速分布、剪切应力分布、泥沙浓度等,进而也改变了原有水沙输移特性及规律。Turcotte等[4]讨论了在河冰形成过程中不同形态的冰对泥沙输移产生的影响;马子普等[5]基于Einstein假定[6],采用指数形式的流速分布公式,推导得到了适用于明渠流及冰盖流的统一的非黏性泥沙颗粒起动流速公式;王军等[7]、陈建国等[8]在室内水槽试验中模拟探讨了冰盖流输沙特性及受力规律;全栋等[9]、宋本辉等[10]利用头道拐水文站资料与野外观测采样相结合,探究了冰下泥沙输移特性及泥沙冲淤规律。

从上述认识出发,本文基于Einstein假定[6],在前人研究的基础上,考虑了冰盖流下冰盖糙率变化对泥沙起动流速的影响,采用泥沙滚动起动模型,运用概率论与力学相结合的分析方法,初步提出了冰盖流泥沙起动流速公式,并探讨了该式在明渠流条件下泥沙起动的适用性。

2 公式推导

2.1 冰盖流泥沙起动标准探讨

一般学者多在力学分析的基础上,结合概率统计理论等来推导泥沙起动的临界条件。目前,对于明渠水流泥沙起动临界条件的研究,公认的是水流底部的瞬时流速u0近似服从正态分布[11]。本文基于Einstein阻力划分原则,即冰盖流流速分布以最大流速线为界将流动沿水深方向分为冰盖区和河床区2层(如图2),并将2层分别等效成明流对数分布[12]。假定冰盖流和明流造床机理相同,故在探讨分析冰盖流泥沙起动过程中,可认为冰盖流河床区与明渠泥沙起动瞬时流速分布理论相似。

注:H0为冰盖下流动的总水深,Yb为河床区水深,Yi为冰盖区水深,nb为河床糙率,ni为冰盖糙率,u为时均流速。图2 冰盖流沿水深方向流速分布Fig.2 Velocity of ice-covered flow along the depth direction

明渠水流底部的瞬时流速分布函数f(u0)为

(1)

类似的,冰盖下水流河床区的瞬时流速分布函数f(ub)为

(2)

同时,窦国仁[11]提出了瞬时底速与平均底速的关系,即

(3)

式中β为与起动标准相应的参数,其中β为3、2、1、0分别表示弱动、中动、普动、强动。

目前,明渠流中大多以弱动(床面上有很少的泥沙在运动,β=3)作为泥沙起动的临界判别标准,基于此,本文以β=3作为判别冰盖流泥沙起动临界标准的临界条件。

2.2 受力分析与起动模式

2.2.1 受力分析

无论是明渠流还是冰盖流,床面上泥沙颗粒起动所受力的种类并无差别。如图3所示,在水流作用下,粗颗粒泥沙起动时一般受力有水流拖曳力FD、上举力FL及水下重力G,计算公式分别为:

注:Fμ为黏结力,FN为附加下压力,a为水流拖曳力FD与泥沙颗粒中心之间的垂直距离,b为水流上举力FL与泥沙颗粒中心之间的水平距离。图3 泥沙颗粒起动受力分析Fig.3 Analysis of forces on incipient sediment particles

(4)

(5)

G=αw(γS-γ)D3。

(6)

式中:αL、αD分别为与水流拖曳力、上举力相应的面积系数,一般取值π/4;αw为与水下重力相应的体积系数,取值π/6;CD为水流拖曳力系数,取值0.4;CL为水流上举力系数,取值0.1[13-14];ρ为水的密度;γS、γ分别为泥沙颗粒及水的重度;D为泥沙颗粒的粒径;g为重力加速度。无论是明渠流还是冰盖流,达到泥沙起动的瞬时流速是相同的。

对于细颗粒泥沙,除受到以上3种力的作用外,还受到颗粒间黏结力Fμ、附加下压力FN的作用[15],如图3所示。FN的表达式为

(7)

此式在以往附加下压力的表达式中引入冰盖糙率与河床糙率的比值n来修正冰凌期不同时段冰盖糙率变化对某一粒径泥沙颗粒附加下压力的影响。在Calkins等[16]构建的糙率比值与水深比值关系的基础上,考虑到实际河道中冰盖厚度相对于水深影响较小,将n表达式简化为

(8)

Fμ的表达式为

式中:αN、αμ分别为与附加下压力与黏结力相应的面积系数;γ′为淤积物的干重度;γ′c为淤积物的稳定干重度,γ′/γ′c≈0.85[3];κ*为黏结力影响系数;δ为薄膜水厚度,按照窦国仁[11]交叉石英丝试验成果取为0.213×10-4cm;υ为水流运动黏滞性系数,与温度变化相关。

基于对多年水温和气温实测数据的分析,就黄河冰封内蒙古段而言,干流各段及其支流的年平均水温都略高于当地的年平均气温,但差值不大,一般为1～2 ℃,冬季多年平均水温和气温差值较大,约在2 ℃以上[17]。黄河内蒙古段畅流期平均气温7.3～7.9 ℃,则水温基本为8.3～9.9 ℃,本次简化取均值9.1 ℃,查水的运动黏滞系数取值表得υ=0.013 cm2/s。同时,根据黄河内蒙古段稳定封冻区域4个主要水文站的实测资料,当旬平均气温<-5 ℃时,开始封河;当旬平均气温<-2 ℃时,开始流凌[18]。基于冀鸿兰[17]的研究,封冻期黄河内蒙古段的水温接近-3～0 ℃,通过温度和黏滞系数的关系计算式[19],可得黄河内蒙古段冰封期黏滞系数υ的取值范围为0.017 7 ～0.019 7 cm2/s,本文简化取均值0.018 7 cm2/s。

2.2.2 起动模式

本文起动模式采用滚动模式,泥沙颗粒绕某支点滚动的平衡方程为

FLLL+FDLD=GLw+FμLμ+FNLN。

(10)

式中LD、LL、Lw、Lμ、LN分别为与水流拖曳力、上举力、水下重力、黏结力及附加下压力相应的力臂。

本文参考张红武[3]对于附加下压力与黏结力作用力臂和面积系数的处理方法,令awLw=C,假定C1、C2为附加下压力及黏结力面积系数和力臂的乘积与重力G的修正系数,代入式(10)得

FLLL+FDLD=αwLw(G+C1FN+C2κ*Fμ) 。

(11)

式中:C1=0.035 2;κ*C2=6.59。

本文研究的水流对床面泥沙作用力及作用力的力臂大小是通过床面泥沙颗粒间的相对暴露度来反映的。韩其为[20]首先提出了相对暴露度Δ′(Δ′=Δ/D)的概念,其中Δ为绝对暴露度,用床面泥沙颗粒的最低点至两颗粒接触点的距离来表示,如图3所示。相对暴露度取值范围的理论研究目前尚无定论,不同学者有着不同的观点,周双等[21]通过试验观察,发现相对暴露度取值范围为0～0.5,与杨奉广等[22]对于相对暴露度的取值范围一致。基于此,考虑到河床泥沙颗粒分布与相对暴露度取值的随机性,为简化计算复杂程度,假定相对暴露度符合均匀分布[21],采用相对暴露度均值代表河床床面泥沙颗粒分布的平均情况,即Δ′=0.25,并以相对暴露度Δ′均值对应的位置关系进行泥沙起动力臂计算。

故绕实际滚动支点A(如图3所示),将作用力力臂表达式代入式(11)得

式中:a为正面推移力至泥沙颗粒中心的垂直距离,参考文献[23]取值,a=D/10;b为上举力至泥沙颗粒中心的水平距离,根据韩其为[20]的研究,b=D/6。

将式(4)—式(7)、式(9)代入式(12),整理得冰盖流泥沙瞬时流速公式,即

(13)

2.3 构建冰盖流泥沙起动公式

冰盖流垂线流速分布可近似等同于紊流粗糙管垂线流速分布,紊流粗糙管流速常用对数分布公式表示,即

(14)

类似的,冰盖流河床的过水断面垂向流速为

(15)

式中:u*为冰盖下摩阻流速(m/s);y0、yb分别为明渠流及冰盖流距离底部的垂向距离(m);κs为渠壁粗糙度(m)。

一般在床面粗糙区,可认为y0=κs=D,又因为冰盖下水流分布与紊流粗糙区(阻力平方区)相似,故可简化认为yb=κs=D。则在床面颗粒的顶端,冰盖流的时均流速与摩阻流速的关系可简化为:

(16)

(17)

式中τ为冰盖下水流切应力。

考虑水流作用在颗粒上的作用点,取底流速作用点的位置为(2/3)D[14],则作用于泥沙颗粒上的时均底流速可表示为

(18)

基于韩其为等[15]的研究, 明渠垂线平均流速ν0与u*之间的关系为

(19)

根据已有研究成果,冰盖流垂线流速分布随冰盖糙率改变而改变。王军等[7]通过室内试验,以中垂线平均流速来表示泥沙起动流速,构建了冰盖流与明流泥沙起动流速的比值与冰盖与河床糙率比值的函数表达式,即

(20)

当ni=0时,冰盖流泥沙起动流速即化为明渠泥沙起动流速。

vb=ξv0。

(21)

类似的,冰盖流垂线平均流速vb与u*的关系可表示为

(22)

由式(3)、式(16)、式(19)、式(22)代入式(13)整理得冰盖流泥沙起动垂线平均流速为

(23)

考虑到天然泥沙颗粒形状对起动产生的影响,引入韩其为等[15]提出的扁度函数f(λ)对上述公式进行修正,其中λ为扁度系数,则底部起动作用流速为

(24)

参考韩其为等[15]给出的f(λ)取值表,本文取f(λ)=1.1,则

(25)

冰盖的存在,改变了垂线上的流速分布。故当泥沙颗粒起动流速以垂线平均流速来表示时,相应会改变。本文公式在推导过程中,考虑了冰盖流糙率变化对垂线流速的影响,及泥沙颗粒间的相对暴露度与细颗粒泥沙间黏结力和附加下压力对泥沙起动的影响,初步探讨提出了适用于明渠流及冰盖流泥沙起动流速公式。

3 合理性分析

王军等[7]通过在泡沫板中插入竹签的方法来率定不同阶段冰盖糙率,进行了冰盖下散粒体泥沙起动流速的试验研究,其中冰盖有光滑冰盖(ni=0.021 2)、稀疏冰盖(ni=0.032 2)、密集冰盖(ni=0.034 7)3种不同类型,相对应通过在河床铺设不同粒径的泥沙颗粒,率定河床糙率依次为0.010 9(D50=0.32 mm)、0.012 8(D50=0.85 mm)、0.013 8(D50=1.32 mm)。

本文根据王军等[7]试验所得数据,将本文公式与马子普公式[5](图4中简写为马式)的计算值绘入图4中。并通过计算泥沙起动流速相对误差,对本文公式的合理性进行检验,见表1。表1给出了本文公式、马子普公式[6]与王军等[7]试验资料起动流速对比分析结果,可以看出本文公式计算值与试验资料的吻合度更好,平均相对误差均<10%,说明本文所引用的部分经验公式和经验取值是合适的。另外,随着河床泥沙颗粒粒径的增大,马子普公式计算值与试验值平均相对误差相比于本文公式更小,这说明本文公式的适用性还有待进一步检验,且马子普公式[6]更适用于粗颗粒泥沙。

表1 冰盖流泥沙起动流速不同公式计算值与试验值的平均相对误差Table 1 Average relative error between the calculated values of different formula and the test values for incipient velocity of sediment under ice-covered flow

图4 冰盖流起动流速公式计算值与试验值对比及相对误差分析Fig.4 Comparison between calculated values and observed data of incipient velocity of sediment under ice-covered flow and analysis of relative error

图4显示,本文公式和马子普公式[6]的计算值与王军等[7]试验值整体上符合良好。本文公式的相对误差均在0%～22%之间,且大多数集中在1%～10%之间,相对误差较小,满足计算精度要求。

与试验实测值相比,本文公式计算值在小水深时较实测值偏大,而大水深时偏小。初步分析原因可能是冰封初期,冰盖形成的结构复杂,冰盖下方的水流常常伴随着冰花运动[24]。由于目前尚不能表述清楚冰花的生长过程,若冰花颗粒被带入紊动水流中,会形成浮冰,波动水流中浮冰之间相互摩擦,改变水流运动参数和流速分布,而本文公式在推导过程中用中垂线平均流速来表征泥沙运动状态,与水深变化关联密切,不同水深下河床近底流速大小不同,而近底水流流速是决定泥沙起动的主要水流参数,具体改进方法还有待在此基础上进一步研究。

同时,限于王军等[7]的试验均是粗颗粒泥沙在冰盖糙率ni均大于河床糙率nb情况下进行的,相当于所得泥沙起动公式仅检验了ni>nb一种情形。关于冰盖下细颗粒泥沙起动的研究,由于目前缺少有关的试验数据,故对于公式的普适性,还有待进一步检验。

4 讨 论

(1)前人的研究大多集中在明渠泥沙起动问题上,对冰下泥沙起动则研究较少。本文公式基于Einstein假定,在窦国仁[11]提出的3种泥沙起动状态间关系、韩其为等[15]提出的相对暴露度概念等前人研究成果的基础上,考虑糙率变化对泥沙起动的影响,通过力学分析与概率统计理论相结合的方法,推导得到了冰盖流泥沙起动公式,并与试验数据对比验证,相对误差较小,能够满足计算要求,补充了已有研究成果对泥沙起动问题研究的不足。

(2)冰盖流与明流造床机理相同,床面上泥沙颗粒所受力的种类并无差别,冰盖的存在仅改变了床面附近以及整个垂线上的流速分布。当本文公式中冰盖糙率ni为0时,即可退化为明渠流泥沙起动流速公式。采用窦国仁整理的平坡水槽试验资料[25]及长江河道实测资料[26]对本文明渠流泥沙起动流速公式进行验证,并与沙玉清公式及武汉水利电力学院公式[27]的计算值一起绘入图5。

图5 明渠流起动流速公式计算值与实测值比较及相对误差分析Fig.5 Comparison between calculated values and observed data of incipient velocity of sediment under open channel flow and analysis of relative error

由图5看出,本文推导的明渠流泥沙起动流速公式计算值与窦国仁等[25]试验值整体上符合良好,相对误差范围为0%～0.86%,满足计算精度要求。但是,在颗粒粒径D≤0.03 mm的范围内,本文公式计算值与实测值相比偏小;分析原因可能是黏性细颗粒泥沙比表面积较大,界面化学效应突出,水体化学条件变化导致其结构变化[28],相比非黏性泥沙颗粒间受力分析更为复杂,本文公式推导引用的部分经验公式和经验取值可能对计算黏性细颗粒泥沙起动流速有一定影响,有待进一步检验。

(3)本文公式在张红武泥沙起动公式[3]基础上,引入糙率变化系数对附加下压力的表达式作了修正。当床面泥沙粒径D=0.32、0.85、1.32 mm时,本文公式计算值与实测值平均相对误差基本<10%,说明本文公式整体结构合理,基本能够描述冰盖下河床泥沙起动规律。不过在验证本文公式计算精度时,限于目前冰盖下泥沙起动的相关研究较少,仅对比了粗颗粒泥沙在ni>nb条件下起动流速与实测数据,缺少对于细颗粒泥沙在冰盖条件下的试验资料和变化规律的验证,故本文公式对于冰盖流细颗粒泥沙起动的适用性和可行性,还需进一步研究和验证。

5 结 论

(1)本文基于Einstein假定,采用滚动起动模型,在前人研究的基础上,考虑了冰盖流糙率变化对泥沙起动流速的影响,提出了冰盖流泥沙起动流速公式。并与王军等[7]试验数据验证,相对误差集中在1%～10%之间,满足计算精度要求,可用于天然河流条件下泥沙起动流速的计算。

(2)冰封初期,冰盖特征复杂,影响冰盖糙率变化的因素有很多,对于本文公式计算值与试验数据的相对误差,现有的研究较难给出理论分析,对其作用机理研究的还不够深入,还需进一步通过野外和室内系统的原型观测和模拟试验分析河床泥沙演变过程,精确率定泥沙颗粒起动公式中的系数,优化表达式,使其更符合实际。

(3)本文探讨了冰盖流条件下床面上泥沙起动的力学机制,与明渠流条件下河床泥沙受力类型相同,冰盖的存在只是改变了床面附近以及整个垂线上的流速分布。当本文公式中冰盖糙率为0时,即转化为明渠流泥沙起动流速公式,并与窦国仁等整理的试验数据对比,相对误差范围为0%～0.86%,满足计算精度要求,补充了已有研究成果对泥沙起动问题研究的不足。