河道整治工程主要包括险工、控导及护滩工程等，其建筑物形式为丁坝、垛、护岸，一般简称坝垛，坝垛由土坝基、护坡和护根三部分组成。护根亦称根石，它是坝垛工程的基础，护根的稳定状况对工程的稳定起着决定性作用。黄河河道整治工程中险工坝岸常用散抛石护根，以维护坝、垛与河岸的安全，坝岸工程稳定与否直接影响到堤防安全，而坝岸根石的稳定与完整程度又直接决定坝岸工程的稳定程度，因此，研究险工坝岸根石走失规律、预测根石走失情况，对于黄河堤防工程设计与防汛抢险具有重要的现实意义。

目前针对根石走失的研究还较少，根石走失分析类似于块体在水流作用下运动，水下运动块体的分析常采用Evett等^[1]提出的水流作用力公式，将淹没在水流中块体所受水流作用力分解为水平推移力和垂直水流方向的上举力；庞启秀等^[2]通过变坡水槽试验测定不同流速强度下不同尺寸块体所受的拖曳力和上举力，提出了含有形状影响的拖曳力计算公式；缑元有^[3]利用Evett的水流作用公式，考虑黄河河道整治工程冲刷坑坡度、水流含沙量、冲刷流速和河床条件等因素，推导出不同含沙量、水深及流速条件下根石起动尺寸，与实测结果基本相同；胡宝琳等^[4]运用PFC2D软件建立了抛石潜堤模型，模拟分析了抛石潜堤的破坏过程，并总结了抛石重量和堤顶宽度等对堤岸失稳的影响；姚仕明等^[5]对护岸抛石开展了水槽试验，解释了抛石护岸工程破坏的原因和形式，并对实际工程中抛石方法做了总结。以上公式推导和物理模型试验方法可得出护岸局部块石的运动及防护措施效果。数值模拟作为一种新的研究手段已越来越多地应用于水利科研中，但目前关于根石走失的数值模拟研究很少，将数值模拟计算结果与实测数据相互对比验证是确定工程安全的最可靠方法，因此有必要建立相关的数值模型对根石走失进行细观模拟。

根石稳定性分析是典型的非连续介质滑移及稳定平衡过程，这一过程适合采用离散单元方法进行模拟，并且还需考虑水流边界对根石走失的影响，基于此本文选用Flow-3D软件对根石走失过程的数值模拟进行初探，为根石走失的数值模拟研究提供技术参考。

1 Flow-3D介绍及数值模型建立

Flow-3D是近年来应用较多的计算流体动力学软件，软件中GMO(通用运动物体)模型可计算物体在水流作用下的运动，物体在运动过程中可拥有6个自由度，物体运动的形式有指定运动与流固耦合运动，流固耦合运动为运动物体与水流相互影响、相互作用，其中物体在运动过程中的相互碰撞也被充分考虑，适用于根石走失的数值模拟。基于N-S水流控制方程，采用软件中的RNG k-ε紊流模型^[6]，该紊流模型经许多学者验证对复杂水流形态模拟精准^[7]，在块体绕流上应用性也很强^[8]。水流液面的描述采用VOF法^[9]，对水流作用下运动的根石采用GMO模型计算，软件独特的FAVOR技术^[10]可对运动的根石进行准确描述，基于上述方法建立根石走失的数值计算模型。

2 数值模型算例验证

根石在水流作用下受到的力有拖曳力、上举力、漩涡力及脉动压力等，其中，水流拖曳力对根石整体的影响较大，本次针对块体在水中的拖曳力进行分析。

2.1 数值模型的建立

图 1为数值试验模型布置，水流冲击置于水槽中部的5 cm立方块体，水槽尺寸为2.0 m×0.5 m×0.3 m，采用200 000个结构化矩形网格对计算区域进行划分，并对块体区域局部加密。水槽中初始水深为0.15 m，入口边界处设定水深为0.15 m的水流，水流流速工况见表 1；出口处设定自由出流边界；侧边和底面均设定为固壁边界；计算区域上方为1个标准大气压。

2.2 拖曳力结果对比分析

对拖曳力的计算采用Evett^[1]提出的公式：

${F_{\rm{D}}} = {C_{\rm{D}}}\frac{1}{2}\rho A{u^2}$

(1)

式中：F_D为块体在水中的拖曳力；C_D为块体拖曳力系数；ρ为流体密度；A为块体的迎流面积；u为水流作用的垂线平均流速。根据Evett^[1]给出的结论，块体形状雷诺数Re>1 000时，当宽高比为1时, C_D取1.05。

数值模型可直接输出作用于物体上的拖曳力，表 1给出了不同流速下拖曳力F_D的公式计算值与数值模拟值，图 2为拟合曲线，由表 1与图 2可以看出，拖曳力的公式计算值与数值模拟值整体拟合效果较好，随着水流流速的增大，拖曳力的公式计算值略微大于数值模拟值，根据Krueger等^[11]进行的拖曳力测量试验，结果显示物体的拖曳力系数C_D随着块体形状雷诺数的增大逐渐减小，如图 2所示，块体形状雷诺数随着水流流速逐渐增大而增大，因而流速逐渐变大后拖曳力的数值计算值会略小于采用固定C_D的公式计算值。拖曳力的公式计算值与数值模拟值对比整体接近，说明数值模型计算水流拖曳力准确度高。

3 工程实例分析

舜帝控导下沿工程位于黄河永济段，工程含13个坝垛，全长1 000 m，设计坝顶高程342.40~342.10 m，平均坝高3.90 m，设计按2020年设计洪水流量4 000 m³/s设防，其主要作用为约束主流摆动范围、护滩保堤，下游坝垛均为险工，散抛根石设计，平面布置及横断面如图 3；工程采用坝、垛结合的平面型式，共设12个垛坝，坝垛间距80 m，迎水面与连坝轴线交角30°，背水面与连坝轴线交角60°，垛头采用半径15 m的圆弧衔接，连坝土体顶宽10 m，垛及连坝迎水面采用块石全裹护。

3.1 实际工程数值模型的建立

实际工程中坝垛段较多，且形状基本一致，选取2个坝垛段进行分析，如图 4所示，按照比例1:1建立数值模型并导入软件中，本次研究为单排散抛根石的走失，根石的布置围绕坝垛均匀排列，坝垛周围分为迎水面、上跨角、坝前头和下跨角区域，对4个区域的根石采用不同颜色标记，迎水面采用黑色，上跨角为红色，坝前头为蓝色，下跨角为粉色，以便对走失后的根石进行观察，并在4个区域的中部设定监测点如图 4，监测计算过程中水流流速变化。根石尺寸为边长50 cm的立方体石块，密度为2 600 kg/m³，每个块体均设置为有6个自由度的碰撞物体，由于石块碰撞之间含有水体，取碰撞恢复系数为0.1，摩擦系数参考抛石基底与粉砂地基土取0.4^[12]。计算河道宽30 m，长182 m，网格划分采用3个连接式网格块，如图 5所示，网格块1和网格块3采用0.50 m×0.50 m×0.50 m的网格，网格块2采用0.25 m×0.25 m×0.25 m的网格，总网格数约为40万个。

根据工程施工水位339.98 m和设计水位341.32 m，取相对水位高度计算。设定计算区域初始水位0.5 m，入口边界水深在180 s内由0.5 m逐渐上升至1.84 m随后保持在1.84 m；出口处设定相应的尾墩控制出流量；侧边和底面均设定为固壁边界；计算区域上方设定一个标准大气压。计算时间为300 s。

3.2 根石走失分析

图 6给出了计算结果中不同时间点的根石走失情况，图 7给出了监测点流速变化。由图 6和7可知，70 s时各处水流流速较小，坝前头根石有微小移动，其余区域根石静止无变化；120 s时由于上游边界处的水位加高，经过坝前头的水流流速增大，部分根石出现明显移动，其中坝前头中部根石开始脱离坝身；150 s时坝前头根石出现大位移，水流流速在2.0 m/s左右时根石起动并顺着挑流主流向下游滚动，迎水面与上跨角处部分根石也出现小位移并相互碰撞；180 s时坝前头根石在水流作用下基本全部走失，大部分走失根石堆积在坝前头与下跨角连接处区域，少数根石走失较远，上跨角处靠近坝前头的根石开始沿着坝身移动；240 s时由于坝前头处根石走失严重，上跨角根石在移动过程中没有根石的阻挡作用，且坝前头水流流速较大，最高达3.6 m/s，因此经过坝前头处根石的运动速度较快，部分根石走失严重；300 s计算结束，此时坝前头处无根石，走失的根石大部分堆积在坝垛裆间，坝前头与上跨角的部分根石走失较远，迎水面与下跨角处的部分根石有较小位移，但在整个过程中并未脱离坝身，受水流影响较小。

3.3 水流拖曳力分析

图 8给出了计算过程中不同区域位置中间根石的拖曳力值变化，可见，迎水面与下跨角两处水流对根石的拖曳力最小，结合图 7可知两处的水流流速最小，并且根石并无较大位移，因此拖曳力的大小随着水流流速的变化趋势改变；上跨角处中间根石在前期拖曳力大小变化与水流流速变化规律相似，结合图 6根石走失情况可知，240 s后该处根石开始沿坝身走失，经过坝前头时水流流速突然增大，根石快速移动并伴有滚动现象，水流对根石的拖曳力大小变化较大，经过坝前头后拖曳力下降；坝前头中部根石在整个计算过程中所受拖曳力较大，由于该处水流流速最早达到根石的起动流速，因而前期拖曳力上升最快，180 s后该处根石已脱离坝身进入坝裆间，后期所受拖曳力较为平稳。

4 结语

本文基于Flow-3D软件建立了根石走失特性研究的数值模型，对实际工程坝垛进行数学建模并计算，计算结果中对根石走失情况与根石拖曳力进行分析，得到以下结论：

(1) 数值模拟结果中根石的走失规律与实际工程相似，该数值模型可定性分析根石走失的运动规律，定量分析根石走失过程中所受拖曳力的变化。

(2) 整个计算过程中坝前头处流速最大，上跨角处流速次之，此两处水流对护岸根石影响较大，其中坝前头处根石起动的流速为2.0 m/s左右，洪水时该处流速最大可达3.6 m/s左右。迎水面与下跨角处流速较小，水流对该两处的根石影响小。

(3) 数值模拟结果显示上跨角与坝前头处少部分根石走失较远，大部分走失的根石集中在坝裆间，迎水面与下跨角处根石基本无走失现象。

(4) 洪水作用下边长为50 cm的正方块体在坝前头与上跨角处走失严重，不能很好地起到保护作用，实际工程中在该两处应采用质量更大的块体护根。