冰盖前缘冰块稳定性研究进展

王 军，桑连升，宋飞虎，程铁杰

(合肥工业大学，安徽 合肥 230009)

1 研究背景

寒冷地区河流冬季常出现冰盖、冰塞和冰坝这些特有的冰情现象[1]，当河流中形成完整冰盖时，水流从明流变成了封闭的暗流，水流湿周增加，上游水位升高，尤其冰塞或冰坝的出现，会导致上游水位短时间升高很多，可能因此诱发凌洪灾害，使河渠建筑物和堤防受到威胁[2-4]、河渠长距离输水流量受到影响等[5]；除此之外，冰盖下的河床演变也不同于明流时期[6-8]。封、开河时的冰塞或冰坝是最容易致灾的[9-10]，上游来冰在到达冰盖前缘时，或停留使冰塞前缘向上游延伸或受水力作用而下潜于冰盖下。由于上游来冰在冰盖前缘下潜与否是冰塞或冰坝演变的重要组成部分和过程，也是冰塞演变数值模型的重要基础，所以自1920年代就有学者关注研究冰盖前缘冰块的稳定性问题[11-12]，但因问题的复杂性，至今其实没有真正的理论解。目前，河冰数学模型虽然已经发展到了可采用离散元模型模拟冰盖的形成与变化[13-16]，但对于冰盖前缘冰块下潜与否仍然是依靠经验公式或水力条件的经验值进行判断。

2 冰盖前缘冰块稳定性研究

一直以来，冰块在冰盖前缘稳定性问题因其重要性受到了众多学者的关注，也已经取得了很多重要成果，这些研究成果按其研究方法大致可分为三类，分别为原型观测研究、理论分析、试验研究等，对该问题的研究属于实际河流应用的需要，所以有些相关研究成果实际上不能单纯的划归其中某类，某种程度上这样的分类只是为了叙述问题的方便。

2.1 原型观测研究方面文献[11]中，有关St. Lawrence Waterway Project的研究发现，当河流中的水流速度V<0.4 m/s时，可能形成光滑冰盖；V<0.7 m/s时，冰块和冰花不会下潜；而V>1 m/s时，不利于河流冰盖的形成，这或许是目前文献中最早有关冰块下潜临界值的描述。由于单一采用流速作为临界标准存在缺陷，水流弗劳德数被引入进来。文献[12]中，孙肇初根据黄河河曲段原型观测资料，发现非下潜的临界流速Vc=0.6 m/s，下潜临界水流弗劳德数Fr=0.09；水利部北京院和西北院根据1962到1965黄河刘家峡至盐锅峡河段的大型冰塞现象观测分析，得出该河段产生冰塞的下潜临界水流弗劳德数也是0.09；Cartier[17]根据小型河渠现场研究得出的下潜临界弗劳德数是0.13；Kivisild[18]的现场观测得出临界下潜弗劳德数是0.06 ～ 0.1之间变化，Oudshoorhm[19]根据在Rhine river下游部分的原型观测，得到的下潜临界弗劳德数是0.06 ～ 0.09；Shen等[20]根据原型观测，建议在一般情况下下潜临界弗劳德数为0.05 ～ 0.09；Beltaos[21]分析了众多原型中封河冰塞的研究文献后，得出了和冰块下潜有关的数据，即温和的冬季，冰塞发展大致稳定后，冰塞下的水流速度一般在0.9 m/s左右，高可以至1.2 m/s，也可以在卡封处降至0.6 m/s；在冰盖形成且冰塞达到一定厚度时，冰塞下流速达到0.8 m/s时冰盖一般不会破裂。上述下潜临界水流弗劳德数的变化范围大致为0.06 ～ 0.13，从现象上看，河流越大，临界值相对越小，但从问题的实质考虑，到达冰盖前缘的冰块下潜与否，不单纯是下潜流速或水流弗劳德数单一能够描述的，冰块下潜与否受众多因素制约，导致这么大的变化范围原因大致分析如下，首先是人为因素，各自判断标准有差异；其次，各个河流上游的来冰量不会相同、弯道所导致的水流变化各异等[22]，使冰块在冰盖前缘的下潜形式不尽相同；再者冰盖前缘的冰块尺寸大小及其形状、冰块自身的性质[23]、冰块与水流密度的差异(可能含沙等)和冰期各河流的水流条件等不会一样；最后还有风速和风向的影响等，如此，造成了不同的河流、不同的水流和冰的条件下，临界水流弗劳德数不尽相同，有一变化范围。

2.2 理论分析方面如前所述，冰块下潜研究很多是采用了多个方法相结合的办法进行研究的，所以，这里介绍的是以理论分析方法为主的这类研究。由于众多学者研究时，采用的相关变量符号有所差异，为叙述问题方便，本文试图统一符号但保持原来的物理意义。Pariset和Hausser[24-25]基于理论分析和试验研究发现，当水流速度较低，冰块遇到障碍停止后，后续来冰逐渐聚积并向上游发展，此时形成单层厚度的冰块层，冰块稳定条件可描述为：

(1)

式中：Vu为冰(块)盖下的水流平均速度；g为重力加速度；ρ和ρi分别为水和冰块的密度；t为冰块厚度；k为反映冰块尺寸和形状的因子。

Pariset和Hausser[24-25]等通过室内木块和塑料块试验以及户外的真冰试验，得到了k值变化范围在0.66 ～ 1.3之间，其中0.66差不多对应立方体的冰块形状，1.3对应于厚度相对薄且扁平的形状。天然河流中，这种以单个冰块厚度形成冰盖情形仅出现在水流速度很小的情况，一般在流凌达到一定密度、流冰向下游运动受阻停留，停留处会有冰块下潜和堆积并形成一定相对稳定的厚度，上游来冰不断聚积、下潜和堆积并向上游发展。针对这种情况，Priset和Hausser根据基本的一维水流伯努利方程和不漫溢条件，得到了维持冰盖前缘冰块稳定的最小冰塞厚度和上游行近水流条件的关系：

(2)

式中：V为上游行近流速；H为上游行近水深。

式(2)可以变化为：

(3)

Michel[26]认为上游水流绕过冰块时，由于局部的水流加速和分离，引起冰块底面压力减少，浮力作用减弱，经分析得出的临界下潜条件为：

(4)

式(4)可重组为：

(5)

比较式(3)(5)，可以看出，对于同样的t/H，式(5)算出的临界下潜速度相对要大。对于式(4)，为考虑冰塞体内冰堆积的空隙影响，Michel在计算式中引入了孔隙率Pj：

(6)

并且Michel得出了类似Pariset和Hausser的式(3)：

(7)

式(7)比式(4)多了孔隙率Pj的影响，只是Michel在推导时，孔隙率是指上游冰盘的孔隙率，而推导过程中对于冰盖下的冰塞体采用了相同的孔隙率，显然两者是有区别的。

根据式(7)，由图1可见，在Fr=0.06 ～ 0.13范围内，可以看出即使对于同样的Fr，孔隙率Pj在一个很大范围内变化，对于较大的孔隙率，可以理解为对应于初始冰盖形成阶段，而较小的孔隙率则可以理解为对应于稳封阶段后期至开河前夕阶段。

图1 Pj与t/H关系图

Uzuner和Kennedy[27]根据伯努利方程和水流连续性方程，推导得出了冰块下的压力分布p(x)的计算式，这是最早对冰块下的压力分布进行的理论研究。

(8)

图2 冰盖前缘冰块稳定性分析示意图

由p(x)对点D产生的力矩为：

(9)

式中Mp为冰(块)盖下的压力产生的力矩。

当冰块向下旋转时，AB面的压力增加，产生一个顺时针力矩，靠近B点的底部分流区BC因压力明显降低产生一个逆时针力矩，因问题的复杂性，Uzuner和Kennedy考虑了将两者的联合作用一并考虑，将力矩Mh表示为：

(10)

式中Cm为力矩系数，与t/H，t/L，ρi/ρ以及转角θ有关。Cm可正可负，取决于上述两个力矩的相对大小。

重力产生的力矩如下式所示：

(11)

初始旋转角度θc是对应不漫溢条件(no spill condition)相对应的旋转角度，由下式计算：

(12)

按照力矩平衡：Mp-Mh-Mw=0，可以计算出下潜临界弗劳德数Fr。

Uzuner和Kennedy[27]对应于上面的理论分析，对比进行了试验研究，其试验条件设置大致为：试验水槽的长宽高分别为3.66 m，0.3048 m和0.6096 m；模拟冰块形状为矩形六面体，ρi/ρ涉及三种取值，分别为0.5，0.67和0.87；模拟冰块有两种厚度，分别为0.03175 m和0.0714 m；t/L变化范围为0.096 ～ 0.773；一些试验中的模拟冰块在上游下表面修成四分之一圆弧，其t/L范围为0.25 ～ 0.26，这些试验点据的临界下潜弗劳德数相对要大一些。Cm与t/L关系曲线如图3所示，Cm在t/L大约等于三分之一处获得最大值，说明了当冰块厚度与长度之比在1/3左右时，冰块较容易下潜。

图3 Cm与ρi/ρ、t/L关系图

而对于较长(t/L<0.1)和较短(t/L>0.8)的模拟冰块，Uzuner和Kennedy经过分析得出了：

(13)

式中β是几何尺寸和流体密度的函数。

一般情况CsF22比1小得多，若在上述括号里面分母部分忽略这一项，可以进一步得到：

(14)

式中β=0则可应用于较长的冰块计算。

Uzuner和Kennedy[27]从一般力学原理出发，对冰块下潜进行了研究，结合试验分析，依据t/L的比值对冰块下潜进行了分类，即下沉下潜和翻转下潜，并分别给出了计算方法，其影响是深远的，迄今为止，冰块下潜分类和计算思路仍然沿用的是Uzuner和Kennedy提出的方法。

Ashton[28]认为t/H比t/L对下潜的影响更大，但t/L也并非无关紧要，因为当t/L增加时，下潜临界值会相应减小。并通过简单的力矩平衡推导，得到了：

(15)

(16)

并通过和Uzuner和Kennedy[27]研究成果比较，认为其研究结果更合理、更符合试验结果。

尴尬的是，无论是Pariset和Hausser[24-25]、Uzuner和Kennedy[27]还是Ashton[28]等，其伯努利方程的推导均是按照下潜冰块宽度与河道水面宽度相同进行的，而实际上，天然河流的上游来冰一般很难和河道水面等宽，所以这些对后人研究影响极大的研究结果其实是建立在了一个似是而非的基础上。

Wong和Beltaos[29]对式(16)进行了重组，目的是能够直接使用流量对冰块下潜进行判断。

Daly和Axelson[30]从考虑下潜力矩MU和抗倾覆力矩MR以及冰块旋转角θ方面出发，认为MR与冰块形状、几何尺寸和冰块及流体的密度有关，而MU与冰块形状、几何尺寸和水流速度与深度以及流体密度有关，并且，MR会在某个θmax处达到最大值MRmax，MU则随θ增加而增加。因此，若θ*≤θmax，有MU(θ*)≤MR(θ*)，此时冰块不会发生下潜，若MU达到并超过MR(θmax)，冰块会发生下潜。但下潜力矩MU和冰块旋转角θ关系复杂，涉及块体分离区、加速区、剪切应力等。因此Daly和Axelson采用了类似Uzuner和Kennedy[27]和Ashton[28]的假定，即下潜力矩MU与作用在块体上的动水压力ρV2和块体上的作用面积成比例。在计算抗倾覆力矩时，Daly和Axelson采用了以下算式：

MR(θ，ρi，ρ，t，L，B)=[ρigVd×rd(θ)]-[(ρgVp(θ))×rp(θ)]

(17)

式中：t、L、B分别为冰块厚度、长度和宽度；Vd为冰块体积；Vp(θ)为冰块排开的水体积；rd(θ)为旋转点到块体重心的距离；rp(θ)为从旋转点到浮力重心的距离。

Daly和Axelson具体给出了以下计算式(参见图2)：

0≤θ≤θ1：

(18)

0≤θ≤θ2：

(19)

θ2≤θ≤π/2：

(20)

式中：θ1为浮动冰块上游角A点淹没的旋转角；θ2为浮动冰块下游角E点淹没的旋转角；t′=ρit/ρ。

(21)

(22)

抗倾覆力矩在θ=0时为零，随着旋转角度的增加而迅速增加，然后缓慢增加到θmax时的最大值，然后随着旋转角度增加而缓慢减少。MRmax和θmax可以通过数值求解式(18)—(22)得出，具体计算结果见文献[30]。

Uzuner和Kennedy[27]较早提出了按照力矩平衡分析冰块在冰盖前缘稳定性的方法，但对于下潜力矩和抗倾覆力矩的计算，依靠力矩系数Cm确定。Daly和Axelson[30]则更全面地提出了以冰块转动角θ为变量的下潜力矩和抗倾覆力矩计算方法，看起来似乎比Uzuner和Kennedy前进了一步，但无论下潜力矩还是抗倾覆力矩，随着转动角的不同，冰块下的分离区以及压力大小和分布都在变化，完全按照力和几何关系写出的表达式是不能反映这种变化的，所以从问题的实质性来看，Daly和Axelson只是在理论计算方法上前进了一步。

McGilvary和Coutermarsh[31-32]结合考虑冰块角动量对冰块在冰盖前缘的稳定性进行了研究，通过对抗倾覆力矩(MR)和下潜转动力矩(MU)分析，给出了冰块转动运动方程可表示为：

(23)

McGilvary和Coutermarsh[31-32]从能量角度考虑冰块的旋转，即流体对冰块累计做功等于冰块的转动能量，得出下式：

(24)

Ambtman等[33-34]通过试验，对冰块下压力分布进行了测量，发现压力分布沿冰块的中心线对称，沿宽度方向基本不变，冰块底面压力分布的减少可分成两个部分，即文丘里效应导致的压力降低和冰块前缘效应导致的压力降低；随着冰块厚度与行近水流深度之比和水流速度的增加，沿冰块中心线的动态压力降低，下潜力和翻转下潜力矩增大；通过测量具有圆形前缘模拟冰块中心线的压力分布，发现冰块前缘形状对流体分离的影响及其压力分布非常重要。在此基础上。Ambtman依据其试验测量结果，结合Daly和Axelson力矩分析方法和计算公式，给出了下潜力和下潜力矩的计算办法，并与 Uzuner和Kennedy[27]、Ashton[28]、Daly和Axelson[30]和Larsen[35]等的试验结果和计算进行了比较，认为其研究成果和这些学者的研究吻合较好。但是，Ambtman等在冰块底面压力分布的试验研究和理论分析中，其模拟冰块和水槽等宽，长度为0.5 m，依靠4个螺纹杆固定，模拟冰块前后都是明流情况，且研究中忽略了模拟冰块的下游效应和边缘效应，天然情况下，上游冰块很少与河道等宽，大部分情况下，冰块宽度和河道宽度相比可能小得多，下潜冰块到达冰盖前缘时其下游不会是完全的明流情况，并且其研究结果不适用于冰塞前缘厚度大于冰块厚度的情况，天然情况下，可能很多时候恰恰相反，冰塞前缘厚度比上游来冰块厚度要大。因此，该研究成果其实有很多尚待提升改进的地方。

2.3 试验研究方面如前所述，冰块稳定性研究很少是采用单一研究方法的，以下的叙述主要是以试验研究为主(包括数值试验)的研究成果。

Beltaos[21]在总结其前面相关学者研究基础上，依据量纲分析，认为冰块下潜临界弗劳德数可表示为：

(25)

式中：ti、li、bi分别是冰块厚度、长度和宽度，其下标i可以理解为冰块尺寸大小不一；si=ρi/ρ。

Beltaos认为Pariset和Hausser[24]、Uzuner和Kennedy[27]和Chee和Haggag[36]的研究结果都是上式不同形式的具体表达。但也正如Beltaos所说那样，有关式(25)分析解的表达式至今未见，式(25)中，明确出现了li/bi因子，迄今为止，众多研究者的计算式中几乎没有出现该因子。

图4 黄河河曲段冰盖前缘处冰块下潜临界条件参数关系图与相关图

(26)

隋觉义等通过研究认为，Fc与t/H关系密切，t/H大，Fc则小，反之则大。

隋觉义等不仅解决了当时河曲段冰情计算的需要，而且是世界上为数不多的原型河道上的真冰试验。

王军[38]根据试验研究认为，下潜临界冰厚弗劳德数可表示为：

(27)

因为试验条件中的si=ρi/ρ为不变量，所以式(27)中没有出现该变量。根据试验数据，回归得出：

(28)

式(28)明确考虑了下潜过程中冰块宽度因子的影响。

练继建等[39-40]基于数值模拟方法，对冰块底面的压力分布进行了精确数值模拟，其研究结论如下：(1)冰厚、流速增大时，前缘效应和文丘里效应增强，冰块前缘的最大压差影响范围增大，压强恢复速度减慢；(2)前缘为圆形断面的冰块较前缘为矩形断面的冰块有着较弱的前缘效应和相同的文丘里效应。并在Ashton[28]公式基础上，提出了自然破碎状态下冰凌下潜的判断准则为：

(29)

式中：k为修正系数，取值范围为1.15 ～ 1.35之间。该式可理解为考虑冰块前缘形状条件下对Ashton公式的一种修正。

Wang等[41]采用四种不同尺寸的冰块，通过考虑冰块运动以及由于水流流动所产生的阻力、碰撞力和压差力，考虑冰块的长宽高三个方向尺寸，对冰块在冰盖前缘的稳定性进行了试验研究，给出了冰块在冰盖前缘是否下潜的临界判别式，该式计算值与试验结果吻合较好。Wang等[42]基于试验研究，通过对冰块的力矩平衡分析，进一步给出了少量冰块下潜和大量冰块下潜的临界判别式。文献[41-42]的研究中包含了对冰块宽度的考虑。

3 结论与建议

冰块下潜临界条件研究是冰塞和冰坝研究的基础，预测冰塞或冰坝的相关数学模型也离不开冰块下潜临界条件的判别。通过对历年来典型研究文献的梳理，相关总结如下：

(1)下潜临界流速和弗劳德数的判别方法至今在工程中仍然广泛使用，但冰块在冰盖前缘的下潜行为受众多因素制约，如水流条件、冰量条件、冰盖厚度变化、气候条件等，因此，不同的河流上乃至同一条河流不同的条件下，临界弗劳德数存在变化范围，区域和局部范围的经验临界值仍然是需要的，如何更加科学合理的给出下潜临界条件的理论表达式，仍然需要今后的努力。

(3)从理论探求的需要，冰块在下潜时其底面的压力和分布变化是非常重要的，直接影响下潜和抗倾覆力矩的正确计算和表达。尽管有些学者对此进行了探索，但如Ambtman[33-34]等研究的那样，无论试验装置还是模拟冰块的状况，和实际冰块下潜发生的情况都相差较大。今后的研究中，应结合天然实际情况，从试验研究和数值模拟研究两个途径对该问题展开研究，进一步从机理层面探寻冰块下潜变化的内在因素和规律。

综上所述，冰块下潜研究是一个古老却又不失现代的课题，还有很多方面需要改进探索，包括冰块的形状本身，就是一个很复杂的问题，现今的研究基本视其为矩形六面体，而天然河道盘型冰块比比皆是，形状各异；另外，冰塞演变过程中，冰盖或冰塞厚度是变化的，而且大多比冰块厚度大，因此，更加努力的贴近实际和自然，是发展至今天的冰块下潜研究应该追求的方向。