桥墩局部冲刷动态模拟及不同截面的冲刷特性

王 飞,姚磊华,张 彬,齐剑峰,李秀芬

(1.河北地质大学勘查技术与工程学院,河北 石家庄 050031; 2.中国地质大学(北京)工程技术学院,北京 100083)

洪水冲刷是桥梁破坏、垮塌的主因之一。桥梁冲刷包括河床的自然演变冲刷、一般冲刷和局部冲刷。其中局部冲刷主要发生在桥墩基础周边,冲刷深度远大于其余两种,所以对于桥梁的破坏最为严重。Lin等[1]对已有的桥梁冲刷破坏事件分析得出,局部冲刷导致的破坏占64%。局部冲刷引起的桥梁承载力的变化对于桥梁的安全稳定及防护意义重大。

桥墩局部冲刷的现场监测难度及不确定性很大,传统的试验研究存在尺度效应影响,对一些点或面的监测存在局限,数值模拟技术作为一种研究冲刷的手段得到了越来越多的关注。Roulund等[2]采用非稳态的RANS模型结合k-ω湍流模型,对圆柱形墩柱周围的流场和冲刷进行模拟,冲刷的发展过程与试验较为吻合,但是平衡冲坑深度低于实际深度15%。Zhao等[3]提出了一种三维有限元方法(FEM)来求解圆柱绕流流场及冲刷,URANS方程通过k-ωSST湍流模型来封闭,计算的冲坑深度小于实际深度10%～20%。韦雁机等[4]基于OpenFOAM开源软件的动网格技术,用输沙率计算床面地形随时间的变化,构建起桩周局部冲刷的动态三维数学模型。祝志文等[5]根据床底泥沙的单宽体积输沙率得到河床高程坐标的瞬时变化,采用边界自适应网格技术修改动边界计算域网格,得到圆柱形桥墩周围局部冲刷坑的演化过程。Khosronejad等[6]对不同横截面形状的桥墩进行了三维动床模拟,采用了流固耦合曲线浸入边界的技术。Ehteram等[7]运用SSIIM软件对桥台的冲刷过程进行了三维模拟,得到了冲刷坑深度和形状并与试验结果进行了比较。

为了克服URANS模型对湍流脉动能量预测的不足,一些学者提出了采用分离涡(DES)和大涡(LES)模拟的方法[8-10]。但计算消耗过大,还很难应用于实际工程,在此基础上进一步考虑泥沙输运的模拟更是具有很大挑战。冲刷达到平衡所需的时间尺度一般为小时或天,而湍流脉动的时间尺度为秒或者更小,如此大的时间悬殊,使得采用DES或LES进行水流-泥沙耦合模拟不切实际。

笔者采用非稳态的RANS模型结合k-ωSST湍流模型,利用CFD软件FLUENT的用户自定义函数功能和动态网格更新技术,实现了桥墩冲刷过程的可视化模拟,并对圆柱形、尖角形和流线形桥墩的冲刷特性进行了分析。

1 动态网格模拟

1.1 计算模型

数值模拟根据Melville试验[11]条件进行建模分析。试验水槽长19 m,宽0.456 m,水深0.15 m,模型桥墩为圆柱形,放置在水槽中心,直径a=0.050 8 m,水槽底部铺设的中值粒径D50=0.385 mm的泥沙,休止角为32°。水流平均来流速度为0.25 m/s。数值模拟的流场尺寸及网格划分如图1所示。网格类型为四面体非结构化网格,网格单元为154 540,河床面和桥墩附近的网格进行局部加密。在床面设置了边界层来提高床面剪应力的计算精度，边界层最小高度0.001 m,增长比率1.2,共设置4层。桥墩附近的网格采用尺寸函数功能进行加密,最小尺寸设为0.002 m,增长比率1.2,最大尺寸为0.02 m。由于流场分布相对于x轴基本对称,所以实际数值计算时以x轴为对称轴取其中一半进行计算以节约计算时间。

图1 流场尺寸及网格划分

床面设置为粗糙壁面边界(粗糙高度2.5D50),桥墩表面为光滑壁面,侧壁、对称面以及顶面均设置为对称边界,水流入口为速度入口边界,出口为自由出流边界。对不同湍流模型的试算表明，k-ωSST模型对湍流涡的模拟最为充分,所以选用该湍流模型。

1.2 临界剪应力的修正

在平床模拟中可以根据Shields曲线得到泥沙的起动剪应力,但是在动床模拟中,随着冲刷的发展,床面逐渐具有了一定的坡度,泥沙沿斜坡方向的重力分力使泥沙更易于起动。随着冲刷的发展,床面坡度会逐渐增大,泥沙在斜坡方向的重力分力逐渐增大,从而加速了泥沙的起动,如果沿用平床试验得到的Shields曲线计算起动剪应力,则将造成较大的误差。通常的做法是引入一个修正系数r对起动剪应力τ0进行修正,因此考虑坡度影响的临界剪应力可以表示为

τcr=rτ0

(1)

图2 床面单元的横向及纵向坡角

关于r的计算方法,Lane公式[12]仅考虑了垂直于流动方向的横向坡角α的影响,平行于流动方向的纵向坡角θ设为0(图2);Ikeda公式[13]引入升阻力比值η,但是仍然没有考虑纵向坡的影响；Dey[14]提出的公式同时考虑了横向坡和纵向坡角的影响,并将升阻力的影响考虑到公式中,其表达式为

(2)

式中：φ为泥沙休止角。

根据式(2)可以画出r分别随横向坡角与纵向坡角的变化曲线(图3)。在坡度相同的情况下,横向坡角α对剪应力的影响小于纵向坡角θ,而已有的冲刷数值模拟中对剪应力的修正往往只考虑了横向坡角的影响,会造成较大的误差。

图3 横向与纵向坡度对修正系数的影响

Bihs等[15]将以上几种方法应用到丁坝冲刷的数值模拟中并与试验结果进行了对比,发现不考虑升力影响的方法计算获得的冲刷深度要小于实际的冲刷深度。冲刷的横向和纵向坡度在冲刷发展中均变得较大,所以考虑两个坡度影响的Dey公式是最为精确的修正系数计算公式，因此选用该公式进行模拟。

在数值模拟中实现变临界剪应力方法,关键点在于正确地确定床面单元的纵向坡角θ和横向坡角α。

图4 床面单元θ和α的确定方法

每个床面单元的法向矢量可以通过FLUENT的宏函数F_AREA()获得,从而可以通过向量运算获得床面单元与xOz平面交线及与yOz平面交线的方向向量(分别设为a和b)。以纵向坡角θ的求解为例(图4),计算公式如下:

(3)

(4)

横向坡角α的求解方法同理。

1.3 推移质输运方程

推移质输运方程总体上都可以用以下无量纲形式表示[16]:

(5)

(6)

式中:μd为动摩擦因数;Ub为推移质泥沙颗粒的平均输运速度;k为常数,取值范围为0.4～0.7。

已有的推移质输运方程多为基于非扰动流(无阻流物)的冲刷试验得到的经验或半理论半经验的公式,实际的桥墩冲刷在桥墩附近会产生湍流涡系结构,增加局部冲刷的深度。而涡的增强效应可以通过修正输运方程中的无量纲沙床剪应力τ*来修正。

由湍流涡引起的向心力可以表示[18]为

(7)

式中:ρs为泥沙的颗粒密度;γ0泥沙颗粒的半径;u为流速;Ω为涡度。

(8)

(9)

作用在泥沙颗粒上的总的无量纲剪应力可以修正为

(10)

1.4 河床动态地貌模型

数值模拟中要得到河床面高程的实时变化,可以通过泥沙的连续性建立河床的动态地貌模型,其统一的形式可以用Exner方程表示:

(11)

式中:p为泥沙的孔隙率;q为泥沙的单宽体积输沙率;Es-Ds为泥沙悬移质下沉通量与上浮通量之差;z为床面高程;t为时间。为了简化计算,模拟时暂不考虑悬移质输沙的影响,所以Es-Ds一项忽略不计。

1.5 泥沙坍塌的调整

桥墩的局部冲刷属于局部的一种较大的变形,所以,随着河床面高程的变化,床面的横向和纵向坡度都会不断地增大。但实际情况是当坡度大于泥沙的自然休止角(泥沙的内摩擦角)后,泥沙会不断地坍塌调整,冲坑的侧壁坡度会小于或等于泥沙的休止角。在数值模拟中,如果只考虑输沙率对床面高程的影响,则会出现冲坑侧壁坡度大于泥沙休止角的情况,这显然与实际不符。

为了考虑坍塌的影响,在数值模拟中,每个时间步网格节点高程根据地貌模型调整后,都需要判断每个单元平面外法向方向与z轴正向的夹角大小(图5),若夹角大于φ-2°,则z向坐标大于单元中心z向坐标的节点(如图5A节点)要下移,而小于的要上移(如图5B、C节点)。为了保证数值稳定,每步调整的移动距离为根据输沙率计算得出的该节点高程变化值,只是移动的方向做了调整。

图5 泥沙坍塌调整方法

1.6 UDF程序及网格更新方法

FLUENT中的动网格模型,可以用于模拟流体域边界随时间改变的问题。动网格中的UDF宏主要有3种,其中DEFINE_GRID_MOTION宏提供了对节点和网格最大限度的操作,因此适用于局部冲刷这样的局部大变形问题。

采用DEFINE_GRID_MOTION宏控制底部边界各个节点运动时,每一个时间步都必须执行。在每个时间步开始计算时,比较床面(动边界)各节点实时剪应力τ与床沙起动临界剪应力τc,若存在超临界剪应力(τ-τc>0),该点表现为冲刷,节点下移;否则表现为静止,节点位置不变。在每个时间步某一节点下移的距离取决于河床动态地貌模型和时间步长的大小。时间步长的选取原则是:在每个时间步内,任一网格运动的最大位移不能超过泥沙颗粒的粒径。

床面的剪应力是随着冲刷的过程而不断变化的,FLUENT中通过提取每一网格单元实时剪应力除以单元的面积得到剪应力。在计算每个单元的变临界剪应力时,有3个角度需要通过编写程序间接得到,即河床的纵向坡角、横向坡角和河床面与床面单元外法向的夹角。这3个角度需要通过提取向量并进行相应的运算得到,具体方法参考1.2节。

推移质输运方程修正时需要计算涡度Ω,而Ω=(Ωx,Ωy,Ωz),这里主要考虑的是Ωy,可以通过FLUENT软件中C_DUDZ函数与C_DWDX函数之差计算出来。

对于冲刷问题,局部运动边界位移大,这里选用光顺法和网格重构法结合的方式来对网格进行重构如图6所示,虽然局部变形较大,但是始终保持了较好的网格质量,确保计算结果的收敛和流场计算的准确度。

2 模拟结果与试验结果[11]对比

2.1 冲坑形态

Melville[11]试验研究表明,冲坑发展试验中前30 min发展较为剧烈,而之后冲刷发展变得缓慢。

30 min时的冲刷深度达到4 cm,为总冲刷深度的75%。为了减少计算消耗,数值模拟取前30 min进行研究。

图7 冲刷坑形态对比

图7给出了数值模拟得到的30 min时冲坑的形态,由于数值模拟以x轴为对称轴取其中一半进行模拟,所以效果图中另一半是由镜像获得。

从图7可以看出:①试验和数值模拟的30 min时的冲坑最大深度均为4 cm,但是试验获得的最大冲刷深度位置从两侧约70°的位置贯通到了桥墩迎水面正前方,而数值模拟得到的最大冲刷深度位于两侧而并没用贯通到正前方,因此在桥墩正前方的最大冲刷深度约比试验结果小25%。造成这一误差的原因可能为:为了节约计算时间,以x轴为对称轴,取一半进行计算,因此x轴附近的网格节点变为边界,受边界效应的影响,桥墩正前方位置可能会造成一定的误差;模型对圆柱形墩正前端湍流脉动的模拟能力所限。②从冲刷等值线以及纵横向截面图对比来看,整体分布形态较为吻合,数值模拟得到的桥墩上游纵向坡度及横向坡度约为30°,这与泥沙的自然休止角非常吻合。桥墩下游坡度较缓,约为20°。试验和模拟的冲坑顶面宽度均为10 cm。

2.2 流场结构

墩周马蹄涡及墩前流线如图8所示,无论是墩前冲坑内的涡系还是延伸到两侧的马蹄涡都与试验情形非常吻合。墩后侧的涡系在现有文献中缺乏形象的描述,所以没有进行对比分析。模拟的尾涡轴线垂直于河床面且涡管直径由下向上不断增大,到接近水面的位置有表面的水平涡,这与很多试验的尾涡描述非常吻合。

图8 墩周马蹄涡及墩前流线

不同高程处冲刷坑内流速矢量图如图9所示,不同高程处的冲坑内的流速矢量分布形态非常相似。在冲坑顶部位置(z=0 cm、z=-1 cm高程处)没有出现明显的马蹄涡系,但在z=-2 cm、-3 cm处出现矢量的明显分离,这表示有了明显的马蹄涡。

以上分析表明,数值模拟结果与试验结果基本吻合,模拟结果可靠。

3 不同截面桥墩的冲刷特性

图9 不同高程处冲刷坑内流速矢量

已有文献对于圆柱形桥墩、方形桥墩以及尖角形的桥墩探讨较多,一致的结论为在相同的阻塞比情况下,冲刷深度由大到小为:方形桥墩、圆柱形桥墩、尖角形桥墩,但是对于流线形的桥墩在已有的冲刷方程和试验研究中较为少见。Arneson等[19]指出,流线形桥墩截面可能会降低桥墩迎水面的流动分离,从而缩减马蹄涡和尾涡的强度,降低冲刷的深度。於刚节[20]提出了一种新型的整体流线形曲面丁坝来改善丁坝附近水流结构,进而减少坝头局部冲刷,通过数值模拟发现在相同流量条件下,最大冲刷深度在曲面丁坝附近比在梯形丁坝附近减少了约20%,且冲刷坑范围也有减小。Al-Shukur等[21]对10种不同形状的桥墩的试验研究得出,矩形截面桥墩的冲刷深度最大,流线形桥墩冲刷深度最小。

利用本文的动床模拟方法对尖角形桥墩和流线形桥墩进行了模拟。水流深度、来流速度等均与Melville[11]试验条件保持一致。尖角形桥墩和流线形桥墩最宽位置的宽度等于圆柱形桥墩直径,长宽比取为3。

图11 不同截面桥墩局部冲坑等值线(单位：cm)

数值模拟得到的冲坑深度随时间的变化曲线如图10所示。可以看出,冲刷速度在开始阶段较快,10 min后变缓,符合试验得到的冲刷逐渐变缓的规律。流线形桥墩在30 min时的冲坑深度比圆柱形桥墩降低约45%,比尖角形桥墩降低约40%。

图10 冲坑深度随时间的变化

根据图11,不同截面的桥墩在冲刷10 min和冲刷30 min时,无论是冲坑的深度还是冲刷的范围,圆柱形最大、尖角形次之、流线形最小。因此,从降低局部冲刷的角度来看,流线形桥墩无疑是一种比尖角形桥墩更理想的桥墩形式。

4 结 论

a. 提出的桥墩局部冲刷过程动态模拟方法对冲刷的发展过程、冲坑的形态以及最大冲刷深度都达到了较为精确的模拟,虽为局部的大变形问题,但是网格质量始终较高。

b. 以下几个方面是动床模拟精度的保证:选用k-ωSST湍流模型,在进行床面剪应力修正的同时考虑纵向和横向坡度的影响,泥沙输运方程考虑涡的增强效应,进行泥沙坍塌的调整及网格的局部加密及重构。

c. 流线形桥墩的冲坑深度比圆柱形桥墩降低约30%,比尖角形桥墩降低约20%,冲坑范围也有所减小,是一种减冲效果较为理想的桥墩截面形式。

参考文献：

[ 1 ] LIN C,HAN J,BENNETT C,et al.Case history analysis of bridge failures due to scour[C]//International Symposium of Climatic Effects on Pavement and Geotechnical Infrastructure, Alaska:2014.

[ 2 ] ROULUND A, MUTLU SUMER B,FREDSØE J,et al.Numerical and experimental investigation of flow and scour around a circular pile[J].Journal of Fluid Mechanics,2005,534(7):351-401.

[ 3 ] ZHAO M,CHENG L,AN H.A finite element model for wave-induced scour below a pipeline[J]. Endocrine Reviews, 2006, 13(2):207-19.

[ 4 ] 韦雁机,叶银灿,吴珂,等.桩周局部冲刷三维数值模拟[J].海洋工程, 2009, 27(4):61-66.(WEI Yanji,YE Yincan,WU Ke,et al.3D numerical modeling of flow and scour around a circular pile[J].The Ocean Engineering,2009,27(4):61-66.(in Chinese))

[ 5 ] 祝志文,刘震卿.圆柱形桥墩周围局部冲刷的三维数值模拟[J].中国公路学报,2011,24(2):42-48.(ZHU Zhiwen, LIU Zhenqing.Three-dimensional numerical simulation for local scour around cylindric bridge pier [J].China Journal of Highway and Transport, 2011, 24(2):42-48.(in Chinese))

[ 6 ] KHOSRONEJAD A,KANG S,SOTIROPOULOS F.Experimental and computational investigation of local scour around bridge piers[J].Advances in Water Resources, 2012, 37(1):73-85.

[ 7 ] EHTERAM M,MEYMAND,A M,Numerical modeling of scour depth at side piers of the bridge[J].Journal of Computational and Applied Mathematics, 2015,280(6):68-79.

[ 8 ] PAIK J,SOTIROPOULOS F,PORTÉ-AGEL F.Detached eddy simulation of flow around two wall-mounted cubes in tandem[J]. International Journal of Heat & Fluid Flow,2009,30(2):286-305.

[ 9 ] ESCAURIAZA C,SOTIROPOULOS F.Initial stages of erosion and bed form development in a turbulent flow around a cylindrical pier[J].Journal of Geophysical Research Earth Surface,2011,116(F3):130-137.

[10] 胡彬,水庆象,王大国,等.特征线算子分裂有限元的圆柱绕流大涡模拟[J].水利水电科技进展, 2017, 37(2):27-32.(HU Bin,SHUI Qingxiang,WANG Daguo,et al.Large eddy simulation of flow around circular cylinder combined with characteristic-based operator-splitting finite element method[J]. Advances in Science and Technology of Water Resources, 2017, 37(2):27-32.(in Chinese))

[11] MELVILLE B W.Local scour at bridge sites[R].Auckland:The University of Auckland,1975.

[12] LANE E W.Design of stable channels[J].Transactions of the American Society of Civil Engineers,1955(1):91-100.

[13] IKEDA S.Closure of “incipient motion of sand particles on side slopes”[J].Journal of Hydraulic Engineering,1983,109:784-786.

[14] Dey S.Threshold of sediment motion on combined transverse and longitudinal sloping beds[J].Journal of Hydraulic Research,2003(4):405-415.

[15] BIHS H,OLSEN NRB.Numerical modeling of abutment scour with the focus on the incipient motion on sloping beds[J]. Journal of Hydraulic Engineering,2012,137(10):1287-1292.

[16] APSLEY D D, STANSBY P K.Bed-load sediment transport on large slopes:model formulation and implementation within a RANS solver[J].Journal of Hydraulic Engineering,2008,134(10):1440-1451.

[17] ENGELUND F,FREDSEE J.A sediment transport model for straight alluvial channels[M].Iwa:Iwa Publishing, 1976.

[18] ABBASNIA A H.Improvements on bed-shear stress formulation for pier scour computation[J].International Journal for Numerical Methods in Fluids,2011,67(3):383-402.

[19] ARNESON L A,ZEVENBERGEN L W.Evaluating scour at bridges[R].Washington,D.C.:Federal Highway Adiministration,U.S.Department of Transportation,2012.

[20] 於刚节.正态曲面丁坝附近三维水流及局部冲刷[D].杭州:浙江大学,2016.

[21] AL-SHUKUR A H K,OBEID Z H.Experimental study of bridge pier shape to minimize local scour[J]. International Journal of Civil Engineering and Technology,2016,7(1):162-171.