苏 磊，蔺建铭，张宝森，谢志刚，田治宗，邓 宇，于国卿

(1.黄河水利委员会 黄河水利科学研究院，河南 郑州 450003;2.鄂尔多斯市水利事业发展中心，内蒙古 鄂尔多斯 017000）

1 引言

高寒地区冬季气温低,河渠易结冰而形成冰盖。冰盖显著改变河渠断面的水流结构,使封冻河道各水力要素变得十分复杂[1］,河道边界条件、水动力条件变化导致冰下水位、流量、流速变化,将影响冬季河渠输水,甚至诱发凌汛灾害,因此研究封冻河道水流特性具有重要的理论及实际意义[2］。

国内外学者对冰下流速分布问题进行了大量的研究,冰凌数学模型方面,国内外开发了一系列一维动/静态模型(RICE、ICEJAM、RIVER1D、HEC-RAS等)、二维动/静态模型(CRISSP2-D、DYNARICEDE[3］等)。王军等[4］、杨开林等[5］、茅泽育等[6］结合水力学、热力学原理改进了河道冰凌数学模型,但针对冰盖下水流垂向流速分布特性的研究相对较少。

本文对N-S方程进行简化、离散,采用k-ε紊流模型建立冰盖下水流垂向二维数学模型,模拟封冻河道冰盖下流速分布情况,并与实测资料进行比较,分析、探讨其变化规律。

2 控制方程与离散

2.1 控制方程

如图1所示,x表示流动方向坐标,z表示垂直于河床方向坐标。对于冰盖下的明渠紊流,考虑静压假定的不可压缩流体的控制方程无量纲化后,可以表示如下。

图1 坐标系示意

连续方程:

x方向动量方程:

在z方向上引入σ坐标,某一坐标处的垂向质点流速w可用垂向场流速ω表示为

式中:u、w分别为x、z方向的时均流速；t为时间；ζ为自由水面高程(ζ＝zb＋H),zb为床面高程；R0为罗斯比数；为科氏系数,g为重力加速度；ρ为密度；p为压强；AH为水平向紊动涡黏系数；AV为垂向紊动涡黏系数；H为水深。

此外,式中定义如下无量纲组合:①Rossby数,R0＝Ur/fXr；Froude数,Fr＝Ur/(gZr)1/2。②密度Froude数,。③水平向Ekman数,EH＝AHr/；垂向Ekman数,EV＝AVr/。其中下标r表示对应变量的参考值。

垂向紊动涡黏系数AV采用垂向k-ε模型进行计算,该模型可以较好地反映风剪切应力、床面剪切应力、流速梯度产生项、耗散、扩散和分层等的影响。k-ε模型的基本思想是将垂向紊动涡黏系数AV与紊流动能k和紊流耗散率ε联系起来,即

紊流动能k和紊流耗散率ε的输运方程表示如下:

式中:cν、c1、c2、σk、σε为常数,取值cν＝0.09、c1＝1.44、c2＝1.92、σk＝1.0、σε＝1.3；源项P为流速垂向梯度引起的紊流产生项；G为浮力产生或耗散项。

河床阻力的确定,以x方向为例,河床切应力在x方向表示为,其中:U、V分别为x、y方向的平均流速；cf为(无量纲)床面摩阻系数,cf＝gn2/H1/3＝g/C2,n为Manning糙率系数,C为Chézy系数；mb反映了河床底坡的影响。

冰盖阻力的确定,采用H.A.Einstein提出的阻力划分方法可将流体分为两层等效明流层(冰面层和床面层),假定上下两层互不影响,分别只受冰盖底部和床面粗糙度的影响,与表示河床阻力相同的方式,冰盖底部的切应力在x方向表示为。

2.2 离散求解策略

将x方向流速u表示为如下形式:

式中:u′为x方向流速相较于垂向平均流速的差值,即摄动流速。

u′所满足的方程可通过将三维控制方程减去二维控制方程得到,这样的处理移除了原始方程中的自由水面重力波项,因此求解u′后再通过式(7)间接计算u,能使模型获得更高的稳定性,保证计算质量。

对于摄动动量方程,除时间项外,只有垂直扩散项被视为隐式的,其余各项(包括科氏力项、对流项、水平扩散项和斜压项)均被视为显式的,也就是说,摄动动量方程可以简写为如下形式:

式中:Rx为x方向的水平源项,包括科氏力项、对流项、水平扩散项和斜压项4部分。

将式(8)采用有限体积法进行离散,将其在一个时间步长Δt和一个控制体积Δσ上进行积分可得

如图2所示,假设变量在中心点的值代表该变量在整个控制体积Δσ上的积分平均值,并认为隐式项在n、n＋1两个时刻的加权平均值代表其在整个时间步长Δt上的积分平均值,显式项在n时刻的值代表其在整个时间步长Δt上的积分平均值,则上式可表示为

图2 垂向网格示意

式中:θ为权重系数；k为垂向分层数。

代入公式(9),得到

上式即式(8)的最终离散形式,并对k＝2,3,…,m-1成立,其中m为垂向分层总数,并可写成通用形式:

其中

2.3 定解条件及方程组的求解

2.3.1 边界条件

(1)床面边界条件(k＝1),如图3所示。

图3 床面边界条件示意

式(8)在床面附近所需满足的边界条件为

代入式(9),得到

式(13)写成通用形式为

其中

(2)冰盖边界条件(k＝n),如图4所示。

图4 冰盖边界条件示意

式(8)在冰盖附近所需满足的边界条件为

代入式(9),得到

式(9)写成通用形式为

其中

2.3.2 初始条件

2.3.3 方程组求解

离散方程[式(12)］结合床面边界条件[式(14)］、自由水面边界条件[式(16)］,最终形成如下三对角方程组,该方程组可用“追赶法”求解。

求解上述三对角方程组获得摄动流速u′后,可由式(7)计算得到每层流速u。整个计算的迭代收敛条件为,其中ε取10-4。

3 计算结果分析

应用实测资料对上述模型计算结果进行了验证,如图5、图6所示,图中h为水深、y为相对水深、为断面平均流速。其中,图5为计算值与文献[7］中现场实测资料(River Moskva上长为400 m的封冻顺直河段第4号量测断面垂线流速分布实测成果)的比较情况,图6(a)为计算值与文献[8］中室内水槽实测资料(槽宽80 cm,水深8.92 cm)的比较情况,图6(b)为计算值与文献[8］中第12号水流条件(水深200 cm,冰盖糙率与床面糙率相同)的比较情况。结果显示,模型的计算值与实测值均吻合较好。同时,流速分布对于冰盖糙率ni和床面糙率nb的影响较敏感,最大流速Umax的位置偏向糙率较小的一侧,且在相同的水流条件下,ni/nb的值越小,最大流速点的位置越靠近糙率较小的一侧；当冰盖糙率和床面糙率相等时,流速分布以最大流速Umax为轴呈对称分布,且在最大流速Umax处流速梯度几乎为零[9］。

图5 计算值与文献[7］流速分布比较

图6 计算值与文献[8］流速分布比较

将冰盖下水流以最大流速为界分为床面层和冰面层两层,分层示意如图7所示,hb为床面层水深,yb为床面层相对水深,hi为冰面层水深,yi为冰面层相对水深。以文献[8］中两组工况为例子,工况1为床面糙率和冰盖糙率不相等时(ni/nb＝0.76),将其计算结果绘于图8所示的对数坐标系下；工况2为床面糙率和冰盖糙率相等时(ni/nb＝1),将其计算结果绘于图9所示的对数坐标系下。由图8可知,工况1当ln(yb/hb)＜-0.4时(即yb/hb＜0.67时),床面层流速分布符合对数分布,当ln(yb/hb)＞-0.4时(即yb/hb＞0.67时),床面层流速分布将不再符合对数分布,且流速值略小于符合对数分布的情况；在冰面层有同样的规律但影响区域略有不同,当ln(yi/hi)＜-0.32时(即yi/hi＜0.72时),冰面层流速分布符合对数分布；ln(yi/hi)＞-0.32时(即yi/hi＞0.72时),冰面层流速分布将不再符合对数分布,且流速值略小于符合对数分布的情况。同时,由图9可知,工况2当ln(yb/hb)＜-0.41时(即yb/hb＜0.66时),床面层和冰面层流速分布符合对数分布；ln(yb/hb)＞-0.41时(即yb/hb＞0.66时),床面层和冰面层流速分布将不再符合对数分布。

图8 对数坐标系下流速分布结果(文献[8］中,ni/nb＝0.76)

图9 对数坐标系下流速分布结果(文献[8］中,ni/nb＝1)

计算结果表明,封冻河道的断面流速分布主要与床面糙率、冰盖糙率及河道的水流条件等影响因素相关。床面糙率与冰盖糙率相等和床面糙率与冰盖糙率不相等两种工况,冰面层和床面层均有70%左右的流速符合对数分布,其余靠近最大流速点约30%范围内的流速不再符合对数分布且流速值略小于对数分布的情况。工况1当床面糙率与冰盖糙率相等时,冰面阻力和床面阻力干扰影响的区域相等,最大流速Umax在水深中点的位置。工况2当床面糙率与冰盖糙率不相等时,冰面阻力和床面阻力干扰影响的区域不相等,且糙率越大阻力越大,阻力干扰影响的区域越大,从而导致最大流速Umax位置偏向相对光滑的一侧；Umax的位置与ni/nb比值相关,且ni/nb的值越小Umax位置越靠近相对光滑的一侧。床面糙率与冰盖糙率相等和不相等两种情况,计算结果均表现出相同的规律,这与前人的试验成果吻合。

4 结论

本文基于浅水方程,采用k-ε紊流模型建立冰盖下水流二维数值模型,并将其应用于封冻河道冰下流速分布的计算。结果表明,冰盖下水流流速分布对于冰盖糙率ni和床面糙率nb的影响较敏感,当床面糙率与冰盖糙率不相等时,冰面阻力和床面阻力干扰影响的区域不相等,从而导致最大流速Umax位置偏向相对光滑的一侧,Umax位置的偏离程度与ni/nb值直接相关；当冰盖糙率和床面糙率相等时,流速分布以最大流速Umax为轴呈对称分布。床面糙率与冰盖糙率相等和不相等两种情况,计算结果均显示在最大流速点附近30%的区域内流速分布并不遵循对数分布规律,且流速值略小于对数分布的情况。