当卡水电站下游局部冲刷三维数值模拟研究

，，　，，

(1.西北农林科技大学 水利与建筑工程学院，陕西 杨凌　712100；2.云南水利水电建设工程技术开发有限公司，昆明　650021；3.渭南市东雷二期抽黄工程管理局，陕西 渭南　714000)

1　研究背景

水利水电工程闸坝下游水流流速高、紊动强烈，易造成下游河道产生严重的局部冲刷及沿程冲刷，对水工建筑物自身安全及河道两岸形成严重威胁[1- 2]，对闸坝下游局部冲刷准确预测预报不仅可以避免巨大经济损失，也可为下游城镇、居民、生态环境和企业提供安全保障。因此，研究预测水工建筑物闸坝下游局部冲刷问题具有重大的现实意义，也是当今水利专家难以解决的重要问题[3-6]。

在计算机出现以前，一些学者基于实际观测资料提出了若干闸坝下游冲刷坑深度经验计算公式，如Lacey公式、Kennedy公式、维兹果公式、加切奇拉捷及毛昶熙公式等，但由于这些公式大多基于作者自己观测资料分析得到，公式的结构形式及计算结果差别较大，公式应用范围也受到较大限制[3]。因此，对于复杂的水力学问题，常常需要借助于物理模型进行试验研究，但物理模型试验也具有较大的局限性，常常受到时间、空间、观测手段及研究经费等因素的限制。随着数值分析理论和计算机软硬件的迅猛发展，数值模拟技术在复杂问题的研究中表现出明显的优越性，因而被广泛应用于实际问题的研究中。

王晓松等[7]利用三维标准k-ε模型模拟计算了挑流冲刷平衡时的水流流态，模拟中考虑了复杂的冲坑底部形态和自由液面对紊流流场的影响，得到的水流流态、流场及壁面剪切区与实测情况基本符合；郭红浩等[8]、高改玉[3]、史志鹏[4]分别研究了中低水头水电站泄水建筑物泄流流量、流场、水位及压力分布，均取得与实测情况符合良好的计算结果；杨建明等[9]利用动网格技术模拟计算了挑流冲刷坑二维水流流场，实现了动网格边界随着冲刷坑的发展而更新移动，冲刷坑中流场信息也随之变化。不足之处在于，为了使网格更新过程具有较好的质量，作者规定了河床形态的曲线形式；网格更新与时间步长不完全同步，而是在若干计算时段后更新。劳尔平[10]对河道建筑物下游冲刷过程进行了二维数值模拟，模拟中采用阶梯形动边界逼近冲刷坑平衡状态，以避免网格发生畸变；邓军等[11]在贴体非正交曲线坐标系中，采用二维数学模型模拟计算了二维散粒体泥沙床面冲刷情况，并基于物理模型试验观测对数值模拟结果进行了验证，结果表明两者符合较好；李艳富[12]基于贴体非正交曲线坐标，建立以冲刷坑中动水压力和冲刷深度为判定条件，通过动网格技术使网格边界随之移动变化，实现了对二维无黏性河床挑流冲刷坑的模拟计算；高改玉[3]利用钱宁斜坡上泥沙起动临界剪切应力公式计算值作为边界网格节点移动的依据，对Fluent软件进行二次开发，对横丹江水电站泄水建筑物下游冲刷过程进行了二维数值模拟，得到了较好的模拟结果；凌建明等[13]采用三维标准k-ε模型模拟计算了桥墩附近绕流冲刷坑中水流流态及三维流场，得到了冲刷坑中水流剪切应力在床面上的分布情况，为分析桥墩附近绕流冲刷奠定了基础；韦雁机等[14]采用k-ε紊流模型，对于水面边界采用了刚盖假定，基于OpenFOAM开源软件模拟计算了圆柱体周围的紊动流场，利用经验公式计算该时刻的推移质输沙率，并利用动网格技术处理床面地形随时间的变化，对于斜坡上的泥沙起动条件，采用了钱宁等人的研究成果，充分考虑了坡度的影响。此外， Olsen[15],Ushijima[16],Nobuhisa等[17],祝志文等[18]在进行桥墩周围局部冲刷三维数值模拟时，也具有类似的思考；熊文等[19]基于Fluent软件动网格自动更新技术，采用刚盖假定追踪自由液面，将床面瞬时剪切应力作为泥沙起动机运输的水动力学条件，计算出不同方向上床底泥沙单宽体积输沙率，以此为基础得到床面网格坐标的瞬时变化，从而对桥墩冲刷以及周边流场三维形态进行了全过程动态跟踪分析，为避免床面坡度接近泥沙休止角时造成的输沙率数值计算的异常，采用了Chen等[20]改进的泥沙起动临界剪切应力公式。

综上所述，诸多学者对桥墩周围冲刷坑水流流态以及局部冲刷过程的三维数值模拟进行了较多的研究，但大多采用刚盖假定来追踪自由液面。而对于闸坝下游冲刷坑冲刷过程的数值模拟研究还局限在二维，不能够反映出冲刷过程冲刷坑内强烈的三维特性。此外，与桥墩周围局部冲刷相比，闸坝下游冲刷过程中自由液面变化比较剧烈，采用刚盖假定进行追踪已经不能满足实际的需求。

鉴于此，本文采用VOF水-气界面追踪方法来处理闸坝下游冲刷坑内复杂的自由液面，利用UDF宏函数对Fluent软件进行二次开发，通过采用Chen等[20]改进的泥沙起动临界剪切应力计算公式计算值作为网格移动的判别标准，对当卡水电站泄水闸下游冲刷坑的冲刷过程进行了三维数值模拟，得到了冲刷坑内冲刷坑深度及形态。

2　河床变形模型

2.1　平坡上泥沙起动

根据Shields在1936年推导出的泥沙起动临界剪切应力式(1)判断平坡上泥沙是否起动，即当数值计算得到的剪切应力大于临界起动剪切应力时，泥沙起动，此时

τb,cr=ρg(s-1)d50θcr。

(1)

式中：s为相对密度，s=ρs/ρ，ρs为泥沙颗粒密度；ρ为水流密度；g为重力加速度 ;d50为泥沙颗粒中值粒径；θcr为临界Shields数。根据Shields曲线分段拟合得到的表达式为[21]

(2)

式中D*=d50[s-1g/υ2]1/3,υ为水的运动黏度。

2.2　斜坡上的泥沙起动

在泄水闸下游冲刷过程中，冲刷坑内各个部位都有泥沙起动。泥沙起动边界条件可以归纳为底部平整床面、侧面边坡、背水坡(正坡)和迎水坡(反坡)。除底部床面具有较小范围的平整外，冲刷坑四周均有一定的边坡，因此，冲刷坑边坡上的泥沙起动条件采用斜坡上的泥沙起动计算公式。对于斜坡上临界泥沙起动剪切应力计算公式的选取，高改玉[3]、祝志文等[18]、韦雁机等[14]均采用钱宁在希尔兹公式的基础上考虑坡度及水流方向的影响而得到的斜坡上泥沙起动剪切应力普遍式。但是，当床面坡度等于泥沙休止角时，采用钱宁公式计算出来的泥沙临界起动切应力会趋近于0[18]，无量纲临界起动剪切应力计算公式为

(3)

(4)

(5)

图1　斜坡上泥沙颗粒受力示意图Fig.1　Sketch of forces acting on a sedimentparticle on sloping bed

2.3　推移质输沙率

当斜坡上的泥沙开始起动以后，输沙率是用来描述泄水闸下游冲刷坑形态发展的一个关键要素。由于试验中选用的泥沙颗粒粒径较大，这里认为推移质运动是引起床面冲刷的主要因素，在计算中不再考虑悬移质的影响。这里采用数值计算中广泛使用的van Rijn[23]输沙率计算公式，即

(6)

(7)

(8)

2.4　河床变形的计算

当判别出网格节点是否需要移动后，接着需要解决的问题是确定网格节点随着单个时间步长进行移动的位移量，在这里采用泥沙质量守恒方程[24]来进行计算，即

(9)

式中：h为动床床面高程；t为时间变量；n为动床泥沙颗粒间的孔隙率，取为0.40[25]；qbx和qby则分别采用式(7)、式(8)进行计算。

3　泄水闸下游三维冲刷数值模拟

3.1　控制方程[4]

连续性方程为

(10)

动量方程为

(11)

湍动能κ方程为

(12)

(13)

式中：t为时间;i,j=1, 2, 3;ui为xi方向的速度分量；ρ1为体积分数加权平均的密度；p为修正压力；μ为体积分数加权平均的分子黏性系数；ακ,α分别为湍动能κ、耗散率下的普朗特数，ακ=α=1.39；μt为湍流黏性系数，为平均梯度引起的湍动能产生项。

3.2　工程概况

当卡水电站是玉树地震灾后恢复重建项目——西航当代水电站移址重建工程，为子曲河流域规划中的第4座梯级水电站。坝址位于玉树县下拉秀乡当卡村子曲河中游，距下拉秀乡20 km。当卡水电站主要任务为发电，总装机容量为12 MW，电站总库容为1 109万m3，正常蓄水位为3 848.00 m。

水电站枢纽主要由单孔泄洪闸、溢流坝、右岸重力挡水副坝、左岸砂砾石坝、河床式厂房等组成。经过模型试验，优化后推荐方案的泄洪闸闸室段上、下游方向长32 m，闸孔净宽为7.0 m，闸室前端底板为平底，高程为3 824.5 m，闸室末端以弧段与消力池底板衔接，消力池段长48 m，底板高程为3 818.5 m，消力池段净宽为11 m，其纵剖面图见图2[26]。

图2　泄洪闸纵剖面图Fig.2　Longitudinal profile of flow dischargeconstruction

王新强[27]为了研究流量和泥沙粒径对冲刷坑形态的影响，在参考当卡水电站水工模型试验的基础上，进行了16个组次的试验，分别观测记录了冲刷坑的形态，本文选取组次分别为6，10，14的试验数据进行数值模拟分析。

3.3　计算域网格及边界条件

运用AutoCAD软件建立当卡水电站下游冲刷坑模型区域图，该模型计算区域总长226.5 m，见图3，分为动区域和静区域。动区域为动床，长140 m；静区域为定床，长86.5 m。

图3　数学模型计算区域及边界示意图Fig.3　Computational domain and boundary

采用有限体积法对水流的三维空间形态进行数值模拟，流场网格的正确划分是关键，利用GAMBIT网格划分软件可以生成结构化的六面体网格以及非结构化的四面体网格。一般而言，结构化的六面体网格划分生成的网格质量较高，在同样网格尺寸情况下，单元数量比非结构化的四面体网格较少，计算时间相对缩短，是精细化冲刷网格划分的首选，但是，对于结构形体较为复杂，且对网格有较大变形要求的情况下，非结构化的四面体网格具有较好的优势。因此，为了满足本文模拟的泄水闸下游冲刷坑形态变化比较大的特点，对于动区域部分采用的是非结构化的四面体网格进行划分;同时为了减少网格数量，保证水流的精细模拟,对于静区域部分采用的则是结构化的六面体网格进行划分。数学模型中网格划分的间距为0.50～0.75 m，网格总数为345万。

3.4　动网格技术的运用

研究中动网格自动更新技术需要用到Fluent中的UDF宏的功能，通过利用DEFINE_GRID_MOTION( )宏命令来控制动床边界上每个网格节点随着时间步长(time step)的移动。为了准确模拟泄水闸下游冲刷坑内实时冲刷过程，利用宏NV_MAG(F_STORAGE_R_3V(f,t,SV_WALL_SHEAR))/NV_MAG(A)命令来提取动床边界每个面单元的平均剪切应力，然后通过式(5)求出修正的作用于床面泥沙颗粒上的剪切应力，进而判定面单元网格上(face)的每个节点(node)是否需要移动。

当判定出网格节点需要移动后，接着需要解决的问题是确定网格节点随着单个时间步长进行移动的位移量，在这里采用泥沙质量守恒方程式(9)来进行计算。

(14)

(15)

求出动床边界每个节点的位移量后，就需要在每个时间步长内对动区域内的网格进行动态更新。由于泄水闸下游冲刷坑局部变形较大，本文动床区域网格更新策略联合采用基于弹簧拉伸变形的光顺法(spring based method)和适用于局部大变形的网格重构法(remeshing method)。

图4　冲刷坑地形测量值和计算值等高线对比Fig.4　Comparison of topographic contours between model test and numerical values

图5　数值计算冲刷坑三维形态Fig.5　Three-dimensional scour topography in numerical simulation

3.5　数值模型计算结果

3.5.1数值冲刷模型验证

为验证本研究所建立冲刷模型的准确性，将计算得到的冲刷坑数据与王新强[27]的物理模型试验数据进行对比分析。

图4(a)为王新强[27]物理模型试验冲刷坑地形等高线图，图4(b)为数值计算得到的冲刷坑等高线云图；图5为数值计算得到的冲刷坑三维形态图。表1为模型试验结果与数值计算结果的比较。由图4、图5和表1可知，粒径分别为300，100，50mm时，冲刷坑深度、长度的数值模拟计算值相对于物理模型试验结果的误差为10%左右。考虑到泥沙输移理论与实际泥沙输移本身存在一定的差异性，可以认为本次冲刷数值模型对于泄水闸下游局部冲刷的模拟具备可靠性。

表1冲刷坑深度、长度对比
Table1Comparisonofscourdepthandlengthbetweennumericalsimulationandmodeltest

组次流量/(m3·s-1)粒径d/mm冲刷坑深度冲刷坑长度试验值/m模拟值/m相对误差/%试验值/m模拟值/m相对误差/%6345．363002．853．201227．5025．38810345．361005．005．10270．4363．301014345．36505．155．35478．8473．507

3.5.2冲刷坑三维形态分析

由于冲刷坑为三维形态，为了更精细地观察，在冲刷坑内进行了截面切片，见图6。

图6　冲刷坑内部截面切片位置示意图Fig．6　Sketchofcross⁃sectionslicepositioninsidethescour

纵向切片平行于x轴，对于粒径为300,100,50 mm的工况分别从动区域中心线位置处沿着y轴负方向每隔3,4,4 m作一个切面；横向切片平行于y轴，对应于以上3个工况分别从冲刷坑最大冲刷深度处沿着x轴正方向每隔3,15,15 m作一个切面。图7为3组不同工况下冲刷坑纵向、横向切片形态分析结果。可以看出：在水流条件一定的情况下，随着泥沙颗粒粒径的减小，冲刷坑的冲刷深度与范围逐渐增大。冲刷坑纵向截面切片图(图7(a)、图7(c)、图7(e))中，动区域中心线位置处冲刷最为剧烈，随着纵向截面切片位置的偏移，由于水流作用强度的减弱，冲刷坑的深度与范围逐渐减小；冲刷坑横向截面切片图(图7(b)、图7(d)、图7(f))中，随着横向截面切片位置的偏移，冲刷坑的深度和范围也逐渐减小。可见，本文模拟出的冲刷坑形态具有较好的三维特性。

图7　各工况下的冲刷坑横向、纵向截面切片形态分布结果Fig.7　Transverse and longitudinal cross-section slices inside the scour under different conditions

4　结　语

基于Fluent动网格技术，采用VOF水-气界面追踪方法来处理闸坝下游冲刷坑内复杂的自由液面，在对Fluent软件进行二次开发的过程中，为避免床面坡度接近泥沙休止角时造成的输沙率数值计算的异常，采用了改进的泥沙起动临界剪切应力计算公式，对当卡水电站泄水闸下游局部冲刷进行了三维数值模拟。在冲刷坑未进行冲刷时，由于泄水闸下游水深较浅，河床床面处得到较大的剪切应力值，该数值大于泥沙起动临界切应力值，进而引起床面泥沙的输移，使得床面高程逐渐下降。随着床面高程的降低，河床床面处的剪切应力值逐渐减小，当该数值等于或小于泥沙起动临界切应力值时，床面泥沙颗粒不再进行运动，河床高程逐渐稳定下来，形成了冲刷坑冲刷后的最终形态。

本文将计算得到的泄水闸下游局部冲刷坑最大深度和形态与物理模型试验值进行了比较，发现两者吻合较好，验证了本文冲刷数值模型的可靠性。

参考文献：

[1]毛昶熙，周名德，柴恭纯.闸坝工程水力学与设计管理[M].北京：水利电力出版社，1995.

[2]严晓达，刘旭东，李贵启，等.低水头引水防沙枢纽[M].北京:水利电力出版社,1990.

[3]高改玉.横丹水电站三维流场模拟及下游二维冲刷坑水力特性研究[D].杨凌：西北农林科技大学,2013.

[4]史志鹏.低水头泄水建筑物消能措施数值模拟[D].杨凌：西北农林科技大学,2010.

再次，你一个小小的意见箱子，有那么金贵么？还要人家安个摄像头保护你？怕谁把你偷回去当枕头枕着，还是当凳子坐着？还是怕谁家缺柴火，偷你回家劈吧劈吧当柴烧了？不怕偷？不怕偷安什么摄像头？还正对着意见箱，明明就是怕人偷么，还狡辩！就冲这，我就对意见箱有意见，摆谱儿摆得也忒大了！

[5]毛荣生，吴飞. 陆水水利枢纽坝下游河床演变规律分析[J].长江科学院院报,1993，10(4):43-49.

[6]刘沛清，冬俊瑞，余常昭. 挑射水流对岩基河床冲刷的探讨[J].长江科学院院报,1994，11(4):1-9.

[7]王晓松,陈璧宏.挑流冲坑三维紊流场的数值模拟[J].水力发电学报,1999，(3):53-61

[8]郭红浩，张根广.九龙峡水电站溢流堰三维数值模拟[J].长江科学院院报,2010，27(3):34-37.

[9]杨建明,吴建华.动网格技术数值模拟挑流冲刷过程[J].水动力学研究与进展，2001， 16(2)：156-161.

[10] 劳尔平.水下淹没建筑物局部流场及河床冲刷坑的数值模拟[D].北京：北京交通大学，2006.

[11] 邓军,许唯临,杨忠超,等.基岩冲刷的数值模拟[J].水科学进展,2005,16(1):47-51.

[12] 李艳富.冲刷坑冲刷过程的数值模拟[J].科学技术与工程, 2009, 9(5): 1176-1180.

[13] 凌建明,林小平,赵鸿铎.圆柱形桥墩附近三维流场及河床局部冲刷分析[J]. 同济大学学报(自然科学版),2007, 35(5): 582-586.

[14] 韦雁机，叶银灿.床面上短圆柱体局部冲刷三维数值模拟[J].水动力学研究与进展(A)，2008，23(6): 655-661.

[15] OLSEN N. Three-Dimensional CFD Modeling of Self-forming Meandering Channel[J]. Journal of Hydraulic Engineering，2003，129(5): 366-372.

[16] USHIJIMA S. Arbitrary Lagrangian-Eulerian Numerical Prediction for Local Scour Caused by Turbulent Flows[J]. Journal of Computational Physics， 1996，125(1): 71-82.

[17] NOBUHISA N, TAKASHI H, TATSUAKI N,etal. Three-dimensional Numerical Model for Flow and Bed Deformation around River Hydraulic Structures[J]. Journal of Hydraulic Engineering, 2005, 131(12): 1074-1087.

[18] 祝志文, 刘震卿. 圆柱形桥墩周围局部冲刷的三维数值模拟[J]. 中国公路学报，2011，24(2): 42-48.

[19] 熊文，汪吉豪，叶见曙.结构形式对桥墩局部冲刷三维性态发展的影响[J].东南大学学报(自然科学版)，2014, 44(1): 155-161.

[20] CHEN Xiao-li, MA Ji-ming, DEY S. Sediment Transport on Arbitrary Slopes: Simplified Model[J]. Journal of Hydraulic Engineering, 2010, 136(5): 311-317.

[21] VAN RIJIN L C. Sediment Transport, Part I: Bed Load Transport[J]. Journal of Hydraulic Engineering, 1984, 110(10): 1431-1456.

[22] 张红武,汪家寅.沙石及模型沙水下休止角试验研究[J].泥沙研究,1989， (3): 90-96.

[23] VAN RIJIN L C. Principles of Sediment Transport in Rivers, Estuaries and Coastal Seas[M]. Amsterdam, the Netherlands: AQUA Publications, 1993.

[24] 钱宁，万兆惠.泥沙运动力学[M].北京：科学出版社，2003.

[25] MELVILLE B. Pier and Abutment Scour-integrated Approach[J]. Journal of Hydraulic Engineering, 1997, 123(2): 125-136.

[26] 王新强，张根广，徐龙飞.当卡水电站消力池的优化试验研究[J].水电能源科学，2014, 32(3): 122-125.

[27] 王新强.当卡水电站下游冲刷坑形态的试验研究[D].杨凌：西北农林科技大学，2014.