橋墩局部沖刷動態(tài)模擬及不同截面的沖刷特性

王 飛,姚磊華,張 彬,齊劍峰,李秀芬

(1.河北地質(zhì)大學勘查技術(shù)與工程學院,河北 石家莊 050031; 2.中國地質(zhì)大學(北京)工程技術(shù)學院,北京 100083)

洪水沖刷是橋梁破壞、垮塌的主因之一。橋梁沖刷包括河床的自然演變沖刷、一般沖刷和局部沖刷。其中局部沖刷主要發(fā)生在橋墩基礎周邊,沖刷深度遠大于其余兩種,所以對于橋梁的破壞最為嚴重。Lin等[1]對已有的橋梁沖刷破壞事件分析得出,局部沖刷導致的破壞占64%。局部沖刷引起的橋梁承載力的變化對于橋梁的安全穩(wěn)定及防護意義重大。

橋墩局部沖刷的現(xiàn)場監(jiān)測難度及不確定性很大,傳統(tǒng)的試驗研究存在尺度效應影響,對一些點或面的監(jiān)測存在局限,數(shù)值模擬技術(shù)作為一種研究沖刷的手段得到了越來越多的關(guān)注。Roulund等[2]采用非穩(wěn)態(tài)的RANS模型結(jié)合k-ω湍流模型,對圓柱形墩柱周圍的流場和沖刷進行模擬,沖刷的發(fā)展過程與試驗較為吻合,但是平衡沖坑深度低于實際深度15%。Zhao等[3]提出了一種三維有限元方法(FEM)來求解圓柱繞流流場及沖刷,URANS方程通過k-ωSST湍流模型來封閉,計算的沖坑深度小于實際深度10%～20%。韋雁機等[4]基于OpenFOAM開源軟件的動網(wǎng)格技術(shù),用輸沙率計算床面地形隨時間的變化,構(gòu)建起樁周局部沖刷的動態(tài)三維數(shù)學模型。祝志文等[5]根據(jù)床底泥沙的單寬體積輸沙率得到河床高程坐標的瞬時變化,采用邊界自適應網(wǎng)格技術(shù)修改動邊界計算域網(wǎng)格,得到圓柱形橋墩周圍局部沖刷坑的演化過程。Khosronejad等[6]對不同橫截面形狀的橋墩進行了三維動床模擬,采用了流固耦合曲線浸入邊界的技術(shù)。Ehteram等[7]運用SSIIM軟件對橋臺的沖刷過程進行了三維模擬,得到了沖刷坑深度和形狀并與試驗結(jié)果進行了比較。

為了克服URANS模型對湍流脈動能量預測的不足,一些學者提出了采用分離渦(DES)和大渦(LES)模擬的方法[8-10]。但計算消耗過大,還很難應用于實際工程,在此基礎上進一步考慮泥沙輸運的模擬更是具有很大挑戰(zhàn)。沖刷達到平衡所需的時間尺度一般為小時或天,而湍流脈動的時間尺度為秒或者更小,如此大的時間懸殊,使得采用DES或LES進行水流-泥沙耦合模擬不切實際。

筆者采用非穩(wěn)態(tài)的RANS模型結(jié)合k-ωSST湍流模型,利用CFD軟件FLUENT的用戶自定義函數(shù)功能和動態(tài)網(wǎng)格更新技術(shù),實現(xiàn)了橋墩沖刷過程的可視化模擬,并對圓柱形、尖角形和流線形橋墩的沖刷特性進行了分析。

1 動態(tài)網(wǎng)格模擬

1.1 計算模型

數(shù)值模擬根據(jù)Melville試驗[11]條件進行建模分析。試驗水槽長19 m,寬0.456 m,水深0.15 m,模型橋墩為圓柱形,放置在水槽中心,直徑a=0.050 8 m,水槽底部鋪設的中值粒徑D50=0.385 mm的泥沙,休止角為32°。水流平均來流速度為0.25 m/s。數(shù)值模擬的流場尺寸及網(wǎng)格劃分如圖1所示。網(wǎng)格類型為四面體非結(jié)構(gòu)化網(wǎng)格,網(wǎng)格單元為154 540,河床面和橋墩附近的網(wǎng)格進行局部加密。在床面設置了邊界層來提高床面剪應力的計算精度，邊界層最小高度0.001 m,增長比率1.2,共設置4層。橋墩附近的網(wǎng)格采用尺寸函數(shù)功能進行加密,最小尺寸設為0.002 m,增長比率1.2,最大尺寸為0.02 m。由于流場分布相對于x軸基本對稱,所以實際數(shù)值計算時以x軸為對稱軸取其中一半進行計算以節(jié)約計算時間。

圖1 流場尺寸及網(wǎng)格劃分

床面設置為粗糙壁面邊界(粗糙高度2.5D50),橋墩表面為光滑壁面,側(cè)壁、對稱面以及頂面均設置為對稱邊界,水流入口為速度入口邊界,出口為自由出流邊界。對不同湍流模型的試算表明，k-ωSST模型對湍流渦的模擬最為充分,所以選用該湍流模型。

1.2 臨界剪應力的修正

在平床模擬中可以根據(jù)Shields曲線得到泥沙的起動剪應力,但是在動床模擬中,隨著沖刷的發(fā)展,床面逐漸具有了一定的坡度,泥沙沿斜坡方向的重力分力使泥沙更易于起動。隨著沖刷的發(fā)展,床面坡度會逐漸增大,泥沙在斜坡方向的重力分力逐漸增大,從而加速了泥沙的起動,如果沿用平床試驗得到的Shields曲線計算起動剪應力,則將造成較大的誤差。通常的做法是引入一個修正系數(shù)r對起動剪應力τ0進行修正,因此考慮坡度影響的臨界剪應力可以表示為

τcr=rτ0

(1)

圖2 床面單元的橫向及縱向坡角

關(guān)于r的計算方法,Lane公式[12]僅考慮了垂直于流動方向的橫向坡角α的影響,平行于流動方向的縱向坡角θ設為0(圖2);Ikeda公式[13]引入升阻力比值η,但是仍然沒有考慮縱向坡的影響；Dey[14]提出的公式同時考慮了橫向坡和縱向坡角的影響,并將升阻力的影響考慮到公式中,其表達式為

(2)

式中：φ為泥沙休止角。

根據(jù)式(2)可以畫出r分別隨橫向坡角與縱向坡角的變化曲線(圖3)。在坡度相同的情況下,橫向坡角α對剪應力的影響小于縱向坡角θ,而已有的沖刷數(shù)值模擬中對剪應力的修正往往只考慮了橫向坡角的影響,會造成較大的誤差。

圖3 橫向與縱向坡度對修正系數(shù)的影響

Bihs等[15]將以上幾種方法應用到丁壩沖刷的數(shù)值模擬中并與試驗結(jié)果進行了對比,發(fā)現(xiàn)不考慮升力影響的方法計算獲得的沖刷深度要小于實際的沖刷深度。沖刷的橫向和縱向坡度在沖刷發(fā)展中均變得較大,所以考慮兩個坡度影響的Dey公式是最為精確的修正系數(shù)計算公式，因此選用該公式進行模擬。

在數(shù)值模擬中實現(xiàn)變臨界剪應力方法,關(guān)鍵點在于正確地確定床面單元的縱向坡角θ和橫向坡角α。

圖4 床面單元θ和α的確定方法

每個床面單元的法向矢量可以通過FLUENT的宏函數(shù)F_AREA()獲得,從而可以通過向量運算獲得床面單元與xOz平面交線及與yOz平面交線的方向向量(分別設為a和b)。以縱向坡角θ的求解為例(圖4),計算公式如下:

(3)

(4)

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

1.3 推移質(zhì)輸運方程

推移質(zhì)輸運方程總體上都可以用以下無量綱形式表示[16]:

(5)

(6)

式中:μd為動摩擦因數(shù);Ub為推移質(zhì)泥沙顆粒的平均輸運速度;k為常數(shù),取值范圍為0.4～0.7。

已有的推移質(zhì)輸運方程多為基于非擾動流(無阻流物)的沖刷試驗得到的經(jīng)驗或半理論半經(jīng)驗的公式,實際的橋墩沖刷在橋墩附近會產(chǎn)生湍流渦系結(jié)構(gòu),增加局部沖刷的深度。而渦的增強效應可以通過修正輸運方程中的無量綱沙床剪應力τ*來修正。

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

(7)

式中:ρs為泥沙的顆粒密度;γ0泥沙顆粒的半徑;u為流速;Ω為渦度。

(8)

(9)

作用在泥沙顆粒上的總的無量綱剪應力可以修正為

(10)

1.4 河床動態(tài)地貌模型

數(shù)值模擬中要得到河床面高程的實時變化,可以通過泥沙的連續(xù)性建立河床的動態(tài)地貌模型,其統(tǒng)一的形式可以用Exner方程表示:

(11)

式中:p為泥沙的孔隙率;q為泥沙的單寬體積輸沙率;Es-Ds為泥沙懸移質(zhì)下沉通量與上浮通量之差;z為床面高程;t為時間。為了簡化計算,模擬時暫不考慮懸移質(zhì)輸沙的影響,所以Es-Ds一項忽略不計。

1.5 泥沙坍塌的調(diào)整

橋墩的局部沖刷屬于局部的一種較大的變形,所以,隨著河床面高程的變化,床面的橫向和縱向坡度都會不斷地增大。但實際情況是當坡度大于泥沙的自然休止角(泥沙的內(nèi)摩擦角)后,泥沙會不斷地坍塌調(diào)整,沖坑的側(cè)壁坡度會小于或等于泥沙的休止角。在數(shù)值模擬中,如果只考慮輸沙率對床面高程的影響,則會出現(xiàn)沖坑側(cè)壁坡度大于泥沙休止角的情況,這顯然與實際不符。

為了考慮坍塌的影響,在數(shù)值模擬中,每個時間步網(wǎng)格節(jié)點高程根據(jù)地貌模型調(diào)整后,都需要判斷每個單元平面外法向方向與z軸正向的夾角大小(圖5),若夾角大于φ-2°,則z向坐標大于單元中心z向坐標的節(jié)點(如圖5A節(jié)點)要下移,而小于的要上移(如圖5B、C節(jié)點)。為了保證數(shù)值穩(wěn)定,每步調(diào)整的移動距離為根據(jù)輸沙率計算得出的該節(jié)點高程變化值,只是移動的方向做了調(diào)整。

圖5 泥沙坍塌調(diào)整方法

1.6 UDF程序及網(wǎng)格更新方法

FLUENT中的動網(wǎng)格模型,可以用于模擬流體域邊界隨時間改變的問題。動網(wǎng)格中的UDF宏主要有3種,其中DEFINE_GRID_MOTION宏提供了對節(jié)點和網(wǎng)格最大限度的操作,因此適用于局部沖刷這樣的局部大變形問題。

采用DEFINE_GRID_MOTION宏控制底部邊界各個節(jié)點運動時,每一個時間步都必須執(zhí)行。在每個時間步開始計算時,比較床面(動邊界)各節(jié)點實時剪應力τ與床沙起動臨界剪應力τc,若存在超臨界剪應力(τ-τc>0),該點表現(xiàn)為沖刷,節(jié)點下移;否則表現(xiàn)為靜止,節(jié)點位置不變。在每個時間步某一節(jié)點下移的距離取決于河床動態(tài)地貌模型和時間步長的大小。時間步長的選取原則是:在每個時間步內(nèi),任一網(wǎng)格運動的最大位移不能超過泥沙顆粒的粒徑。

床面的剪應力是隨著沖刷的過程而不斷變化的,FLUENT中通過提取每一網(wǎng)格單元實時剪應力除以單元的面積得到剪應力。在計算每個單元的變臨界剪應力時,有3個角度需要通過編寫程序間接得到,即河床的縱向坡角、橫向坡角和河床面與床面單元外法向的夾角。這3個角度需要通過提取向量并進行相應的運算得到,具體方法參考1.2節(jié)。

推移質(zhì)輸運方程修正時需要計算渦度Ω,而Ω=(Ωx,Ωy,Ωz),這里主要考慮的是Ωy,可以通過FLUENT軟件中C_DUDZ函數(shù)與C_DWDX函數(shù)之差計算出來。

對于沖刷問題,局部運動邊界位移大,這里選用光順法和網(wǎng)格重構(gòu)法結(jié)合的方式來對網(wǎng)格進行重構(gòu)如圖6所示,雖然局部變形較大,但是始終保持了較好的網(wǎng)格質(zhì)量,確保計算結(jié)果的收斂和流場計算的準確度。

2 模擬結(jié)果與試驗結(jié)果[11]對比

2.1 沖坑形態(tài)

Melville[11]試驗研究表明,沖坑發(fā)展試驗中前30 min發(fā)展較為劇烈,而之后沖刷發(fā)展變得緩慢。

30 min時的沖刷深度達到4 cm,為總沖刷深度的75%。為了減少計算消耗,數(shù)值模擬取前30 min進行研究。

圖7 沖刷坑形態(tài)對比

圖7給出了數(shù)值模擬得到的30 min時沖坑的形態(tài),由于數(shù)值模擬以x軸為對稱軸取其中一半進行模擬,所以效果圖中另一半是由鏡像獲得。

從圖7可以看出:①試驗和數(shù)值模擬的30 min時的沖坑最大深度均為4 cm,但是試驗獲得的最大沖刷深度位置從兩側(cè)約70°的位置貫通到了橋墩迎水面正前方,而數(shù)值模擬得到的最大沖刷深度位于兩側(cè)而并沒用貫通到正前方,因此在橋墩正前方的最大沖刷深度約比試驗結(jié)果小25%。造成這一誤差的原因可能為:為了節(jié)約計算時間,以x軸為對稱軸,取一半進行計算,因此x軸附近的網(wǎng)格節(jié)點變?yōu)檫吔?受邊界效應的影響,橋墩正前方位置可能會造成一定的誤差;模型對圓柱形墩正前端湍流脈動的模擬能力所限。②從沖刷等值線以及縱橫向截面圖對比來看,整體分布形態(tài)較為吻合,數(shù)值模擬得到的橋墩上游縱向坡度及橫向坡度約為30°,這與泥沙的自然休止角非常吻合。橋墩下游坡度較緩,約為20°。試驗和模擬的沖坑頂面寬度均為10 cm。

2.2 流場結(jié)構(gòu)

墩周馬蹄渦及墩前流線如圖8所示,無論是墩前沖坑內(nèi)的渦系還是延伸到兩側(cè)的馬蹄渦都與試驗情形非常吻合。墩后側(cè)的渦系在現(xiàn)有文獻中缺乏形象的描述,所以沒有進行對比分析。模擬的尾渦軸線垂直于河床面且渦管直徑由下向上不斷增大,到接近水面的位置有表面的水平渦,這與很多試驗的尾渦描述非常吻合。

圖8 墩周馬蹄渦及墩前流線

不同高程處沖刷坑內(nèi)流速矢量圖如圖9所示,不同高程處的沖坑內(nèi)的流速矢量分布形態(tài)非常相似。在沖坑頂部位置(z=0 cm、z=-1 cm高程處)沒有出現(xiàn)明顯的馬蹄渦系,但在z=-2 cm、-3 cm處出現(xiàn)矢量的明顯分離,這表示有了明顯的馬蹄渦。

以上分析表明,數(shù)值模擬結(jié)果與試驗結(jié)果基本吻合,模擬結(jié)果可靠。

3 不同截面橋墩的沖刷特性

圖9 不同高程處沖刷坑內(nèi)流速矢量

已有文獻對于圓柱形橋墩、方形橋墩以及尖角形的橋墩探討較多,一致的結(jié)論為在相同的阻塞比情況下,沖刷深度由大到小為:方形橋墩、圓柱形橋墩、尖角形橋墩,但是對于流線形的橋墩在已有的沖刷方程和試驗研究中較為少見。Arneson等[19]指出,流線形橋墩截面可能會降低橋墩迎水面的流動分離,從而縮減馬蹄渦和尾渦的強度,降低沖刷的深度。於剛節(jié)[20]提出了一種新型的整體流線形曲面丁壩來改善丁壩附近水流結(jié)構(gòu),進而減少壩頭局部沖刷,通過數(shù)值模擬發(fā)現(xiàn)在相同流量條件下,最大沖刷深度在曲面丁壩附近比在梯形丁壩附近減少了約20%,且沖刷坑范圍也有減小。Al-Shukur等[21]對10種不同形狀的橋墩的試驗研究得出,矩形截面橋墩的沖刷深度最大,流線形橋墩沖刷深度最小。

利用本文的動床模擬方法對尖角形橋墩和流線形橋墩進行了模擬。水流深度、來流速度等均與Melville[11]試驗條件保持一致。尖角形橋墩和流線形橋墩最寬位置的寬度等于圓柱形橋墩直徑,長寬比取為3。

圖11 不同截面橋墩局部沖坑等值線(單位：cm)

數(shù)值模擬得到的沖坑深度隨時間的變化曲線如圖10所示?？梢钥闯?沖刷速度在開始階段較快,10 min后變緩,符合試驗得到的沖刷逐漸變緩的規(guī)律。流線形橋墩在30 min時的沖坑深度比圓柱形橋墩降低約45%,比尖角形橋墩降低約40%。

圖10 沖坑深度隨時間的變化

根據(jù)圖11,不同截面的橋墩在沖刷10 min和沖刷30 min時,無論是沖坑的深度還是沖刷的范圍,圓柱形最大、尖角形次之、流線形最小。因此,從降低局部沖刷的角度來看,流線形橋墩無疑是一種比尖角形橋墩更理想的橋墩形式。

4 結(jié) 論

a. 提出的橋墩局部沖刷過程動態(tài)模擬方法對沖刷的發(fā)展過程、沖坑的形態(tài)以及最大沖刷深度都達到了較為精確的模擬,雖為局部的大變形問題,但是網(wǎng)格質(zhì)量始終較高。

b. 以下幾個方面是動床模擬精度的保證:選用k-ωSST湍流模型,在進行床面剪應力修正的同時考慮縱向和橫向坡度的影響,泥沙輸運方程考慮渦的增強效應,進行泥沙坍塌的調(diào)整及網(wǎng)格的局部加密及重構(gòu)。

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 ] 韋雁機,葉銀燦,吳珂,等.樁周局部沖刷三維數(shù)值模擬[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 ] 祝志文,劉震卿.圓柱形橋墩周圍局部沖刷的三維數(shù)值模擬[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] 於剛節(jié).正態(tài)曲面丁壩附近三維水流及局部沖刷[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.