20世纪50年代末，滦河开始兴建护岸丁坝工程。丁坝具有导流、护岸、防冲和稳定河势的作用，是保滩护岸工程中常见的水工建筑物^[1]。但滦河上的丁坝群建设没有很好规划，丁坝长度、结构形式各异，而且大部分丁坝出现了不同程度的水毁。其中马良子段以险著称，若丁坝失效，滦河主流将直冲小埝，危及马良子村的安全。因此，模拟马良子丁坝群的局部冲刷，并基于丁坝冲刷规律对其优化十分必要。

目前，丁坝的局部冲刷是众多专家学者研究的重点，丁坝建成后由于丁坝上游壅水，会在上游坝根出现一小的立轴回流；至坝头附近，出现下潜流，产生立轴漩涡并相互碰撞、破碎、向下游运移导致坝头局部冲刷剧烈。丁坝下游则出现较大的立轴回流及由此于下游坝根处引发次生反向小回流。这些部位会产生较大的切应力，导致河床发生局部冲刷。喻涛^[2]测量了丁坝横向冲坑的发展过程。Bey等^[3]认为丁坝冲坑的产生是由于丁坝的存在使得丁坝附近湍流对河床底部有冲蚀作用，并认为最大冲深和冲坑形状同样值得研究，因为冲坑形状一定程度上也会影响丁坝水毁失稳。吴学文等^[4]认为丁坝局部冲刷主要是坝头的绕流作用，由此建立了坝头绕流挤压流动物理图示，考虑了非均匀沙的起动问题，得出了非均匀沙河床上丁坝局部冲深公式。事实上，天然河道中的丁坝三维水流情况十分复杂，许多参数难以计量，马永军^[5]分析以往研究结果发现，各专家学者们都只是把回流长度、水流边界、冲深直接联系起来，并没有找出各参数之间的关系。苏伟等^[6]对不同结构形式丁坝进行物理模型试验，发现丁坝的结构形式对水毁影响明显。计算机技术的发展，为复杂水流问题的求解提供了计算手段，Tingsanchahi等^[7]利用二维平均水深模型并引入校正因子改善k-ε模型，以求解河床底部应力。Molls等^[8]通过结合恒定的涡黏性湍流模型改进了不稳定二维平均水深模型。Jia等^[9]利用二维平均水深模型计算丁坝附近水流流线和流速。Duan等^[10]用改进二维模型对天然河道中泥沙输运过程进行仿真。结果表明，二维模型只能近似计算泥沙水平运动，无法较好地模拟泥沙在河道中的运输。越来越多的学者发现丁坝附近水流表现出很复杂的三维特性^[11]，因此，开始丁坝的三维数值模拟研究。Nagata等^[11]利用三维水利模型计算了雷诺平均方程和斯托克斯方程，成功模拟了丁坝周围的泥沙输运过程。Ouillon等^[12]利用三维k-ε模型研究了丁坝周围三维水流结构和自由液面形状，求解了坝后冲刷区水流的水力特性。杨兰等^[13]通过使用Flow-3D软件，对上挑丁坝群附近流场及局部冲刷进行了三维数值模拟，并得出了可适用于计算丁坝群的湍流模型和推移质输沙率模型。但目前对丁坝群的三维数值模拟或物理模型模拟大部分针对的是直槽式水槽，而非实际工程中复杂的天然河道。本文首先基于不同丁坝结构形式的三维水槽模型模拟，对滦河上各种丁坝的冲刷规律进行探索，作为调整丁坝的依据。然后对天然河道的丁坝群进行三维数值模拟，对比分析冲刷数据与实测地形，其结果可为工程实际提供一定的参考。

1.   理论模型

1.1   控制方程

丁坝群流场三维流动控制方程为连续性方程和动量方程：

式中：u_i为i方向的分速度；u′_i为i方向的脉动速度；P为压力；S_ij为应变率张量；$\overline {u' {_i}{u' {_j}}}$为雷诺应力张量；ρ为流体密度；v为动力黏度；v_t为湍流黏度；k为湍动能；δ_ij为克罗内克符号(δ_ij=1，i=j，δ_ij=0，i≠j)。

1.2   RNG k-ε方程

标准k-ε模型用于强旋转流或带有弯曲壁面流动时容易出现失真。RNG k-ε模型通过修正湍动黏性系数和在ε方程中增加了反映主流的时均应变率^[7]，可以很好地模拟高应变率和流线弯曲程度较大的流动，而天然河道的复杂地形以及丁坝与水流之间的相互作用会带来剧烈的水流变形和破碎，且许多专家学者认为RNG k-ε模型对丁坝得到的结论更为可靠。

RNG k-ε控制方程：

式中：k为湍动能，$k = \dfrac{1}{2}\overline {{u' _i}{u' _j}} $；ε为湍动能耗散，$\varepsilon = v\overline {\left( {\dfrac{{\partial {u_t}' }}{{\partial {x_k}}}} \right)\left( {\dfrac{{\partial {u_t}' }}{{\partial {x_k}}}} \right)} $；μ_t为涡黏性系数，μ_t=C_μρk²/ε；τ_ij为雷诺应力，${\tau _{ij}} = \dfrac{{{\mu _{\rm{t}}}}}{\rho }\left( {\dfrac{{\partial {u_i}}}{{\partial {x_j}}} + \dfrac{{\partial {u_j}}}{{\partial {x_i}}}} \right) - \dfrac{2}{3}\rho k{\delta _{ij}}$；模型常数C_μ=0.085，C_1ε=1.42，C_2ε=1.68，σ_k=σ_ε=0.717 9。

1.3   泥沙输运模型

目前，研究丁坝局部冲刷的规律主要有5个理论途径，分别是流速、切应力、能量平衡、统计法则及沙波运动，而湍流流场和泥沙输运对于丁坝群局部冲刷坑的发育影响在理论层面还不是十分成熟。基于前人对丁坝冲刷规律的试验以及泥沙运动理论分析，选用以希尔兹(Shields)数为基础的泥沙推移质输沙率模型作为此次丁坝群冲刷三维数值模拟的理论模型。

与杨兰等^[13]的泥沙输运模型相同，模型计算中不考虑悬移质对冲刷坑发育的影响，只分析河床底部推移质的冲刷作用，并认为推移质为均匀的泥沙颗粒。采用的推移质输沙率公式以切应力为主要参数，且切应力与输沙率的关系成正比。

水流切应力:

泥沙起动的临界切应力:

希尔兹数:

临界希尔兹数:

式中：U_*为摩阻流速；U_*c为临界摩阻流速；ρ和ρ_s分别为水流和泥沙密度；d为泥沙平均粒径。

假定单宽推移质输沙率计算式为：

式中：q_b为单宽推移质输沙率；u_b为推移质的平均输运速度；p为泥沙起动概率。对于实际工程地质情况的土工试验，可以根据上式得到u_b和p，即可求解q_b。

动摩擦力为:

F_D=f_D并联立(7)，得到:

式中：aU_*为推移质运动的水流速度，当离沙床较近时，a=6～10；C_D为推力系数。当θ₀=0时，u_b=0，θ₀相当于止动相对切应力，应小于临界相对切应力θ_c。由此，式(15)可写为:

取泥沙临界切应力和运动颗粒所受切应力之和为切应力，即:

式中：n为泥沙颗粒总数，而单位面积河床泥沙颗粒总数1/d²与p的关系为：p=n/d^-2。联立式(9)和(10)，得到床面泥沙颗粒的起动概率为

综上所述，影响单宽推移质输沙率的两个重要因素为u_b和p，将u_b和p的表达式代入(11)，并取动摩擦系数β=0.08，常数a=9.3。可得如下推移质输沙公式：

式(19)即为所采用的泥沙运动模型。可以发现，对于某些泥沙颗粒，推移质输沙率只与希尔兹数有关，只有当希尔兹数大于临界希尔兹数时，泥沙才会起动被冲刷，采用θ_c=0.05。

2.   模型验证

实际工程丁坝段如图1所示，可见11#丁坝相比周围的丁坝明显较长，不够合理。这是由于丁坝有束窄河床的作用，河床被束窄后，河道断面面积减小；虽然丁坝前有一定程度的壅水，但断面面积仍有大比例收缩，导致相同流量下流速迅速增大，特别是河道底部的切应力会迅速增大，下切河床造成冲刷。同时，该河段自9#丁坝后河道宽度明显缩窄，流速更大，由于滦河段河床地质条件差，基本是粉细砂组成，更容易遭到冲刷。

图  1  马良子段丁坝群布置

Figure  1.  Layout of spur dike group along Maliangzi reach

下载: 全尺寸图片 幻灯片

对比河底实测地形图发现，11#丁坝后有很深的冲刷坑，最深处超过13 m，因此选择11#丁坝作为主要的研究对象，研究不同丁坝长度对丁坝群的冲刷影响。

计算模型长1 400 m，宽1 000 m，高30 m，其中11#丁坝上游400 m至下游300 m为动床计算区域，其余为定床。河床主槽糙率为0.02，平均水深为5 m，为非淹没丁坝群，计算丁坝群位于左岸，计算河道模型平面轮廓与天然河道一致，考虑到河道宽度远远大于水深，因此河道断面形状取矩形。原始11#丁坝长为34.5 m，动床所铺泥沙厚度为20 m，泥沙平均粒径d=0.018 mm，密度ρ_s=2 650 kg/m³。计算区域采用矩形网格划分为A，B，C区，A区和C区为定床，B区为动床。对动床部分加密网格，总网格数量约133万，网格分块划分见图2。

图  2  计算区域网格剖分

Figure  2.  Meshes of computational domain

下载: 全尺寸图片 幻灯片

计算时采用有限体积法离散控制方程，湍流模型为RNG k-ε模型。上游设为速度进口，根据上游实测流速和平均深度，给定模型流量为8 430 m³/s，平均水深5 m，出口设为自由出流，上表面采用VOF法捕捉液体自由面，固体边界为无滑移面，近壁面采用标准壁面函数进行处理。对比11#坝头附近的冲刷深度与实测地形图，实测与模拟冲刷地形拟合较好，如图3所示。可见，在复杂天然河道形状下，可以利用RNG k-ε模型和推移质输沙率模型模拟丁坝的局部冲刷。

图  3  坝头附近冲刷深度横纵断面对比

Figure  3.  Comparison of transverse and longitudinal sections of the scouring depth near the spur dike head

下载: 全尺寸图片 幻灯片

3.   基于丁坝坝长的调整方案

滦河上马良子段丁坝结构形式简单，多为下挑式丁坝，挑角约70°~85°，变化幅度不大，且坝头形式单一，为直立式丁坝。通过丁坝设计资料以及地形勘察资料发现，马良子段11#丁坝设计明显过长，设计长度34.5 m，比相邻丁坝长约10 m，而且10#~12#丁坝在河流收缩段。测量数据也表明11#丁坝坝后是冲刷最严重的区域，最大冲深达13.98 m。因此丁坝长度不合理很可能是造成严重冲刷的原因。本文将原坝长方案定为方案1，将11#丁坝缩短5 m和10 m，对应方案2和3，对丁坝周围水力特性及冲刷特性进行研究。

3.1   丁坝绕流流场与流速分布

方案1~3计算域水流流场与流速如图4所示。可见，随着11#丁坝坝长的减小，水流随之平顺，丁坝后端坝田附近形成的漩涡减小。同时，在原有丁坝群间距不变的情况下，缩短11#丁坝依然可以保持丁坝群对水流的有效影响，即坝前水流流入坝后坝田区与坝田内水流相互作用，与岸边不直接接触，主河槽流速较两岸大，方案1中过11#坝后主河槽流速明显增加，方案2和3增速不明显。

图  4  11#丁坝附近绕流流场及流速分布

Figure  4.  Flow field and velocity distribution around No.11 spur dike

下载: 全尺寸图片 幻灯片

3.2   冲刷分析

丁坝群冲刷形态对比如图5所示，可以发现调整方案1和2中丁坝群最大冲深坑出现在11#丁坝坝头附近，并且在主流区形成一狭长冲刷带；而方案3丁坝群则没有明显冲刷坑，冲刷深度较为一致。这说明缩短丁坝群中过长的丁坝可以有效减轻其坝头的局部冲刷，丁坝群调整后河床冲刷明显减轻。

图  5  丁坝群冲刷形态对比

Figure  5.  Comparison of scour patterns around spur dikes

下载: 全尺寸图片 幻灯片

根据实测冲刷地形在11#丁坝坝头附近设置测点，其中测点1是原坝后冲坑最深处，测点2和4分别是测点1上、下游100 m位置，测点3为测点1的右边5 m处，测点布置如图6所示。丁坝群冲刷数值对比如表1所示。可见，缩短过长的丁坝长度可以有效减小丁坝局部冲刷深度，从而降低丁坝的水毁风险。

表  1  3种方案各测点数值对比

Table  1.  Comparison of values of measuring points in three schemes

	测点	水面高度/m	流速/(m·s-1)	冲深/m	冲深减小值/m
方案1	1	25.45	6.40	13.98
	2	26.20	6.00	7.00
	3	25.63	6.49	7.00
	4	25.30	6.50	10.70
方案2	1	25.20	6.90	11.00	2.98
	2	26.00	5.60	6.50	0.50
	3	25.80	6.60	7.80	0.80
	4	25.90	6.50	9.50	1.20
方案3	1	25.90	5.40	5.30	8.68
	2	26.00	3.50	5.00	2.00
	3	25.80	5.00	5.40	1.60
	4	25.80	4.00	4.70	6.00

下载: 导出CSV 

| 显示表格

图  6  测点布置

Figure  6.  Layout of measuring points

下载: 全尺寸图片 幻灯片

4.   结　语

基于天然河道形状，利用推移质输沙率模型和VOF自由表面处理法，对非淹没丁坝群的绕流流场、流速、局部冲刷进行三维数值模拟研究，并与实测资料进行对比，得到以下主要结论：丁坝群中某一丁坝长度较周围丁坝明显偏长易引起河道水流集中，坝后漩涡加剧，且丁坝长度愈长，阻水能力愈大，坝头冲刷越剧烈。最大冲刷发生在坝头处，且在冲坑内略偏坝轴线下游，上游部分坡度较陡，范围较小，下游部分坡度较缓，范围较大。说明丁坝群设计时应合理规划丁坝长度，丁坝群内各个丁坝长度不协调可能会形成新的水害。

由于实际工程中丁坝的冲刷坑形成非常复杂，本次模拟的冲刷坑下游冲刷与实际还有一些偏差，对丁坝局部冲刷机理和实际水流过程有待进一步研究。