滦河下游腰庄段丁坝群水流特性及局部冲刷三维数值模拟

李大鸣,肖程之,李佩瑶,李彦卿,卜世龙

(1.天津大学 水利工程仿真与安全国家重点实验室,天津 300350; 2.生态环境部海河流域北海海域生态环境监督管理局生态环境监测与科学研究中心,天津 300170)

0 引 言

滦河下游丁坝群是滦河流域防洪体系中主要组成部分,为减轻滦河下游沿岸的洪水威胁以及为河道两岸人民群众的生产生活提供安全保障。20世纪50年代末开始在滦河修建丁坝工程,由于当时丁坝整体性较差、间距大、防冲效果差,造成的丁坝损坏情况也比较严重。由于流场的改变,丁坝附近泥沙冲淤状况将重新调整,从而产生局部冲刷和淤积现象。丁坝的局部冲刷直接影响水工建筑物的安全,因此分析局部冲刷一直是丁坝研究的重点之一。

为提高丁坝的稳定性与安全性,许多学者从丁坝间距、丁坝长度、丁坝挑角及丁坝材料等方面研究丁坝的局部冲刷。由于物理模型试验形象、直观等优点得到了广泛的应用,学者们开始通过设计模型试验模拟丁坝局部冲刷的实际情况。Gu等[1]通过比较不同渗透性堤防的水流与输沙特性,发现透水堤防悬沙通过堤防的纵向运输远大于横向运输。Li等[2]探讨了清水条件下直壁丁坝局部冲刷坑几何结构的空间特征。Pandey等[3]确定了砂砾河床条件下主要影响冲刷坑深度的参数,并以此推导出新的丁坝冲刷预测模型。Özyaman等[4]对泥沙级配、水流深度、丁坝角和丁坝长度对冲刷过程的影响进行了研究,并提出4个方程以预测冲刷深度与冲刷量。Hajibehzad等[5]设置了弯曲水槽中4组不同间距丁坝的冲刷试验,分析了冲刷坑深度以及冲刷位置与丁坝间距之间的关系。喻涛等[6]通过水槽概化物理模型试验,研究了丁坝在非恒定流作用下的水力特性和冲刷机理。张林忠等[7]借助模型试验分析了不同工程条件下不同间距丁坝群局部冲刷情况。

计算机运算能力的飞速提高为复杂水沙模型的构建提供了基础。Liao等[8]为解决普通二维模型不足以模拟三维流场和局部冲刷的缺点,通过将休止角公式和河床几何调整机制整合到二维动床模型中,以改进局部冲刷坑的数值模拟。Haider等[9]将丁坝布置形式分为“L”型与“T”型两种形式,采用k-ε模型(k为湍动能,ε为湍动能耗散值)研究了两种形式在矩形明渠水流相应的流动和湍流特性。胡志毅[10]通过对不同挑脚的丁坝附近流场及冲刷情况进行了三维数值模拟;发现挑角与最大冲刷深度、坝前及坝后冲坑长度之间的关系。为研究弯道处不同长度丁坝对水流特征的影响,张岩等[11]对60°弯道内的3种长度的丁坝绕流进行了三维数值模拟,得出坝后湍动能及湍流黏度作用强度及范围随丁坝长度增加而增大。戚蓝等[12]通过三维数值模拟方法对马良子段丁坝群冲刷坑发展进行模拟,提出了丁坝长度调整方案。吕庆标等[13]对无丁坝保护的河道岸线崩溃机理进行了研究,说明了无保护岸线在特定水流条件下可能发生破坏的危险性。张小雅等[14]研究了液压坝群对河道冲淤变化特性的影响,为不同水流条件下河床稳定和堤坝安全提供了可供参考的条件。

近年来对于丁坝局部冲刷的研究取得了很大进展,但以往研究大多未考虑试验与天然河道的差别,如河道断面不均匀、流量过程为非恒定流等。同时大多数学者通常仅研究单丁坝影响水流和冲刷特性的因素,但对丁坝群的研究较少。而在天然河道中丁坝工程往往采用多个丁坝相互配合以提高整体效果。与单个丁坝相比,丁坝群内各丁坝流场相互影响,冲刷问题复杂多变。为研究丁坝群附近流场及局部冲刷情况,本文根据实际丁坝群河道整治工程构建三维水流泥沙数学模型,通过模拟丁坝群附近水力特性,探讨丁坝局部冲刷情况。

1 三维水流泥沙模型数值模拟理论基础

1.1 基本控制方程

基本控制方程如下:

(1)

(3)

(4)

1.2 RNG k-ε模型

k方程为

ε方程为

(6)

式中:ε为湍动能耗散值;μt为涡黏性系数,μt=Cμρk2/ε,Cμ是模型常数,Cμ=0.085;τij为雷诺应力;C1ε、C2ε、σk均为模型常数,C1ε=1.42,C2ε=1.68,σk=0.717 9。

1.3 泥沙冲刷模型

模型的计算过程通过考虑泥沙的悬浮和沉积2种存在状态进行。悬浮泥沙通常浓度较低,随流体流动而发生平流。堆积沉积物的淤积量可以通过临界填充率(默认值为0.64)的调节进行处理。只有堆积的沉积物表面薄层颗粒(厚度仅为几倍粒径)才能以推移质输运的形式移动。

在剪切力和小涡流作用下,堆积界面上的泥沙被挟起和再悬浮。由于不可能计算每个沉积物颗粒的流体动力,因此必须使用经验模型。本文使用的经验模型基于Mastbergen和Van den Berg公式[15],通过Soulsby-Whitehouse方程[16]预测临界希尔兹数。

计算临界希尔兹数的第一步是求解无量纲参数d*,i,即

(7)

式中:ρi为沉积物i的密度;ρf为流体的密度;d50为泥沙中值粒径;μf为流体的动力黏滞系数;‖g‖为重力加速度。

再由Soulsby-Whitehouse 公式[17]可计算临界希尔兹数θcr,i为

(8)

临界希尔兹数可以通过修改倾斜面的休止角进行调整,因为在倾斜界面处,填充的沉积物不太稳定,更容易被沿斜坡向下移动的流体夹带。

调整后的临界希尔兹数θ′cr,i表达式为

(9)

式中:β为河床的坡度角;Ψ是流体和上坡方向之间的夹角,若水流直接向坡上流动则Ψ=0°;φi是用户定义的不同沉积物种类的休止角(默认值为32°)。采用θcr,i=0.05。

然后计算泥沙的起悬速度ulift,i为

(10)

式中:αi为起悬参数,建议值为0.018[18];di为泥沙分组级别粒径;d*为di对应的无量纲参数(由公式(7)d50换为di确定);θi为d*对应的希尔兹数(由公式(8)d*,i换为d*确定);ns为沉积界面外法线方向余弦。

沉降速度usettling,i计算采用Soulsby提出的方程[18],即

式中vf为流体的运动黏滞系数 。

推移质输运是由于泥沙在沉积床表面滚动或跳跃而产生的泥沙输运方式。对于推移质输运,目前使用的模型由Meyer、Peter和Müller等人提出[19],该模型预测了沉积床表面单位宽度的泥沙体积流量,即

Φi=βi(θi-θ′cr,i)1.5cb,i。

(12)

接着需要对床载厚度δi(即跃变沉积物的厚度)进行估算。用来估算厚度的关系式为

(13)

为了计算泥沙在每个计算单元中的运动情况,将qb,i的值转换为床载速度ubedload,i,即

(14)

式中fb为沉积物的临界堆积分数。

2 模型的建立与验证

2.1 研究区域概况

滦河流域位于39°10′N—42°35′N, 115°40′E—119°20′E,控制流域面积4.47×104km2。滦河干流河道长度约900 km,河道平均坡降为2.68‰。滦河发源于河北省丰宁满族自治县大滩镇孤石山村东南部小梁山南麓,多伦以上属高原,地势平坦,河道坡降约为0.5‰;郭家屯以下至潘家口河段穿行于燕山峡谷间河谷宽为200～300 m,河道坡降约为3.33‰～1.67‰,河道蜿蜒曲折;潘家口水库以下河宽200～500 m,过桑园峡谷进入迁安盆地,河谷宽至1 000 m以上。京山铁路桥以下进入平原区,经主河道干流宣泄入海,平原区河道坡降为0.66‰。流域内建有潘家口、大黑汀、滦河支流伊逊河上庙宫水库、青龙河上桃林口水库等大型水库及其它防洪措施。

2.2 丁坝群概况

丁坝群位于滦河下游腰庄险工处,对应堤防桩号LHY68+200—LHY70+100。该工程布置丁坝10道(如图1所示),各坝头连线与河道主流向平行,坝身长度为50 m,均为下挑式,倾角60°～85°,丁坝间距为坝长的2～3倍,其中丁坝B1—B6改变主流流向,间距取100 m,为2倍坝长;丁坝B7—B10对主流方向改变不大,坝间距取150 m,为3倍坝长。坝顶宽度均为2 m。

丁坝坝头采用格宾笼抛石结构,坝头顶部宽5.0 m,长6.0 m,临水侧为圆弧形按坝头高度不同砌筑2～4层台阶,每阶高0.5～1.0 m,宽2.0 m,坡比1∶2。坝头基础分为干场和有水两种情况,干场采用明挖,开挖后用50 cm厚格宾石笼进行防护;有水时先用抛石回填至水面以上0.3～0.5 m,再铺设50 cm 厚格宾石笼护垫,如图2所示。

图2 腰庄丁坝坝头结构Fig.2 Structure of groin head

2.3 模型建立

模型计算范围由滦河马良子段丁坝下游至入海口,地形数据来源于2015年河北省水利水电勘测设计研究院的实测值。由CAD软件提取高程点后,在Mike21软件中划分二维网格后沿垂直方向从水面到河底进行等距插值得到三维网格,最终在Flow3D中完成编译计算。

计算模型长2 518 m、宽1 518 m、高20 m,选定河床主槽糙率为0.025,泥沙颗粒半径为0.062 mm,密度为2 650 kg/m3,临界希尔兹数为0.05[12]。丁坝群之间的水流流态较为复杂,网格大小会对计算结果产生直接的影响。因此在划分网格时通过矩形网格划分,将整个计算区域划分为A、B两区以减少不必要网格的数量,A区定床,B区为动床。采用嵌套网格的方式对坝体附近进行局部加密以保持更加完整的丁坝外观结构及附近水流形态。如图3所示,总网格数量约160万。

图3 计算区域网格划分Fig.3 Meshing of computational domain

2.4 模型验证

2015年为丁坝群初试阶段,工程初试阶段布置了3道丁坝。初试工程丁坝建成后,经历了2015年滦河洪水过程的冲刷,因此本文以初试工程丁坝为研究对象对三维泥沙模型进行验证。通过对2014年洪水过后初试工程附近实测地形提取生成地形数据,采用相同参数构建模型,如图4所示。

计算时采用有限体积法离散控制方程以及重整化群(Renormalization Group,RNG)k-ε湍流模型。上游设为流量入口,根据2015年行洪情况给定模型3 d洪水过程如图5所示,平均水深为5 m,下游出口设为自由出流,上表面采用流体体积(Volume of Fluid,VOF)法捕捉自由液面,侧面及下壁面设为无滑移边界[12,21]。

图5 2015年洪水过程Fig.5 Flood process in 2015

根据2015年洪水过后的实测地形图发现,丁坝群坝头处冲刷严重,丁坝B1坝头处有较深的冲刷坑,冲坑最深处达到了-7 m。因此,选取丁坝B1附近最深冲刷坑,将过流断面方向视为横断面,与之垂直的断面视为纵断面,如图6所示,将数值模拟的冲刷坑断面地形与实测地形值进行对比。

图6 丁坝B1附近冲刷坑断面选取Fig.6 Section of scour pit near spur dike B1

由图7(a)可以看出丁坝B1横断面冲刷坑模拟值与实测冲刷地形之间拟合较好,图7(b)冲刷坑纵断面30～50 m处误差较大。

图7 丁坝B1坝头附近冲刷坑横纵断面形态对比Fig.7 Comparison between cross and longitudinalsectional shapes of scour pit near the head of spur dike B1

产生误差较大的原因有:

(1)数值模拟中网格划分较大,计算结果不够精确。

(2)模型所计算的点为计算网格的平均值,并不能与实际值一一对应,导致两者整体上存在一定误差。

(3)模型中底沙分布较为理想,与实际情况存在一定的差异会对水流结构与泥沙运动造成一定影响。

虽然模型在定量验证中存在一定误差,但是模型基本能反应实际工程中水沙运动及冲刷情况,具有较好的可靠性,该模型可用于对复杂的天然河道下丁坝群附近局部冲刷的模拟。

3 丁坝群水流流速特性

3.1 模型计算方案

模型针对研究区域内5 a一遇、10 a一遇来流过程下恒定流和非恒定流情况进行模拟计算,计算工况见表1,不同重现期非恒定流过程见图8。为保证丁坝水力特性和局部冲刷达到稳定状态,对每种工况模拟时间为24 h。

表1 计算工况Table 1 Computation conditions

3.2 丁坝群流场分析

3.2.1 丁坝群垂向流场分析

为了分析丁坝群垂向水流流态,以工况1为例选取靠近河床附近层面、自由液面附近层面及中间层面进行分析,如图9所示。

不同层面坝头处均形成了缓流区(在挑流坝作用下,河道主流出现急流区,在坝头迎水面和背水面形成缓流区。水力学定义急流为干扰波不能向上游传播的流速,缓流为干扰波能向上游传播的流速。),对比3个截面可以发现随着高程的增加,缓流区面积也逐渐增大。同时各层面均能观察到坝田区内的漩涡。总体来说不同层面处水流流态分布大致相似,因此对其他工况仅选取中间层面进行分析。

3.2.2 丁坝群平面流场分析

图10为不同工况下最大流量时丁坝群中间层面流场图。从图10可以看出左岸(图中偏上位置)流速大于右岸流速,说明丁群坝工程对河岸起到保护作用。水流流经首坝B1时,由于丁坝的挡水作用在坝前形成一个水流滞留区,同时产生了一个较小的顺时针漩涡。丁坝B1坝体起到了明显的挑流作用,在丁坝B2—B5处产生了一块较大的缓流区。丁坝B1—B5的坝田区产生了较大的顺时针方向漩涡,漩涡形状接近方形,漩涡中心位于相邻两坝中间。丁坝B6较丁坝B2—B5位置略往主流区突出,根据分析比较各丁坝坝头处流线疏密程度,可以得到丁坝B6坝头处流线更为密集的结果,说明丁坝B6起到的挑流效果较明显。而丁坝B6—B9之间的间距为三倍坝长,丁坝群的整体性较丁坝B2—B5有所减弱,坝头处的缓流区面积较小,而丁坝B2—B5坝田区漩涡较B6—B9坝田区漩涡更为规则。由于B10位置向主流区突出,也直接受到主流冲击;同时B10下游没有丁坝,水流能够充分扩散,坝后无漩涡产生。

图10 丁坝群中间层面流场Fig.10 Flow fields at middle layer under different conditions

3.3 丁坝群水流流速变化

3.3.1 丁坝群流速分布

以工况1为例分析坝头处流速沿程变化情况,如图11所示。由图10流场图及图11可以发现,丁坝B1坝头处流速最大,为5.0 m/s。丁坝B2—B5处于丁坝B1掩护范围内,坝头处流速均小于B1。随与B1距离增加,B2—B4坝头流速在3.65 m/s左右,至B5增大至3.99 m/s,说明水流到B5处开始恢复。随水流逐渐恢复,丁坝B6坝头受主流冲击,坝头处流速明显增大,流速为4.55 m/s。丁坝B7—B9在缓流区,坝头流速均小于B6,由于丁坝B7—B9间距较大,坝头处流速大于丁坝B2—B5,坝头流速在4.05 m/s左右。这一规律与前人分析研究一致[22]。

图11 工况1各坝头处流速变化Fig.11 Variation of flow velocity at groin heads in condition 1

3.3.2 横断面流速分布

河道中布置丁坝后,水流受到丁坝的影响将重新调整流态,横断面的流速分布随之发生变化。因此研究丁坝横断面的流速变化,对于分析丁坝群的水流特性及周围河床冲刷具有重要意义。由于丁坝群共布置10道丁坝,丁坝数量较多,因此仅选择其中5道进行分析。丁坝B1、B6、B10坝头处流线集中流速较大,因此选取这三道典型丁坝。丁坝B1—B6改变主流方向,坝间距为两倍坝长,以B3为代表;丁坝B6—B10间距为三倍坝长,以B8为代表。恒定流情况下各丁坝横断面流速分布如图12所示。同一断面不同工况流速对比如图13所示。

图12 恒定流情况下各丁坝横断面流速分布Fig.12 Flow velocity distribution in cross section of spur dikes under constant flow

图13 不同工况下各丁坝横断面流速分布Fig.13 Flow velocity distribution in cross section of spur dikes under different conditions

通过对比图12及图13,两种工况横断面流速整体上呈现中间大两侧小分布,且坝头处流速小于主流流速。工况1与工况2主流流速分别约为6.4 m/s和8 m/s,工况1主流最大流速区在距右岸125—175 m处,工况2主流最大流速区在距右岸100—150 m处。由于河道地形左高右低,当洪水来临时,河道右侧积蓄大量洪水,流速较大。两种工况横断面流速梯度均为右岸大于左岸,这是由于在右岸布设丁坝群导致水流流速小于左岸。丁坝B1、B6、B10在距离右岸50—100 m处水流流速均呈“V”形分布,由于漩涡的产生水流紊动较为剧烈,坝头处流速较大。丁坝B3、B8由于受上游丁坝掩护,坝头流速均小于B1、B6、B10。工况2丁坝B1、B6、B10坝头流速相比工况1分别增大1.71、1.82、1.74 m/s,与主流增大的流速相近;而B3、B8坝头流速分别增大0.93、0.95 m/s。可以发现当流量增大时,直接受主流冲击的丁坝坝头流速变化与主流一致;中间丁坝受上游丁坝的掩护,流速变化比主流小。

受主流直接冲击的丁坝B1、B6、B10,坝头附近流线较为集中,流速较大。其余丁坝根据丁坝间距与坝长的比例分为两组,第一组为B2—B5坝间距是坝长的2倍,第二组为B7—B9坝间距是坝长的3倍,B7—B9坝头整体流速均大于B2—B5,说明坝头流速随丁坝间距增大而增大。当上游来流增大时,主流最大流速区开始向右岸偏移,各丁坝坝头流速增幅不一致,B1、B6、B10增幅较大, B7—B9及B2—B5次之。

4 丁坝群冲刷特性

4.1 丁坝群冲刷状况

通过模拟,对不同工况下丁坝群冲淤产生的地形变化与2017年原始地形进行对比分析,各丁坝具体测点位置及高程变化如图14及表2所示。结合图14和表2可以看出两种工况下各坝坝头发生了不同程度的冲刷或淤积。由图14(b)可以看到在工况3时,从B1处开始形成了一段狭长的冲刷带1,一直延伸至B8坝头处;B10下游也形成了一小段冲刷带2。在工况4时,冲刷带1延伸至B9坝头处,同时冲刷带2开始向上游延伸。通过对比图14(b)

表2 各测点高程对比Table 2 Comparison of height among measuring points

图14 非恒定流丁坝群冲淤变化Fig.14 Change of erosion and deposition of spur dike group in unsteady flow

和图14(c)及表2可发现,工况4冲刷带范围与深度均大于工况3,其中工况3条件下冲刷带1在B7处开始变宽,工况4条件下冲刷带1在B5处开始变宽。通过表2可以发现,在两种工况下丁坝B1所受冲刷均最为严重。由于B1上游没有其他丁坝的存在,坝头处直接受到主流冲击。丁坝B6、B10坝头上游有丁坝掩护,受到部分主流冲击,坝头处也形成较深的冲刷坑。而坝田区各丁坝则发生淤积,B1坝后淤积尤为明显,并且随着流量的增大B1与B2之间的坝田区河床高程有明显抬高。由于行进水流在绕过丁坝流动时分为两部分,一部分水流绕过丁坝流动引起泥沙起动,在坝头附近形成冲刷坑;而一部分水流发生转向,在坝田区形成了漩涡,部分泥沙颗粒在水流的作用下被卷进漩涡从而产生淤积。

4.2 丁坝群冲刷分析

由于丁坝群中丁坝B1、B6、B10 坝头附近冲刷较为严重,因此对这3道典型丁坝坝头附近冲刷情况进行分析。选取冲刷坑沿过流方向的最宽断面为横断面;与之垂直的最长断面为纵断面,对工况3条件下典型丁坝坝头处到达最大冲刷深度时冲刷坑的横纵断面冲刷深度变化进行分析,如图15所示。

图15 工况1典型丁坝坝头冲刷深度变化Fig.15 Variation of scour depth at typical groin head in condition 1

从纵断面图可以看出冲刷坑下游的冲刷明显减弱,这与水流流速、挟沙能力的减小及泥沙颗粒的沉积有关。其中丁坝B1冲刷坑深度为3.08 m,纵、横断面最大宽度分别为29.98、43.97 m;B6冲刷坑深度为2.69 m,纵、横断面最大宽度分别为26.03、35.92 m;B10冲刷坑深度为1.98 m,纵、横断面最大宽度分别为69.87、22.13 m。从上述数据来看,丁坝B1处冲刷最深,范围也较大,综合来说冲刷最为严重。丁坝B6由于受到上游丁坝一定的掩护作用,因此冲刷严重程度次于B1。丁坝B10下游水流能够自由扩散,因此不受壅水的影响,导致B10坝头处沿水流方向形成了一段狭长的冲刷带,但冲刷深度最小。

表3列出了不同工况下3种典型丁坝最大冲刷深度及稳定后冲刷深度。可以发现2种工况下3种丁坝稳定后冲刷深度略小于最大冲刷深度,在工况3条件下,丁坝B1、B6、B10坝头冲刷坑的稳定深度较最大深度减小了0.21、0.27、0.16 m;工况4的分别减小0.19、0.17、0.20 m,说明模拟后期冲刷坑有一定程度的淤积。通过对比2种工况,工况4较工况3的丁坝B1坝头处最大冲刷深度及稳定深度分别增加0.98、1.0 m;B6分别增加0.85 m和0.95 m;B10分别增加0.98 m和0.94 m。

5 结 论

(1)通过对流速特性进行分析可知:受主流直接冲击的丁坝B1、B6、B10,坝头附近流线较为集中,流速较大;其余丁坝根据丁坝间距分为B2—B5(2倍坝长)和B7—B9(3倍坝长)两部分,B7—B9坝头整体流速均大于B2—B5,说明坝头流速随着丁坝间距增大而增大;当上游来流增大时,主流最大流速区开始向右岸偏移,各丁坝坝头流速增幅不一致,增幅大小为B1、B6、B10较大, B7—B9及B2—B5次之。

(2)分析丁坝群进行冲刷特性可知:2种工况中丁坝B1、B6、B10坝头处冲刷均较为严重,冲刷严重程度B1>B6>B10,说明丁坝群中受冲刷的严重程度与是否受主流冲击有关。在实际丁坝群工程中,应当加强丁坝群中受主流冲击丁坝的维护。