泥沙局部冲淤二维数值模拟仿真

张芝永，拾兵

(中国海洋大学工程学院，山东青岛266100)

泥沙冲淤是河流、海洋中比较普遍的一种自然现象，在近床面位置存在结构物情况下，床面泥沙运动剧烈，局部地形变化强烈，局部冲淤问题严重.对于此类问题，国内外许多学者主要是通过模型实验手段来进行研究分析，其研究对象主要包括桥墩［1］、丁坝［2-3］、海底管道［4-5］等结构物周围的局部冲淤问题.

随着计算机技术的发展，局部冲淤的数值模拟方法越来越引起研究者的重视.Ouillon等［6］采用标准的三维k-ε紊流模型模拟了丁坝周围的流场，并对平整床面时的床面切应力沿河床的分布与最终冲刷坑的床面形态进行了定性比较.彭静等［7］将线性与非线性三维k-ε模型应用于丁坝绕流模拟，并比较了各种模型的模拟效果.Brφrs［8］提出了新的二维冲淤模型采用标准的k-ε紊流模型和同时考虑悬移质和推移质输运，对海底管道局部冲刷过程进行了模拟.Liang和Cheng［9-10］基于有限差分法发展了一个可以精确地预测海底管道冲刷时间历程的二维模型.Liu［11］利用VOF和动网格技术对水平射流情况下的底床冲淤进行了数值模拟，该模型考虑了水面的变化情况，冲坑变化过程模拟值与试验值吻合较好.韦燕机［12］利用OpenFOAM和动网格技术对桩周局部冲刷进行了三维数值模拟.Zhi-wen Zhu［13］利用网格技术对桥墩冲刷进行了数值模拟，验证了网格较好的地形边界拟合效果.综上所述，以上大部分模型均是通过动网格技术来模拟泥沙的局部冲刷，这说明动网格技术要比其他如两相流技术来模拟泥沙更为准确.但上述模型也存在一些不足，如大部分局部冲淤模拟均是通过非通用计算程序来完成的，这限制了其数值模拟方法的推广.

FLUENT是一款目前较为流行的商用CFD软件，它在航空航天，能源、水利、环境等领域得到了广泛的应用，是目前全球功能最强大的计算流体力学商用软件.该软件在水动力计算模块有着广泛的通用性，但该软件并没有利用动网格来计算泥沙冲淤的模块.有鉴于此，作者充分利用FLUENT的二次开发接口，即自定义函数功能，通过C++编写了泥沙输运模块，并将其嵌入到FLUENT中，实现了泥沙局部冲淤的二维数值模拟.

1 水动力模型

1.1 控制方程

水动力模块中，其控制方程为连续性方程和动量方程:

式中:ui为流体速度在i方向上的分量;ρ是流体的密度;μ是流体粘度;p是作用在流体微元上的压力;v是流体的运动粘滞系数，Sij是平均应变速率张量.ui'是i方向速度脉动值为雷诺应力张量，vT紊流运动粘滞系数，δij是克罗内克尔符号.

1.2 湍流模型

对于不可压缩流体且不考虑用户自定义的源项时，标准k-ε模型方程为k方程:

ε方程:

式中:Gk是由于平均速度梯度引起的湍动能k的产生项;C1ε和C2ε为经验常数，取C1ε=1.44，C2ε=1.92;σk和σε分别是与湍动能k和耗散率ε对应的 Prandtl数，取 σk=1.0，σε=1.3.对于床面近底层，采用壁面函数法将壁面上的物理量与湍流核心区求得物理量联系起来.

2 泥沙输运模型

2.1 推移质输沙模型

床面泥沙的起动可以通过临界希尔兹数来确定，其希尔兹数的表达式为

式中:τ为床面切应力，ρ为水的密度，g为重力加速度，s为泥沙比重，d50为泥沙中值粒径.

水平床面上泥沙临界起动希尔兹数为

其中，D*是无量纲泥沙颗粒尺寸，定义为

在有坡度的床面上泥沙临界起动希尔兹数 按式(10)进行修正［14］:

式中:α是砂床面和水平面的坡角，规定上坡为正角，下坡为负角;φ是泥沙的休止角.

平面单宽体积输沙率q0由式(11)求得

推移质单宽体积输沙率可用下式表示:

式中:q0为平面单宽体积输沙率，τ为床面剪应力，τx为床面剪应力x方向分量，h为床面高程，经验系数C取值范围为 1.5 ～2.3［8］.

2.2 床面变形模型

根据输沙平衡，床面高程h的变化可以用冲淤方程表示为

式中:p0为床沙孔隙率，qbx为x方向推移质单宽体积输沙率.

3 动边界处理方法

3.1 动边界处理方法

动网格模型可以用来模拟流场形状由于边界运动而随时间改变的问题.其广泛应用于活塞、阀门及柔性体等有运动边界工况的模拟.床面的冲淤变形形式类似于柔性体变形，需要对各个节点的运动情况进行描述.因此本文采用FLUENT的DEFINE_GRID_MOTION宏命令来给出边界节点的运动方式.在床面节点位置更新后，需要对区域内部网格进行调整.FLUENT提供了3种网格更新方式:1)光顺网格法，2)动态层网格法，3)局部网格重画法.非结构化网格比结构化网格更加适用于不规则边界的变形问题，因此本模型采用非结构化网格，其网格更新方法为光顺网格法和局部网格重画法相结合的方式，这2种方法的网格更新效果可见图1.图1为下文算例中水流下坡面处局部网格变化情况.从图中可以看出，在底部边界发生变化后，区域内网格进行重新划分，各区域网格尺寸尺度基本保持不变，靠近下部变形边界的网格仍然比较密集，这对于确保近壁面流场计算的准确性至关重要.

图1 动网格更新效果Fig.1 Mesh deformation

3.2 沙滑模型

由于泥沙传输的复杂性和与流场相互作用的非线性，在海床高程更新过程中会发生反射，导致网格畸变，造成不存在的冲淤形态和数值的不稳定.为解决此问题，采用Liang提出的沙滑模型［10］:当局部冲淤斜坡角度超过泥沙休止角时，对相应的海床节点进行调节，使冲淤斜坡角等于休止角.

4 数值解法及计算流程

水动力控制方程及紊流方程使用基于单元中心的有限体积法离散，利用非稳态求解器进行求解，瞬态项采用二阶隐格式，对流项采用二阶迎风差分格式，扩散项采用中心差分格式;压力与速度的耦合使用SIMPLE算法.

床面变形方程式(13)通过有限差分法进行离散，其离散后形式为

式中:p代表床面节点序号;n为当前时间步;n+1为下一时间步;Δtb为床面更新时间步长.

为节省计算时间，流场计算的时间步长和床面更新的时间步长采用不同值，床面更新时间步长Δtb要远大于流场计算时间步长Δtflow.其数值模拟过程可概括为以下几个过程:

1)计算多个步长的速度、压力、湍流数及切应力等;

2)利用最后一步得到的床面切应力数据，调用泥沙输运模块，计算推移质输沙率和床面高程变化情况;

3)采用动网格技术，改变床面节点位置，同时重新调整计算区域内部网格;

4)返回步骤1)，重新上面的计算，直至到达指定时间.

5 模型验证

为验证本文所建模型的可靠性，通过2个算例对其进行了验证.

5.1 水渠冲淤

Van Rijn［15］对有坡度水渠的水流情况进行了试验研究，李昌良［16］在其基础上又对其冲淤情况进行了进一步的研究.其试验布置方案如图2所示.左边进口平均流速0.51 m/s.水深0.39 m，沙床泥沙平均粒径d50=0.16 mm.在数值模拟中，其计算区域与试验布置相同，左边进口设置为速度入口，右边为自由出流边界，水面边界设置为对称边界.底面边界设置为wall边界，其粗糙度为2.5d50.整个区域采用非结构化网格进行离散.

图2 试验布置Fig.2 Layout of experiment

图3列出了在初始地形情况下4个垂直断面的流速分布情况，从图中可以看出数值模拟结果与试验结果是较为一致的，这说明了数值模拟中水面采用刚盖假设的合理性，同时为进一步的泥沙冲淤过程研究奠定了基础.

图4为不同时段的网格情况及床面地形变化情况，可以发现在床面地形发生变化后，采用光顺网格法和局部网格重画法更新后的非结构化网格很好地拟合了床面地形的变化，在变化过程中，网格的疏密程度也得到了较好的保持.图5给出了3 h后床面冲淤情况.从图中可以看出，在沟内由于水深的增加，水流流速减小，进而导致床面切应力减小，泥沙在此淤积.在水流下坡区域，床面切应力减小，泥沙淤积较为严重，而在水流上坡区域，下游流速加大，床面切应力增大，从而导致该区域冲刷严重.数值结果与前人结果比较一致，准确地反映了沟内床面的变化规律.

图3 各断面流速分布情况Fig.3 Velocity distributions at different crossvsections

图4 不同时刻网格更新情况Fig.4 Mesh deformation in different time

图5 3 h后床面地形Fig.5 Bed profile at t=3.0 h

5.2 海底管道局部冲刷

本算例对近床面的海底管道局部冲刷的冲刷历程进行详细的分析.

Liang在文献［17］中对有间隙海底管道的局部冲刷进行了详细的数值模拟研究.在这里采用其中一计算算例进行了数值计算，算例中管道直径D为10 cm，管道与床面之间垂直距离G为0.5D，来流流速为0.5 m/s，泥沙中值粒径为 0.36 mm，水深 3.5D.选择距管道中心上下游各20D距离的区域为计算区域.左边界为速度进口边界，右边界为自由出流边界，上部边界采用对称边界，下部边界为wall边界(如图6所示).整个区域采用非结构化网格进行离散，在近壁面处和管道周围进行加密.

图6 计算区域Fig.6 Computational domain

图7 不同时刻x方向流速等值线(单位:m/s)Fig.7 Contours of x veolicity in different time(unit:m/s)

由于动网格计算耗时较长，本文只模拟了80 min内的床面冲刷过程.图7列出了冲刷过程中，不同时刻的x方向流速分布图，该图同样直观地描述了床面的冲刷变形情况.在t=5 min时，由于初始冲刷阶段间隙内流速较大，管道局部冲刷程度较为严重，管道正下方出现沙坑，下游区域出现沙丘.在t=15 min时，上游冲坑深度加大，下游区沙丘逐渐向下游移动并并且沙丘高度逐渐减小.在t=80 min时，下游沙丘消失，整个过程与实际的冲刷过程规律是较为一致的.从图中还可以发现，管道下游并未出现涡旋脱落现象，这可能是由网格尺度的限制以及启动动网格模型而导致的计算精度下降等原因造成的.而在实际试验中，管后区域是存在涡旋脱落现象的.

Liang［17］的研究结果表明，有涡旋脱落时，管道后尾流区范围比较大而且变化比较剧烈，大范围的尾流区使得更多水流通过管道下方孔道流向下游，进而导致近床面流速加大，管下和管后床面切应力也相应增大，同时，涡旋的不断脱落又使得管道下方及后方床面的切应力波动变化，波动幅值有可能会是平均值的数倍.因此，有涡旋释脱落时的冲刷深度要大于无涡旋脱落情况.图8列出了t=80 min时，Liang的模拟结果与本文模拟结果的对比情况.从图中可以看出，管道上游的床面冲刷模拟结果较为接近，管道下游区域床面冲刷模拟结果与Liang的无涡旋脱落时的冲刷结果比较接近，而与有涡旋脱落时的模拟结果则有较大偏差.其原因在于本文模拟过程中并未出现涡旋脱落现象，从而造成管道下方及后方床面切应力计算值偏小，冲刷程度较小.

图8 t=80 min时床面地形Fig.8 Bed profile at t=80 min

6 结论

通过对商用CFD软件FLUENT的二次开发，建立了局部冲淤二维数值模型.该模型的泥沙输运模型通过FLUENT的DEFINE_GRID_MOTION宏命令嵌入，床面边界的冲淤变化利用动网格技术来实现.

数值计算的结果表明，泥沙冲淤模型可以比较准确的预测不同条件下泥沙冲淤过程，可用于模拟以推移质输沙为主的二维泥沙冲淤问题，具有广阔的应用前景.但该模型同样存在一些问题，如对于有涡旋脱落时的泥沙冲淤模拟偏差较大以及计算时间过长等问题.这都有待于进一步优化和研究.

［1］MELVILLE B W.Pier and abutment scour:integrated Approach［J］.Journal of Hydraulic Engineering，1997，123(2):125-136.

［2］彭静，河原能久.丁坝群近体流动结构的可视化实验研究［J］.水利学报，2000，31(3):44-47.PENG Jing，YOSHIHISA Kawahara.Visualization of flow structure around submerged spur dikes［J］.Journal of Hydraulic Engineering，2000，31(3):44-47.

［3］ROGER A，KUHNLE C V，ALONSO F，et al.Local scour associated with angled spur dikes［J］.Journal of Hydraulic Engineering，2002，128(12):1087-1093.

［4］SUMER B M，JENSEN H R，FREDSØE J.Effect of leewake on scour below pipelines in current［J］.J Waterway，Port，Coastal and Ocean Engineering，1988，114(5):599-614.

［5］CHIEW Y M.Prediction of maximum scour depth at submarine pipelines［J］.Journal of Hydraulic Engineering，1991，117(4):452-466.

［6］OUILLON S，DARTUS D.Three-dimensional computation of flow around groyne［J］.Journal of Hydraulic Engineering，1997，123(11):962-970.

［7］彭静，河源能久.线性与非线性紊流模型及其在丁坝绕流中的应用［J］.水动力学研究与进展:A辑.2003，18(5):589-594.PENG Jing ，KAWAHARA Yoshihisa.Application of linear and non-linear turbulent models in spur dike flow［J］.Chinese Journal of Hydrodynamics，2003，18(5):589-594.

［8］BRØRS B.Numerical modeling of flow and scour at pipelines［J］.Journal of Hydraulic Engineering，1999，125(5):511-523.

［9］LIANG D F，CHENG L，LI F.Numerical modelling of scour below a pipeline in currents.part I:flow simulation［J］.Coastal Engineering，2004，52(1):25-42.

［10］LIANG D F，CHENG L，LI F.Numerical modelling of scour below a pipeline in currents.partⅡ:Scour simulation［J］.Coastal Engineering，2004，52(1):43-62.

［11］LIU Xiaofeng，GARCÍA M H.Three-dimensional numerical model with free water surface and mesh deformation for local sediment scour［J］.J Waterway，Port，Coastal and Ocean Engineering，2008，134(4):203-217.

［12］韦雁机，叶银灿.床面上短圆柱体局部冲刷三维数值模拟［J］.水动力学研究与进展:A辑，2008，23(6):655-661.WEI Yanji，YE Yincan.3D numerical modeling of flow and scour around short cylinder［J］.Chinese Journal of Hydrodynamics，2008，23(6):655-661.

［13］ZHU Zhiwen，LIU Zhenqing.CFD prediction of local scour hole around bridge piers［J］.Journal of Central South University of Technology，2012，19(1):273-281.

［14］ALLEN J R L.Simple models for the shape and symmetry of tidal sand waves:statically stable equilibrium forms［J］.Marine Geology，1982，48(12):31-49.

［15］VAN RIJN L C.Principles of sediment transport in rivers，estuaries and coastal seas［M］.Amsterdam:Aqua Publications，1993:1245-1247.

［16］李昌良.泥沙运动与底床变形的数值模拟［D］.青岛:中国海洋大学，2008:58-61.LI Changliang.The numerical simulation of sediment transport and bed evolution［D］.Qingdao:Ocean University of China，2008:58-61.

［17］LIANG Dongfang，CHENG Liang，YEOW K.Numerical study of the Reynolds-number dependence of two-dimensional scour beneath offshore pipelines in steady currents［J］.Ocean Engineering，2005，32(4):1590-1607.