考虑悬移质效应的桥墩动床冲刷精细化分析方法

熊文+姚浩+C.S.Cai+叶见曙

摘 要：基于欧拉流固两相流理论并考虑悬移质效应进行桥墩动床冲刷精细化分析.针对动床冲刷与清水冲刷的主要区别，采用不平衡输沙法建立床面变形方程，考虑悬移质与推移质质量交换从而间接影响床面变形.通过对计算流体力学（CFD）软件ANSYS Fluent二次开发，利用床底切应力以及底面网格泥沙浓度分别计算推移质输沙率以及悬移质与推移质的交换通量，由此得到床底网格瞬时变化值，从而实现边界网格的动态更新，以此进行桥墩周围河床动床冲刷过程的动态模拟.通过将数值计算结果与经典理论分析及部分试验结果进行对比，在悬移质分布、局部冲刷深度以及冲刷坑形态等方面充分验证所提出动床冲刷计算方法的合理性与准确性.最后通过参数分析证明了悬移质浓度会明显影响动床冲刷的深度与形态，得出了选择动床冲刷模型并采用流固两相流进行桥墩冲刷精细化分析具有显著必要性的结论.

关键词：桥墩；两相流模型；悬移质；动床冲刷；计算流体力学（CFD）

中图分类号：U442.59 文献标识码：A

文章编号：1674-2974（2016）05-0052-09

Abstract：Based on the Eulerian-Eulerian two-phase flow theory， the live-bed pier scour was carefully simulated with the consideration of suspended load. The riverbed variations were obtained using the non-equilibrium sediment transport model by calculating the mass exchange between the suspended load and traction load. By redeveloping a Computational Fluid Dynamics （CFD） software， i.e.， ANSYS Fluent， the sediment transport rates and exchange flux between the suspended load and traction load were calculated using the shear stress of sediment and sediment concentration. By doing this， the riverbed boundary can be real-timely updated according to the calculated riverbed variations to conduct the live-bed scour simulation. The accuracy and rationality of the proposed simulation was fully verified by comparing with the classic theory and several experimental results from the viewpoints of the suspended load distribution， scour depth， and scour hole profile. The significant influence of the sediment concentration on the scour performance was finally proven by a parametric study. It can be concluded that using the live-bed scour models based on the two-phase flow theory should be very necessary for an accurate simulation.

Key words：piers； two-phase flow model； suspended load； live-bed scour； computational fluid dynamics（CFD）

桥墩冲刷（包括桩基冲刷）病害是当今桥梁结构功能失效、丧失其安全性能的最主要原因之一.以美国为例，从1966～2005年，全美倒塌桥梁中（1 502座）58%的破坏桥梁与桥梁基础结构冲刷病害相关，美国交通部已将桥梁基础结构冲刷看做桥梁结构功能及安全性能失效的最常见原因之一[[1-4[].一般来说，桥墩冲刷深度预测主要基于日常检查以及主观经验判断，准确性不高.尽管针对个别特大型桥梁进行冲刷模型实验，但模型相似比难以确定，且实验中人力物力投入较高.利用相关规范中的冲刷深度计算公式可对桥墩冲刷病害进行快速预判，但其计算假设条件苛刻，参数单一，难以准确模拟冲刷三维整体形态，直接用于实际复杂环境中的桥墩冲刷计算时，其结果的准确性显然难以保证[[3].

为准确地进行桥墩冲刷精细化分析与发展趋势预测，计算流体动力学（CFD）数值方法在冲刷分析中得到越来越多的应用.韦雁机采用k-ε湍流模型并考虑推移质输运模拟了短圆柱体周围局部冲刷过程[[5[]；祝志文等对桥墩局部冲刷进行三维仿真模拟，考虑了床面坡度对临界起动切应力的综合影响[[6[]；叶桢采用两相流VOF法模拟丁坝局部冲刷形态[[7[].但这些研究均仅限于桥墩或丁坝发生清水冲刷，并未考虑冲刷过程中悬移质输运与推移质输运的耦合影响以及由此引起的局部冲刷特有性质，即没有进行动床冲刷分析.事实上，众多桥墩基础多在洪水期间发生淘空从而严重影响桥梁安全；此时由于水流流速较大，水流中悬移质浓度较高，局部冲刷显然已不属于清水冲刷范畴，而必须考虑上游来沙对桥墩冲刷以及周边流场的影响，即属于动床冲刷的范畴[[8-9].因此，对于桥墩动床冲刷的精细化分析以及探讨悬移质对动床冲刷形态的影响显得尤为重要.

本文基于欧拉流固两相流理论并考虑悬移质效应进行桥墩动床冲刷精细化分析.针对动床冲刷与清水冲刷的主要区别，采用不平衡输沙法建立床面变形方程，考虑悬移质与推移质层质量交换从而间接影响床面变形.通过对计算流体力学（CFD）软件ANSYS Fluent进行二次开发，利用床底切应力以及底面网格泥沙浓度分别计算推移质输沙率以及悬移质与推移质的交换通量，由此得到床底网格瞬时变化值，从而实现边界网格动态更新，以此进行桥墩周围床底动床冲刷过程的动态模拟.同时选取合适的物理模型利用以上方法进行数值建模，将数值计算结果与经典理论分析及局部试验数据进行对比，在悬移质分布、局部冲刷深度以及冲刷坑形态等方面充分验证所提出的动床冲刷计算方法的合理性与准确性.最后，通过参数分析得到悬移质浓度显著影响动床冲刷深度及形态的结论，证明了欧拉流固两相流动床冲刷模型取代常规清水冲刷模型的强烈必要性.

1 动床冲刷

当水流平均流速v小于河床泥沙起动平均流速v0时，非墩台周边床面泥沙处于静止状态.但墩台附近水流由于阻水效应流速迅速增大并在局部形成绕流漩涡，漩涡剧烈淘刷桥墩周边泥沙，引起泥沙推移、跳跃等运动，并跟随水流移动，最终带走泥沙形成明显 的局部冲刷坑.随冲刷床面不断下切，坑内绕流流速逐渐降低至坑内泥沙起动流速，此时局部冲刷停止.这一冲刷过程称为清水冲刷.在清水冲刷条件下，主要考虑推移质输运，即床面变形方程仅包括推移质单宽输沙率的影响.

而事实上，一旦水流平均流速v大于床沙起动平均流速v0时，非墩台周边床面泥沙也会随水流运动上浮，部分泥沙运动一段时间后回落至床面，而另外一部分悬浮于水流中.显然局部冲刷坑内将得到上游来沙的补给.一旦单位时间内上游落入冲刷坑内泥沙量与漩涡卷走泥沙量相等，局部冲刷即停止.这一冲刷过程称为动床冲刷.动床冲刷中由于床面泥沙大量起动，水流中悬浮泥沙量很大，悬移质泥沙输运以及其对河床变形的影响将不可忽略.

1.1 欧拉流固两相流模型

按上文所述，冲刷坑形成过程中，水流挟带泥沙共同运动；而泥沙的存在又反过来改变水流的物理性质与紊动结构，从而影响其能量、流速等分布.传统数值计算方法忽略这种流（水流）固（泥沙）相互耦合作用，仅采用单相流进行计算.事实上，动床冲刷条件下，床面泥沙大量起动，相互耦合作用不可忽略，必须采用多相流模型进行模拟.即将河床中泥沙同样视作一种“流体”，认为其与水流共同存在于同一空间（河床水流界面附近）并相互渗透，却又有着各自不同的速度、温度、密度以及不同的体积分数.显然采用多相流模型来模拟桥墩局部三维流场，将更接近实际情况.

常用多相流模型有VOF模型、Mixture模型以及欧拉模型.其中欧拉模型是最为复杂的多相流模型，可以对各相单独进行计算，且每相均有各自独立的守恒方程，力学概念清晰，计算精度较高，一般适用于包括气泡柱、上浮、颗粒悬浮和流化床的模拟.因此，本文采用欧拉流固两相流（泥沙：固相；水流：流相）模型来分析桥墩周边局部冲刷.具体来说，将颗粒（泥沙）作为拟流体，从而认为颗粒（泥沙）与流体（水）是共同存在且相互渗透的连续介质，并同时在欧拉坐标系下进行处理.计算机对这种连续流体模型场的求解主要遵循质量、动量和能量3个守恒原理进行，与此相对应的控制方程包括质量守恒方程、动量守恒方程以及能量守恒方程.因在本问题中水流为不可压缩流体，故不必考虑状态方程作为补充条件.

1.2 冲淤模型

1.2.1 河床变形方程

动床冲刷的河床变形方程基于不平衡输沙法得到，即根据推移质层内不平衡输沙来计算河床变化值.一般认为悬移质在水流中扩散运输对河床变形不产生直接影响，而是通过与推移层之间的质量交换来影响底部河床变形，见式（1）.

1.2.2 推移质输运

单位宽度单位时间内推移质泥沙的输沙体积，为推移质单宽体积输沙率，可按如下经验公式得到[[11[]：

1.3 动网格技术

本文局部冲刷所引起的床底网格变形采用计算流体动力学（CFD）软件Fluent中动网格技术实现，变形通过自定义UDF控制.具体来说，首先计算得到床底切应力以及悬移质浓度，以此得到推移质单宽体积输沙率以及悬移质质量交换.利用推移质单宽输沙率值和悬移质上浮下沉通量之差，并结合离散化床面变形方程计算网格节点瞬时变化量.具体来说，床面变形方程（见式（1））中面单元中心高程变化梯度h/t利用式（12）完成离散化处理.

最终基于式（12）得到的各网格节点高程瞬时变化计算值实时改变网格节点纵坐标，以此实现桥墩周围床底局部冲刷过程的动态模拟.

2 数值仿真模型

2.1 模型描述

对CFD计算软件Fluent进行二次开发，按本文提出的建模理论与方法，建立相应的桥墩冲刷精细化空间数值仿真模型.

基于经典Melville冲刷试验，数值仿真模型区域选为长101.6 cm，宽45.6 cm，深15 cm的水槽（图1），水槽中放置直径为5.08 cm的圆柱作为桥墩，圆柱型桥墩中心距水槽两侧距离为22.8 cm[[12[].床底泥沙平均粒径d50为0.385 cm，水流平均速度v为0.30 m/s，泥沙休止角为32°，泥沙密度为2 680 kg/m3.由于桥墩周围水流流动梯度变化较大，需要对桥墩周围一定范围内网格进行局部加密，相应网格划分后的数值仿真模型见图2，其中最大网格尺寸0.015 m，底层网格厚度0.005 m，总网格数34 793个，计算压力速度耦合采用PISO算法，方程空间离散格式采用quick格式，网格与计算时间步均已满足CFD模拟要求.这种网格加密划分方法可明显减少外侧网格倾斜程度及数量，有助于将倾斜网格线的影响降低；采用数学方法检测多种划分网格的畸角以及最大、最小网格尺寸等指标，并进行对比，本文所采用的划分网格质量远优于其他划分方法.

另外，本研究判别清水冲刷、动床冲刷的关键参数床沙起动平均流速v0参考张瑞瑾公式确定[[13[].

2.2 边界条件

基于欧拉流固两相流理论的动床冲刷模型，其边界条件设置的关键在于入口处的水流流速分布、悬移质速度分布以及悬移质浓度.

采用与图1相同的数值仿真模型，首先假设入口水流速度分布为断面平均流速0.30 m/s，不加载动态网格更新程序，计算获得出口处水流速度分布，然后将此时该出口速度分布作为实际冲刷计算时的模型入口水流速度分布，如图3所示，其中横坐标z表示水深.这样做可以使得入口流速分布更加符合实际情况.

水流中悬移质泥沙分别具有水平和垂向2个运动速度.其中水平运动速度与水流运动速度一致，即在模型入口处两者设置为相同.对于垂向运动速度，可近似认为入口处紊动水流的扩散作用尚未完全发生，入口处泥沙仅仅具有重力下沉的沉速，该速度采用式（10）计算得到，即s=0.054 286 m/s.显然这样设置是符合实际情况的，另外经过扩散的泥沙计算分布将在下节中作为流场模拟准确性的验证对象.

对于悬移质在入口处的浓度分布，本文参考K. Debnath试验中悬移质浓度范围为993～1 332 g/m3[[14[]，按泥沙密度2 680 kg/m3换算为泥沙体积分数为0.000 370～0.000 497.据此，本研究选取入口处泥沙体积分数为0.000 4.

参考桥墩冲刷流场普遍模拟方法，模型中其他边界的设置如图4所示，其中流场与桥墩以及河床的界面分别设置为光滑与粗糙壁面边界[[15].

3 准确性验证

通过将上文提出的基于欧拉流固两相流理论的动床冲刷模型数值计算结果与经典理论分析及局部试验数据进行对比，在悬移质分布、局部冲刷深度以及冲刷坑形态等方面进行该模拟计算方法的准确性验证.

3.1 悬移质分布验证

悬移质模拟方法对局部冲刷数值模拟的准确性起着至关重要的作用.基于钱宁、万兆惠针对悬移质垂向分布的扩散理论[[11[]，对本文提出的悬移质模拟方法进行准确性验证.

悬移质分布扩散理论认为，当悬移质含沙量的垂线分布达到平衡状态时，泥沙扩散成为一个均匀、恒定的问题，据此建立微分方程为：

根据上述扩散理论以及模型数值计算中某一已知点体积比含沙量Sva，便可分别得到墩前、墩后以及下游靠近出口处的A（x=0.25 m），B（x=0.35 m）和C（x=1.0 m）3点（具体位置如图5所示）悬移质体积分数垂向分布理论曲线，并将其与本文所提出的动床冲刷数值模型计算结果进行比较，如图6所示.

从图6（a）～（c）可看出，墩前位置A，墩后位置B以及出口处C的悬移质体积分数数值模拟计算结果与扩散理论分析结果，无论在数值上还是分布规律上均非常相似.深度较小时，体积分数较大，随着深度增大泥沙扩散，体积分数迅速降低并趋于一稳

定数值.仅在A和B位置低高度处略有差别，其主要原因在于A处离水流入口处较近，悬移质分布仍未达到完全平衡分布状态；而B处由于存在桥墩阻挡，悬移质分布相对受到影响，使其与理论计算值出现局部偏差.显然对于水流及悬移质均充分发展的位置C，数值模拟计算结果与扩散理论分析结果已

极为相近，仅在个别位置存在微小差别.因此，可以认为本文提出的基于欧拉流固两相流理论的动床冲刷模型能够准确模拟流场中悬移质分布.

3.2 局部冲刷深度验证

本研究采用悬移质体积分数为0.000 4，对应质量浓度为1 072 mg/L，数值模拟得到30 min内局部冲刷深度与时间的关系如图7所示.

从图7可看出，计算结果具有较强的收敛性和稳定性（初始冲刷较快，后期逐渐放缓，并最终达到平衡冲刷深度），30 min时最大局部冲刷深度为6.485 cm.而在同样参数条件下清水冲刷计算得到的最大局部冲刷深度仅为5.855 cm，两者冲刷深度比值为1.11倍.该比值与既有试验数据相比，虽没有达到K. Debnath等提出的1.54倍那么大[[14[]，但鉴于K. Debnath的试验是基于黏性土床面，而本文基于沙性土床面，因此动床冲刷的计算准确性仍可得到定性验证.

另外可以发现，清水冲刷深度计算结果要明显小于代表真实情况的动床局部冲刷计算深度.其主要原因在于动床模型考虑悬移质影响，而悬移质运动会导致水流冲刷河床能量的增强.定性地说，1）与清水冲刷相比，动床冲刷中水流有一部分空间被悬移质占据，减少了通过黏性转化为热能的流体有效空间，使得水流能量损失减少，同时由于悬移质质量大于水流，导致总体动能增大，冲刷河床能力增强；2）悬移质存在使挟沙水流黏性增加，导致近壁区层流层厚度增大，占据水力周界的水力粗糙区域，使得周界变得更为光滑，从而减少水流流速的降低或能量的损失，使得冲刷河床能力增强[[11[].也就是说，悬移质的存在会使水流系统能量增加，局部冲刷能力增强，最终导致动床冲刷深度增大.显然以上两点解释均与水流与泥沙之间的相互耦合作用紧密相关，必须通过两相流计算模型才能考虑到在连续流体模型场中泥沙水流之间的相互耦合行为与现象.所以，必须建立两相流动床冲刷模型才能确保冲刷深度/形态计算结果与实际情况的一致性.

3.3 局部冲刷坑形态验证

动床冲刷与清水冲刷的根本区别在于桥墩周围泥沙被绕流漩涡带走而形成的冲刷坑内有无上游来沙的补给.事实上，冲刷过程中局部冲刷坑内泥沙的损失和补给是同时发生的，直至达到平衡状态，因此很难直观地从计算结果中反映冲刷坑内泥沙补给状况.但可以肯定的是，上游来沙沉积于局部冲刷坑内主要源于冲刷坑的不断下切，进而流场空间增大，水流强度与流速下降，导致水流挟沙能力减弱，最终悬移质逐渐沉降.也就是说，在桥墩周围水流流动强度和速度较低的区域将出现悬移质的堆积.

本节将基于悬移质堆积这一形态特点对本文所提出的动床冲刷模型计算准确性进行验证.图8与图9分别给出利用本文动床冲刷数值仿真模型计算得到的流场床面桥墩剖面图以及局部冲刷三维形态图（计算时刻：90 s），可以明显看出桥墩前后均存在高于原床面高程的凸起区域.具体来说，桥墩前方悬移质堆积是由于水流遇桥墩阻挡使局部流速降低甚至流向反向，从而水流挟沙能力减弱导致桥墩前方泥沙堆积.而水流绕过桥墩同样会引起桥墩后方局部流速降低，挟沙能力减弱，进而导致悬移质与少量推移质在桥墩后方沉积.显然，本文提出的基于欧拉流固两相流理论的动床冲刷模型能够准确模拟冲刷过程中悬移质堆积现象，也证明了本文悬移质模拟方法的正确性.

综上，无论从悬移质分布、局部冲刷深度或者局部冲刷坑形态方面与经典理论分析及局部试验数据进行对比，本文所提出的基于欧拉流固两相流理论的动床冲刷模型均能够准确模拟出桥墩冲刷及相应流场特征，计算结果具有较高的可靠性，与常规计算方法相比更加符合桥墩冲刷实际情况.

4 悬移质浓度影响分析

正如前文所述，悬移质是动床冲刷区别于清水冲刷的最重要因素，而悬移质浓度又是影响动床冲刷形态的最关键物理参数，有必要对这种影响进行更加深入的研究.

仍然采用上文冲刷模型，参考K. Debnath等试验中悬移质浓度范围为993～1 332 mg/L[[11[]，按容重2 680 kg/m3换算为泥沙体积分数范围0.000 370～0.000 497.进而计算入口处悬移质各平均浓度下30 min局部冲刷最大深度，并针对入口处悬移质平均浓度（体积分数/体积含沙量）c（10-4）与最大冲刷深度h （cm）的关系进行公式拟合，如图10所示.

从图10可以看出，随悬移质浓度不断增加，桥墩最大冲刷深度也随即增长.这一规律也与洪水期间冲刷明显发展较快的事实相一致，即洪水期间水流流速与挟沙能力增大导致悬移质浓度迅速上升，是洪水期间桥墩容易发生淘空的重要原因.所以，可增铺较大粒径床沙或设置隔离物使水流悬移质浓度降低，从而减小冲刷.值得注意的是，最大冲刷深度在理论上并不一定会随悬移质浓度增加而无限增大，这是因为当其浓度增加到一定程度（一般远超于常规河流可能达到的浓度），水流与悬移质的混合流体已再不属于牛顿流体范畴，其基本性质与特征已无法采用本文中数学表达式描述，由于该部分已超出本文研究范畴，此处不再赘述.需要说明的是，图10中所涉及的悬移质浓度已包含大部分可能出现的泥沙浓度值，结论已具备一定的普遍性.

另外，图11还给出悬移质体积分数为0.000 370， 0.000 420以及0.000 497时最大局部冲刷深度随时间的变化曲线.从图中明显看出，冲刷初期，悬移质浓度越高其冲刷深度越小；而在冲刷后期，悬移质浓度越高反而其冲刷深度越大.这主要是由于在两相流模型中冲刷初期悬移质泥沙的存在会减弱水流的紊动形态，悬移质浓度越高，水流紊动形态越不明显，此时局部冲刷程度也就越弱.而在冲刷发展后期，如3.2节所述，悬移质的存在减少水流能量损失逐渐成为主要因素，且该损失减少会随时间推移更为明显，因此此时悬移质浓度越高，局部冲刷深度也越大.

综上可见，悬移质浓度不仅与最大冲刷深度明显相关，而且同时影响整个冲刷过程的形态发展；也就是说，悬移质在整个冲刷分析中不仅不可忽略，而且不同悬移质浓度甚至会改变冲刷进程特征.所以，利用本文提出的基于欧拉流固两相流理论的动床冲刷模型取代传统的清水冲刷模型势在必行.

5 结论与展望

1）针对动床冲刷与清水冲刷的主要区别，采用不平衡输沙法建立床面变形方程，得以考虑悬移质与推移质层质量交换对床面变形的影响.

2）利用床底切应力以及底面网格泥沙浓度分别计算推移质输沙率以及悬移质与推移质的交换通量，由此得到床底网格瞬时变化值，从而实现桥墩周围床底动床冲刷过程的动态模拟.

3）将数值计算结果与经典理论分析及局部试验数据进行对比，在悬移质分布、局部冲刷深度以及冲刷坑形态等方面充分验证了所提出动床冲刷计算方法的合理性与准确性.

4）通过参数分析得到悬移质浓度显著影响动床冲刷深度及形态，证明了动床冲刷模型并采用欧拉流固两相流分析的必要性，使得计算结果与实际情况更加一致.

5）未来拟利用本文精细化动床冲刷模型对现有规范所采用的计算公式进行准确性修正，同时将在后续研究中积累实验数据，最终完成更为细致的定量对比与验证分析.

参考文献

[1] DENG L， CAI C S. Applications of fiber optic sensors in civil engineering[J]. Structural Engineering & Mechanics， 2007， 25（5）： 577-596.

[2] HUNT E B. Monitoring scour critical bridges[R]. Washington， DC： Transportation Research Board， National Research Council， 2009.

[3] XIONG W， CAI C S， KONG X. Instrumentation design for bridge scour monitoring using fiber bragg grating sensors[J]. Applied Optics， 2012， 51（5）： 547-557.

[4] KATTELL J， ERIKSSON M. Bridge scour evaluation： screening， analysis， and countermeasures[R]. Washington， DC： US Forest Service， 1998：1-10.

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

[6] 祝志文， 刘震卿. 圆柱形桥墩周围局部冲刷的三维数值模拟[J].中国公路学报， 2011， 24（2）： 43-48.

[7] 叶帧. 丁坝附近紊流的三维数值模拟及局部冲刷[D]. 杭州：浙江大学建筑工程学院， 2007： 6-7.

[8] 刘春嵘， 邓丽颖， 呼和敖德. 复杂流动下泥沙起动概率的图像测量[J]. 湖南大学学报：自然科学版， 2008， 35（3）： 24-27.

[9] 周健， 杜强， 于仕才. 泥石流启动试验的数值模拟研究[J]. 湖南大学学报：自然科学版， 2015， 42（9）： 96-103.

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

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

[12]MELVILLE B M. Local scour at bridge sites[D]. Auckland， New Zealand： School of Engineering， University of Auckland， 1975：86-102.

[13]高冬光， 王亚玲. 桥涵水文[M]. 4版. 北京：人民交通出版社， 2008： 150-152.

[14]DEBNATH K， CHAUDHURI S. Effect of suspended sediment concentration on local scour around cylinder for clay-sand mixed sediment beds[J]. Engineering Geology， 2011， 117（3/4）： 236-245.

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