丁坝是河流和海岸工程的常见构筑物, 主要用于保护河海岸免遭侵蚀, 改善航道水深条件以维护通航.丁坝会显著改变所处位置附近的水体运动, 引起强三维流动, 造成坝头处较深的局部冲刷坑, 影响坝基稳定和工程安危, 冲刷深度是丁坝设计的重要参数.

水动力条件和冲刷计算模型是丁坝研究的2个主要方面.水动力条件以丁坝周围马蹄涡系或整体流速场为研究对象.早期马蹄涡系发现于桥墩周围, 后来发现丁坝周围也存在马蹄涡系, 许多在桥墩周围流动研究中的有益成果被借鉴到丁坝周围流动的研究中来.烟气法、染料法、荧光粉法、PIV等^[1-3]流动显示方法及数值模拟被用于三维马蹄涡系的视觉呈现.多普勒、ADV等^[4-6]仪器广泛用于三维流速测量, 进而分析时均流动状态.

冲刷计算模型可分为经验方法和半经验模型.因不明确丁坝周围流动或因三维流动过于复杂, 直至现在仍有很多研究者采用经验方法^[7-9].半经验方法有流量束窄法^[10-11], 下降流法^[12], 涡量法和切应力法^[8]等.

虽然不同学者提出的丁坝局部冲刷模型和计算公式大为不同, 但多数公式最终可归纳为一种形式, 即以正挑矩形丁坝在均匀沙床面清水冲刷情况下的局部冲刷深度为基础值, 其他情况下的局部冲刷深度通过对基础值进行水流强度、岸坡、坝型、含沙浓度、挑角和泥沙级配等方面的修正而得到.

本文通过对各种非淹没直立正挑矩形丁坝三维流场的实验数据和可视化研究成果分析后发现：1)水流遇到丁坝后的流动在时均状态下存在明确的不同流动组成部分, 马蹄涡间存在明显的强弱, 可用于简化流动状态；2)坝头处床面附近流速的3个分量中, 横向(坝根向坝头方向)流速占优, 由于切应力与横向流速的平方成正比, 可知坝头处横向切应力对局部冲刷起着重要作用, 但目前尚无研究者明确指出横向切应力的作用并运用于局部冲刷计算.

本文以非淹没正挑矩形丁坝为对象, 分析丁坝周围流动的马蹄涡系图像, 建立考虑横向流速作用的局部冲刷深度的半经验公式.

1 物理原型

本文以非淹没丁坝周围流动情况作为研究对象：正挑矩形丁坝位于矩形槽内, 槽底铺设均匀沙, 远处行进水流平均流速U, 且U小于槽底泥沙的起动流速(见图 1).行进水流遇到丁坝后, 流动主要分为2个趋势, 表层流体上涌绕过丁坝向下游流动, 另一部分流体向下运动形成马蹄涡系.流动造成坝头处形成较大冲刷坑(见图 2).本文采用统一坐标系, 行进水流方向为x, 即纵向；丁坝坝根指向坝头方向为y, 即横向；按照右手螺旋确定z垂直向上；取未冲刷床面、丁坝及与丁坝相交的侧壁面三平面的交点为坐标原点.

2 马蹄涡系

1962年, Schiwind采用烟气法首先研究了层流情况下位于水平底面上的楔形物前的马蹄涡系随雷诺数的变化.随后, Kwan等^{[1-3, 13-14]}采用各种流动显示方法进一步加大加深了对马蹄涡系的探究：从层流扩展到紊流, 从楔形物扩展到圆柱、方体, 障碍物从设置在正中间扩展到位于侧壁, 装置底部由水平面扩展到存在冲刷坑.这一系列的研究结果可以归纳出与本文有关的几点：a)在层流情况下, 随着雷诺数Re的增大, 主涡个数增加, 但目前最多发现为3个, Re数再增大时, 主涡间会出现自我合并；b)在紊流情况下, 床面附近主涡为1~3个, 但在冲刷过程中随冲刷坑的增大, 主涡数量减少；c)当只有2个主涡时, N₂明显弱于N₁(离障碍物最近的主涡编号为1, 向上游侧主涡编号逐渐增大), 而当有3个涡时, N₃明显弱于N₁和N₂, 即存在明显的强弱之分.

马蹄涡是由源源不断的上部来流、无法穿越的障碍物与足够大的流动动能差转化成的压差造成.当主涡N₁形成后, 对于它上游的流体, N₁相当于一个新的障碍物.层流在垂向上流速分布为线性, 由于动能与流速的平方成正比, 当Re数增加时, 上游流体能量的增大造成第2个主涡N₂形成, 但又由于能量差, N₂要弱于N₁.等速度再进一步增大时, 主涡N₃形成.N₂与N₃的形成, 分流了一部分原本流向N₁的流体, 且由于N₂与N₃位于上游, 更容易得到流量补给, 造成N₂和N₃的壮大, N₁一定程度的减弱.当流动从层流向紊流过渡, 流速分布由线性向对数分布过渡, N₂的能量进一步增大到能够穿透N₁的壁垒时, 两涡合并形成新的N₁, 如此循环.但为什么目前发现的主涡数量最多只有3个？估计与一个主涡形成所需的能量及一个涡体冲破另一个涡体壁垒所需能量有关, 雷诺数的增大可导致第4个主涡出现, 但同时又使得N₁和N₂出现合并, 从而造成了主涡个数的稳定.雷诺数越大, N₁与N₂越容易出现合并, 这就意味着在稳定的紊流条件下, 主涡的个数不多于层流情况下的3个.而存在冲刷坑的情况下, 由于越靠近障碍物冲刷深度越大, 从而相对于其他主涡, N₁处于较低位置, 随着冲刷坑深度加大越容易被上游侧的涡体截流合并.从而在考虑存在稳定冲刷坑时的丁坝周围马蹄涡系内的主涡个数时, 本文认为可以采用2个主涡简化(见图 3), N₁作用较强, N₂相对较弱.图 3中其余3个涡N₀为壁涡, N₃和N₄为副涡, 伴随主涡出现.以此5个涡中心点为结点, 相应流动边界上产生8个半鞍点S₀~S₇, 流动内部产生1个鞍点S₈, 该情况也满足Hunt等^[15]提出的结点与鞍点的数量关系要求.

按照马蹄涡系简化流动, 在立面I-I上可将不同高度的行进来流自上而下分成A~D这4个流层, 分别供给马蹄涡系不同部分, 并以a~d表示相应流层的代表迹线.A层流体不参与马蹄涡系, a在流向丁坝过程中上抬, 然后绕过丁坝向下游运动, 在丁坝表面形成半鞍点S₀；B层流体主要供给壁涡N₀, b在与丁坝作用过程中, 贴壁面斜向下向外运动, 在坝面上和冲刷坑底部形成半鞍点S₁和S₂；C层流体主要供给主涡N₁和副涡N₃, c₁绕过涡N₁向床面运动, 在冲刷坑剖面上形成半鞍点S₃和S₄, c₂和c₃缠绕形成主涡N₁, 且沿程不断有同层流体加入；D层流体主要供给主涡N₂和副涡N₄, d₁在冲刷坑底部形成半鞍点S₅和S₆, d₂和d₃缠绕形成主涡N₃, 且在剖面上形成半鞍点S₇.由此构成马蹄涡系的三维图像, 图 4给出了三维图像中的部分迹线.

3 分层流体冲刷作用

B~D层流体分别占据冲刷坑内不同区域sec1~ sec3, 产生了不同的局部冲刷作用.

B层流体因受向外向下的压力梯度作用, 沿丁坝表面斜向下向坝头底部汇聚(见图 5(a)), 流量和流速沿程增大, 在坝头处sec1区域横向流速达到最大, 但因纵向流动区域受限, 纵向流速较小.该层流体在丁坝底部绕过丁坝向下游运动, 从而还存在绕流效应.聚流和绕流作用使坝头处床面最先受到冲刷, 随后受冲刷区域向下游及坝根方向发展.随着冲刷坑加深加大, 因聚流作用增大的横向流速随过水面积的增大而逐渐减小, 降至起动流速时冲刷逐渐停止.绕流作用同样经历逐渐减弱过程, 当冲刷平衡时几可忽略.坝头处B层流体流速尤其是横向流速是影响坝头处冲深的关键因素.

C层流体的局部冲刷作用主要体现在冲刷坑坡面上.冲刷初期流速增大, 在sec2区域内近底形成较大的横向和纵向流速.随著冲刷坑进一步增大, 靠近床面的C层流体缓慢减速, 冲刷逐渐停止.

D层流体在sec3区域内近底纵向和横向流速均较小, 冲刷作用较弱.冲刷坑发展过程中, 部分泥沙因水下休止角小于坡面倾角, 沿坡面下滑被水流带走.

各部分流体即是一个整体又有明确分工, 对计算局部冲刷坑深度而言, 下文将重点关注B层流体.

4 最大冲深作用流量

由第3部分可知, 坝头最大冲刷深度达到的平衡条件应重点考虑B层流体, 即坝头处B层流体流量可视为最大冲深的作用流量.

B层流体位于流动上层, 除壁面附近外受黏性作用较小, 从而可用势流理论来求解这部分水平面流体的流量.根据镜像理论, 可将丁坝绕流看作是无界流体内的垂直平板绕流(见图 6), 由圆柱绕流的保角变换得到.圆柱绕流的复势函数为

$ F\left( \xi \right) = {U_0}\left( {\xi + \frac{{r_0^2}}{\xi }} \right) = \phi + {\rm{i}}\varphi . $

(1)

式中：ξ为复平面坐标；U₀为均匀来流流速；ϕ为势函数；φ为流函数；r₀=L/2为圆柱半径；L为坝长.

采用儒可夫斯基变换

$ Z = \xi + {\left( {{\rm{i}}L/2} \right)^2}/\xi = x + {\rm{i}}y, $

(2)

则Z复平面坐标中垂直平板绕流复势为

$ f\left( Z \right) = \left\{ \begin{array}{l} {U_0}\sqrt {{Z^2} + {L^2}} \;\;\;\;\;\;,x > 0;\\ - {U_0}\sqrt {{Z^2} + {L^2}} \;\;\;\;,x < 0. \end{array} \right. $

(3)

从而流函数为

$ \varphi = {\mathop{\rm Im}\nolimits} \;\;\;\;\left( {f\left( Z \right)} \right). $

(4)

式中：Im(*)表示取复函数虚部.取壁面φ=0, 则φ(Z)即为壁面到Z点处的单位厚度过流流量.这种无界流体内的垂直平板绕流是平面流动, 需适当修正才能符合实际三维流动.

B层流体在冲刷坑上部存在很小的垂向流速, 三维迹线b如图 5所示.图 5中第1条迹线b₁向右偏转是由于受到指向侧壁的压力梯度的作用.需要指出的是, 由于流体斜向下运动, 通过冲深坑的流体只有b₁~b₅并不是b₁~b₆全部(图 5(c)), 对应横向长度为(L-δ_y).δ_y是b层进入冲刷坑的流体在考虑向下运动之前能够到达距离丁坝头部最近的y方向距离, 即b₅迹线在考虑垂向运动前离坝头部的最近距离.

因丁坝迎流面附近B层流体的垂向流速随与迎流面距离的加大而减小, 设δ_x为考虑明显垂向作用时该位置与丁坝表面的纵向距离, 即当距离大于δ_x时流动可视为平面势流, 其流线用公式(4)表示, 从而通过MN段的单位厚度流量为

$ \varphi \left| {_{x = - {\delta _x},y = L - {\delta _y}}} \right. = {\mathop{\rm Im}\nolimits} \left( { - {U_0}\sqrt {{{\left[ { - {\delta _x} + {\rm{i}}\left( {L - {\delta _y}} \right)} \right]}^2} + {L^2}} } \right). $

(5)

δ_x, δ_y与L相比为小量, 略去2阶小量δ_x²、δ_y²和δ_xδ_y得

$ \begin{array}{*{20}{c}} {\varphi \left| {_{x = {\delta _x},y = L - {\delta _y}}} \right. = \sqrt {2L} {U_0}{\mathop{\rm Im}\nolimits} \left( { - \sqrt {{\delta _y} - {\rm{i}}{\delta _x}} } \right) = }\\ {\sqrt {2L} {U_0}{{\left( {\delta _x^2 + \delta _y^2} \right)}^{1/4}}\sin \left( {0.5\arctan \frac{{{\delta _x}}}{{{\delta _y}}}} \right).} \end{array} $

(6)

若设h_B为B层流体的垂向厚度, δ_y=mδ_x, m为比例系数, 从而B层流体对丁坝局部冲刷起主要作用的流量可表示为

$ \begin{array}{l} {q_{V,B}} = \varphi \left| {_{x = - {\delta _x},y = L - {\delta _y}}} \right.{h_B} = \\ \sqrt {2L{\delta _x}} {U_B}{h_B}{\left( {1 + {m^2}} \right)^{1/4}}\sin \left( {0.5\arctan {m^{ - 1}}} \right). \end{array} $

(7)

式中：U_B为B层流体流速.

B层流体对应行进来流中的上层流体, 故U_B值与表面流速相当.若采用对数流速分布, 则U_B与平均流速U和摩阻流速u_*的关系为

$ {U_B} = U + {u_ * }/\kappa . $

(8)

式中：κ为卡门常数.

将式(8)代入(7)可得

$ \begin{array}{l} {q_{V,B}} = \sqrt 2 {\left( {1 + {m^2}} \right)^{1/4}}\sin \left( {0.5\arctan {m^{ - 1}}} \right) \times \\ \;\;\;\;\;\;\;\;\;\;\sqrt {L{\delta _x}} {h_B}\left( {U + {u_ * }/\kappa } \right), \end{array} $

(9)

式(9)是本文得到的第一个重要公式, 表示作用于最大冲刷坑深度所在位置的B流层流量, 简称最大冲深作用流量q_{V, B}.

由式(9)可知, q_{V, B}与B层流体流速(U+u_*/κ)和厚度h_B成正比, 与丁坝长度L和垂向作用距离δ_x的1/2次方成正比.

5 丁坝局部冲深公式 5.1 函数关系推导

q_{V, B}将通过坝头处冲刷坑内断面sec1, 由流量守恒可得

$ {q_{V,B}} = {k_1}{d_{{\rm{se\_b}}}}{L_1}{v_{{\rm{secl}}}}. $

(10)

式中：v_sec1为通过sec1的平均横向流速；d_{se_b}为sec1断面最大冲刷深度；k₁为sec1形状系数；L₁为sec1断面宽度.当冲刷平衡时, v_sec1等于起动流速v_c, 从而

$ \begin{array}{l} \sqrt 2 {\left( {1 + {m^2}} \right)^{1/4}}\sin \left( {0.5\arctan {m^{ - 1}}} \right) \times \\ \;\;\;\;\;\;\sqrt {L{\delta _x}} {h_B}\left( {U + {u_ * }/\kappa } \right) = {k_1}{L_1}{d_{{\rm{se\_b}}}}{v_c}. \end{array} $

(11)

根据实测数据, 坝头处sec1内横向流速的垂向分布如图 7所示, 床面附近αk_s处存在最大值v_ymax, k_s为粗糙高度, α为参数.近底区域横向流速的垂向分布近似以下:

$ \frac{{{v_y}}}{{{v_{y,\max }}}} = {\left( {\frac{{z + {d_{{\rm{se\_b}}}}}}{{\alpha {k_{\rm{s}}}}}} \right)^n},\;\;\;\;\;z \in \left[ { - {d_{{\rm{se\_b}}}}, - {d_{{\rm{se\_b}}}} + \alpha {k_s}} \right], $

(12)

式中：v_{y, max}为近底最大横向流速；n为分布指数；z=0为未冲刷床面位置；z=-d_{se_b}+αk_s为最大横向流速位置.

对近底区域积分得平均横向流速为

$ \overline {{v_y}} = \frac{1}{{\alpha {k_s}}}\int_{ - {d_{{\rm{se\_b}}}}}^{ - {d_{{\rm{se\_b}}}} + \alpha {k_{\rm{s}}}} {{v_0}{{\left( {\frac{{z + {d_{{\rm{se\_b}}}}}}{{\alpha {k_{\rm{s}}}}}} \right)}^n}dz} = \frac{{{v_{y,\max }}}}{{1 + n}}. $

(13)

z=-d_{se_b} +βk_s时可得床面泥沙作用流速

$ {v_b} = {v_{y,\max }}{\left( {\frac{{\beta {k_{\rm{s}}}}}{{\alpha {k_{\rm{s}}}}}} \right)^n} = \left( {1 + n} \right)\overline {{v}}_y {\left( {\frac{\beta }{\alpha }} \right)^n}, $

(14)

式中：β为常数.

坝头附近冲刷坑三维形状类似倒置圆台, 圆台坡度约等于泥沙休止角.从冲刷坑最深处底部水平的角度及泥沙颗粒受力平衡角度出发, 起动底流速可表示为

$ {v_{bc}} = {k_2}\sqrt {\left( {{\rho _s}/\rho - 1} \right)gd} , $

(15)

式中：ρ_s为泥沙密度；ρ为清水密度；g为重力加速度；d为泥沙粒径；k₂为常数.分析认为sec1内横向流速的垂向平均值与$\overline{{{v}_{y}}}$相当, 即冲刷平衡时有${{v}_{c}}={{k}_{3}}\ \overline{{{v}_{y}}}$, k₃为系数.由此可得

$ {v_c} = \frac{{{k_2}{k_3}}}{{\left( {1 + n} \right)}}{\left( {\frac{\alpha }{\beta }} \right)^n}\sqrt {\left( {{\rho _s}/\rho - 1} \right)gd} . $

(16)

将式(16)代入式(11)可得

$ {d_{{\rm{se\_b}}}} = \vartheta \frac{{{h_B}}}{{{L_1}}}\frac{{\sqrt {L{\delta _x}} \left( {U + {u_ * }/\kappa } \right)}}{{\sqrt {\left( {{\rho _s}/\rho - 1} \right)gd} }}. $

(17)

其中

$ \vartheta =\text{ }\sqrt{2}(1+n){{(1+{{m}^{2}})}^{1/4}}\text{sin}\ \text{ }(0.5\text{ arctan }{{\mathit{m}}^{-1}}) $

(β/α)ⁿ/(k₁k₂k₃)为综合常数.鉴于L₁为B流层在冲刷坑内厚度, 与B流层厚度h_B直接相关；δ_x直接影响B层流体在床面上的落点, 与L₁成正比, 可知h_B、L₁和δ_x均与水深h相关, 从而式(17)可写成

$ \frac{{{d_{{\rm{se\_b}}}}}}{h} = f\left( h \right)\sqrt {\frac{L}{h}} \frac{{\left( {U + {u_ * }/\kappa } \right)}}{{\sqrt {\left( {{\rho _s}/\rho - 1} \right)gd} }}. $

(18)

式(18)是本文得到的第2个重要公式, 表示最大冲刷深度与丁坝尺度及流动参数间的关系.可知最大相对冲刷深度d_{se_b}/h与相对坝长L/h的1/2次方, 和无量纲近表层流速 $\left( U+{{u}_{*}}/\kappa \right)/\sqrt{gd}$的1次方成正比.多位研究者^{[8, 16-17]}根据量纲分析得到冲刷深度与坝长的0.46~0.635次方成正比, 本文得到的0.5与之一致.无量纲数f(h)的关系由实测资料确定.

5.2 参数率定与验证

表 1列出了多位研究者^{[7-8, 10, 18-22]}的直立正挑矩形非淹没丁坝在均匀沙(粒径分布的标准差σ_g＜1.5)床面及清水冲刷情况下的局部冲刷深度实验数据, 各参数变化范围为：中值粒径d₅₀为0.26~18 mm, L/h值为1.0~36.47, 流速U为0.183~1.38 m/s, (ρ_s-ρ)/ρ=1.65, σ_g=1~1.46, u_*=0.012~0.173 m/s, d_{se_b}/h值0.53~6.75, U/U_c=0.5~1, 这里U_c为与行进流速和水深相对应的起动流速.表中, 文献[7]的u_*由对数流速分布公式和(gχJ)^0.5计算所得值取几何平均得到(χ为湿周, J为坡度), 其余u_*由对数流速分布公式计算得到.

取表 1内冲刷数据对式(18)进行参数率定, 取f(h)=0.402.式(18)可进一步写为

$ {d_{{\rm{se\_b}}}}/h = 0.402\frac{{\sqrt {L/h} \left( {U + {u_ * }/\kappa } \right)}}{{\sqrt {\left( {{\rho _s}/\rho - 1} \right)gd} }}. $

(19)

将表 1中156个实验室实测数据与采用式(19)得到的计算值比较(见图 8), 各点聚集在45°理论线附近, 可见理论与实测符合良好.

5.3 工程实例运用

万艳春等^[23]给出了北江一个现场正挑丁坝的实测局部冲深数据, 其中丁坝有效长度L= 25 m, 行进水深h=9.3 m, 行进流速U=1.73 m/s, 边坡M=1：2.5, 非均匀沙中值粒径d₅₀=1 mm, 平均粒径d_m=2.36 mm, (ρ_s-ρ)/ρ=1.65, 局部冲刷深度d_se=8.3 m.由对数流速分布公式可得u_*=0.064 m/s, 由沙莫夫公式可得U_c=0.66 m/s, 可知U/U_c＞1.

该情况属于浑水非均匀沙且带边坡的丁坝冲刷情况, 由前言部分可知, 计算此局部冲刷深度需考虑边坡M, 非均匀沙参数σ_g及水流强度U/U_c＞1的影响, 参考常规处理方法, 即

$ {d_{{\rm{se}}}} = {K_M}{K_I}{K_\sigma }{d_{{\rm{se\_b}}}} $

(20)

式中: K_M, K_I, K_σ分别为岸坡、水流强度、泥沙级配对冲深的修正系数.取d=d_m, 并将L, h, U, u_*代入式(20), 同时结合常用的边坡对局部冲深d_se的影响系数$1/\sqrt{{{M}^{2}}+1}$, Lim等^[24]给出的非均匀沙参数σ_g与(d_se+h)_非均匀沙/(d_se+h)_均匀沙的关系图(因文献[23]未给出σ_g数值, 这里取d_m/d₅₀值2.36)和Gill^[7]给出的行进流速大于起动流速造成的冲刷深度影响关系(d_se+h)_U≥Uc/(d_se+h)_U＜Uc=(U_c/U)^6/7, 最后计算得到局部冲刷深度为8.1 m, 与实际测量值偏差为2.3%, 符合良好.

6 结语

丁坝作为河流和海岸工程的常见构筑物已有不少研究成果, 但以往计算模型的建立往往基于一、二维流动或经验方法.本文从三维流动出发, 以5个马蹄涡的流动图像刻画非淹没丁坝周围马蹄涡系的流动状态.通过对实验室及数值模拟三维流动数据的分析, 提出了坝头横向流速在冲刷过程中起关键作用这一新结论.在提出的流动图像的基础上, 针对正挑矩形丁坝周围的流动情况, 通过考虑竖向流速, 修正了横向平板绕流的势流解, 导出最大冲深的作用流量公式(9)；通过流量守恒和冲刷平衡条件, 建立了最大冲刷深度与坝长、行进流速、水深及床面摩阻流速的函数关系, 表明最大相对冲深d_{se_b}/h与相对坝长L/h的1/2次方和无量纲流速$\left( U+{{u}_{*}}/\kappa \right)/\sqrt{gd}$呈线性关系.根据室验数据率定参数, 获得局部冲刷深度计算公式(19), 计算值与实测值符合良好, 表明本文提出的流动图像及依此建立的计算公式具有理论与应用价值.

本文的局部冲深公式基于非淹没正挑矩形丁坝周围的流动图像建立, 而现场丁坝的情况相对复杂, 丁坝挑角、坝型、岸坡等其他参数对水流流动形态有一定影响, 流动图像可能会有所变化, 在后续研究中将进一步考虑由这些因素引起的流动图像变化, 对现有的计算模型进行修正, 以便更符合工程实际.