河床中值粒徑對橋墩沖刷坑演化過程的影響

高 虎，王秋生

(北京工業(yè)大學(xué) 城市與工程安全減災(zāi)教育部重點實驗室，北京 100124)

1 研究背景

水流遇到橋墩時，河道橫截面由于橋墩的存在而減小，會使水流發(fā)生擾動，在橋墩周圍形成非常復(fù)雜的三維流場。伴隨著同樣復(fù)雜的泥沙運動，橋墩周圍發(fā)生局部沖刷使橋墩基礎(chǔ)的埋深減小。這使得橋梁承載力減小，進(jìn)而導(dǎo)致橋墩的坍塌，造成不可估量的損失。根據(jù)Shirolet[1]的研究，美國在過去的30 a中，有超過1 000座橋倒塌，其中大約有60%的橋梁水毀與橋墩的沖刷破壞有關(guān)。易仁彥[2]收集了國內(nèi)在2000年之后15 a的106起橋梁坍塌事故，發(fā)現(xiàn)有半數(shù)以上的破壞是由橋墩的局部沖刷導(dǎo)致的,局部沖刷嚴(yán)重影響了橋梁的安全?，F(xiàn)階段，對于局部沖刷作出準(zhǔn)確的預(yù)測是迫切且重要的。

基于泥沙運動的動理學(xué)理論和流體力學(xué)的橋墩局部沖刷的數(shù)值模擬理論得到了廣泛的應(yīng)用。Olsen等[3]最先運用三維數(shù)值模擬的方法研究了橋墩的局部沖刷，將泥沙的輸運方程耦合到納維斯托克斯方程方程中，進(jìn)而計算出了橋墩的最大沖刷坑深度，印證了橋墩沖刷的三維模擬是可以實現(xiàn)的。Hoffmans和Booij等[4]在假定流體靜壓分布和拋物線解法的情況下，驗證了流速、泥沙濃度的測量值與計算值是一致的。凌建明等[5]通過網(wǎng)格動態(tài)劃分技術(shù)、標(biāo)準(zhǔn)的k-ε湍流模型和流體體積函數(shù)模型對三維的橋墩周圍流場進(jìn)行數(shù)值模擬，得到了橋墩附近的切應(yīng)力分布圖。Khosronejad等[6]研究了URANS模型對各種橋墩 (圓柱形、方形和菱形)的預(yù)測能力，結(jié)合曲線浸入邊界法對k-w紊流閉合的控制方程進(jìn)行了數(shù)值求解，計算結(jié)果表明不同形狀的橋墩在沖刷深度和沖刷速度上是存在差別的，且數(shù)值模擬結(jié)果的精確度也存在差別。李玲等[7]對圓柱擾流的流態(tài)進(jìn)行了研究，得出了圓柱數(shù)量不同，流動流態(tài)不同，且多排圓柱不利于尾流渦結(jié)構(gòu)的穩(wěn)定。祝志文等[8]使用標(biāo)準(zhǔn)k-ε湍流模型在計算流體動力學(xué)軟件Fluent中對圓柱形橋墩周圍的流場進(jìn)行了數(shù)值模擬，并且通過二次開發(fā)Fluent軟件使得沖刷坑的演變形態(tài)可以直觀地觀察。李紹武等[9]使用大渦模擬湍流模型，得出網(wǎng)格的細(xì)化程度對局部沖刷的預(yù)測存在較大影響，并研究了橋墩的直徑對橋墩局部沖刷深度的影響。龍慶[10]通過數(shù)值模擬與王順意等[11]的試驗進(jìn)行結(jié)果對比，最終得到了挾沙水流對局部沖刷演變規(guī)律的影響。李紹武等[9]通過數(shù)值模擬與Roulund等[12]的物理模型試驗結(jié)果對比，最終得到了圓柱直徑對局部沖刷深度的影響。張曙光等[13]通過數(shù)值模擬與Melville等[14]的經(jīng)典沖刷試驗資料進(jìn)行對比驗證，認(rèn)為大渦模擬可以更準(zhǔn)確地預(yù)測橋墩局部沖刷。

綜上所述，橋墩局部沖刷的數(shù)值模擬已經(jīng)有了較多的研究，但是對于水流流速及河床顆粒中值粒徑的改變對沖刷坑的影響問題卻很少涉及。并且大部分模擬的模型都忽略了由于非穩(wěn)態(tài)、非平衡過程出現(xiàn)的擬序結(jié)構(gòu)，將導(dǎo)致最終的結(jié)果偏差較大，而大渦模擬可以很好地解決這個問題。本文將在大渦模擬下結(jié)合水流動力學(xué)及泥沙動理學(xué)對橋墩在水流沖刷下的流場變化及河床形態(tài)變化進(jìn)行數(shù)值模擬。

2 理論分析

2.1 大渦模擬下的控制方程

流體的運動是通過非線性、瞬態(tài)、二階微分方程來描述的，通常使用近似解來代替運動方程的解析解來解決原始的問題。大渦模擬是求解這組方程的一種重要方法。它的原理是通過濾波將湍流運動分解為大尺度運動與小尺度運動，其中大尺度運動直接進(jìn)行計算求解，小尺度運動將會等效替換成為亞格子雷諾應(yīng)力項。這樣會使得運算量減小并簡化計算過程。流體的主要控制方程為連續(xù)方程(1)和納維斯托克斯方程(2)，具體表達(dá)式為：

(1)

(2)

運用渦旋黏性模型和Smagorinsky模型[15]來求解運算τij，即

(3)

將流體的本構(gòu)方程代入動量守恒方程中得到納維斯托克斯方程。本構(gòu)方程的表達(dá)式為

(4)

式中：δij是克羅內(nèi)符號(Kronecker symbol)，當(dāng)i≠j時，δij=0，當(dāng)i=j時，δij=1；sij是應(yīng)變速率張量。

2.2 泥沙運動方程及底床變形

泥沙顆粒在水流作用下離開床面的過程，在泥沙運動學(xué)中被稱為泥沙的起動。泥沙顆粒從床面起動與其受力狀態(tài)是緊密聯(lián)系的[16]。最重要的假定是根據(jù)統(tǒng)計學(xué)理論，單個顆粒的行為可以替代顆粒群體的統(tǒng)計規(guī)律。這里采用Mastbergen和Van den Berg的經(jīng)驗輸沙模型[17]。泥沙處于臨界狀態(tài)時，作用在泥沙顆粒上的剪切力與有效重力的比，即臨界Shields數(shù)是判斷泥沙起動的重要參數(shù)。當(dāng)泥沙顆粒實際Shields數(shù)大于臨界Shields數(shù)時，泥沙起動。采用Soulsby-Whitehouse公式[18]計算無量綱臨界Shields數(shù)θcr，即

(5)

式中：β為床面的坡度；φ為泥沙的休止角；d*為無量綱化后的顆粒粒徑，計算公式為

(6)

式中：‖g‖為重力加速度；ρ為流體的密度；ρs為沉積物的密度；d為顆粒粒徑。

局部Shields數(shù)θ根據(jù)床層剪切應(yīng)力τ計算得出，即

(7)

由于泥沙顆粒與水的密度存在差異，在重力和水流作用下，泥沙顆粒與水體之間會發(fā)生相對運動，表現(xiàn)為水流對泥沙攜帶和泥沙自身重力的沉降。河床泥沙在水流中的存在形態(tài)分為推移質(zhì)和懸移質(zhì)2種。推移質(zhì)基于現(xiàn)象的描述通常是指那些在床面以滾動、滑動、跳躍為主要運動形式的泥沙顆粒。懸移質(zhì)粒徑相對較小，或水流強度相對較大，泥沙顆粒往往彌散于水體中以懸浮狀態(tài)運動。對于推移質(zhì)來說，通常使用Meyer-Peter和Muller公式[19-20]，即

Φ=8(θ-θcr)3/2。

(8)

式中Φ為推移質(zhì)無量綱輸沙率，定義為

(9)

式中qb為推移質(zhì)單位寬度輸沙率。

對于懸移質(zhì)來說，其輸移由基于紊流擴散理論和Fick定律的對流擴散方程來描述[21]，即

(10)

式中：cp為懸移質(zhì)的濃度；D為泥沙的紊流擴散張量；us為懸移質(zhì)的移動速度，根據(jù)式(11)及Soulsby公式[19]計算得到，即

(11)

2.3 Tru-VOF自由液面

Tru-VOF(True Volume of Fluid)自由液面是通過在歐拉網(wǎng)格的表面進(jìn)行追蹤，從而判定網(wǎng)格的屬性進(jìn)而確定邊界。其核心是先定義一個流體體積函數(shù)，即目標(biāo)流體體積所占目標(biāo)網(wǎng)格體積的大小。由于在不同流體的交界面處仍然使用一套動量的方程，則通過這一函數(shù)便可以對網(wǎng)格的屬性進(jìn)行研究并確定自由面[22]。其核心思想是[23]在整個網(wǎng)格空間中，所有網(wǎng)格的體積分?jǐn)?shù)均為1，定義第n相流體的流體體積函數(shù)Wn，則當(dāng)此函數(shù)Wn=0時，則該網(wǎng)格內(nèi)不存在該流體；當(dāng)Wn=0～1時，該網(wǎng)格內(nèi)存在該流體并且存在其他流體或者存在空氣氣泡，此時該網(wǎng)格應(yīng)為不同流體或者氣液交界的界面網(wǎng)格；當(dāng)Wn=1時，則該網(wǎng)格被該流體占滿。Wn表達(dá)為時間和空間的函數(shù)，應(yīng)該滿足在該系統(tǒng)內(nèi)dWn/dt=0的輸運關(guān)系，即

式中u、v、w為速度在3個方向的分量。

求解式(12)，便得到每個網(wǎng)格內(nèi)的流體分布情況，從而可以確定自由界面的分布情況。

Tru-VOF方法對氣體區(qū)域不考慮其壓力、速度等變化，液面設(shè)定為相同的壓力。這樣可以更好地捕捉到氣液分界面，計算時間縮短，精度也有較大提升。

3 模型的驗證

3.1 計算條件

本文使用Melville等[14]的沖刷試驗結(jié)果與本文所做的沖刷模擬結(jié)果進(jìn)行驗證對比。模型的簡圖見圖1，具體條件為：計算域尺寸為0.6 m×0.2 m×0.2 m，在計算域內(nèi)設(shè)置網(wǎng)格的尺寸為0.002 m×0.002 m×0.002 m。選用直徑為5.06 cm的圓柱位于計算域中央。

圖1 模型簡圖Fig.1 Diagram of experimental layout

計算域底部兩側(cè)分別有一層長度為0.01 m的固定床，其余部分鋪滿均勻無黏性土顆粒，兩者的厚度都取0.06 m，土顆粒的中值粒徑d50=0.385 mm，密度ρs=2 650 kg/m3，休止角φ=32°，床面的粗糙度取為河床土顆粒的中值粒徑的2.5倍。床面向上是流動的溫度為20 ℃的水，水的密度ρ=1 000 kg/m3，黏度μ=0.001 Pa·s，水體的初始流體速度取0.25 m/s，水的深度為0.05 m。模型中設(shè)置重力加速度為9.81 kg/m3，模擬的總時間為100 min。

計算域的底部為壁面邊界條件；上部為壓力邊界條件，設(shè)置流體占比量(Fluid friction)為0；兩側(cè)邊界采用對稱邊界條件；入口處采用速度邊界條件，設(shè)置入口速度大小為0.25 m/s，流體的高度為0.05 m；出口處采用流出邊界條件。

3.2 模擬結(jié)果與試驗結(jié)果比對

圖2給出了水平狀態(tài)下流場的數(shù)值模擬結(jié)果與水槽試驗結(jié)果的對比[24]。可以觀察到水流在觸碰到圓柱形橋墩之前發(fā)生了流場的擴散。在此處流線變得松散，并向圓柱橋墩兩側(cè)發(fā)散使得圓柱前側(cè)的水流沿x方向的速度急速下降甚至出現(xiàn)碰撞回流，圓柱橋墩兩側(cè)的水流流線則變得更加密集，流速在這里急速上升。由于橋墩的存在，橋墩后側(cè)的水流流線變得特別稀疏，水流流速也很小，加之水壓力不平衡，在此處會出現(xiàn)對稱的尾渦結(jié)構(gòu)，并在橋墩背水側(cè)延續(xù)很長的距離才會逐漸消散。

圖2 水平狀態(tài)下t=2 min時的流場對比Fig.2 Comparison of flow field between test andsimulation at t=2 min in horizontal state

圖3給出了垂直狀態(tài)下的流場對比。在數(shù)值模擬進(jìn)行到80 min時，橋墩前側(cè)已經(jīng)形成沖刷坑，水體在觸碰到橋墩時會變成向上和向下的兩股水流：向下形成了下潛流，在橋墩前側(cè)底部接觸到河床又會發(fā)生水流方向的改變形成馬蹄形旋渦；向上形成了上升流，在水面處形成了涌波，數(shù)值模擬與試驗結(jié)果基本吻合。在數(shù)值模擬結(jié)果圖中，給出的橋墩前側(cè)和后側(cè)的沖刷坑的具體形態(tài)也比較符合實際，但是橋墩后側(cè)由于尾渦的作用，水流的流向在此二維圖中無法展現(xiàn)。

圖3 垂直狀態(tài)下達(dá)到平衡狀態(tài)時的流場對比Fig.3 Comparison of flow field at equilibrium invertical state between test and simulation

圖4給出了在局部沖刷達(dá)到穩(wěn)定后，橋墩周圍局部沖刷深度的高程?？梢钥吹綌?shù)值模擬的高程分布是與試驗結(jié)果較吻合的。

圖4 平衡狀態(tài)下的沖刷坑形態(tài)Fig.4 Shape of scour pit at equilibrium

表1對比了試驗結(jié)果和不同數(shù)值模擬結(jié)果的最大沖刷深度，可以發(fā)現(xiàn)本文采用的大渦模擬方法在與試驗結(jié)果的比對中誤差最小,證明該數(shù)值計算方法得到的模擬結(jié)果比較準(zhǔn)確。

表1 不同數(shù)值模擬模型最大沖刷深度的結(jié)果比對Table 1 Comparison of the results of maximum scourdepth among different numerical simulation models

4 數(shù)值試驗過程、結(jié)果和討論

4.1 數(shù)值試驗方案

影響橋墩沖刷的因素眾多，其中水流因素和泥沙因素是最重要的兩類。保持模型的數(shù)據(jù)不變，只改變水流流速與河床的顆粒中值粒徑大小。通常水流的流速在0.8 m/s左右，本文也選用這個區(qū)間的水流流速來進(jìn)行研究；考慮到本文研究的是砂土對沖刷的影響，因此選用的河床顆粒大小應(yīng)包括細(xì)砂(顆粒粒徑≤0.25 mm)、中砂(0.25 mm<顆粒粒徑≤0.5 mm)和粗砂(0.5 mm<顆粒粒徑)，具體模型設(shè)定見表2。通過控制變量法來觀測某一條件的改變對整個沖刷過程的影響。

表2 模擬工況參數(shù)設(shè)置Table 2 Setting of test condition parameters

4.2 觀測點的沖刷深度演化規(guī)律

在局部沖刷剛開始時，沖刷坑的形成是迅速且劇烈的。在迎水側(cè)由于水流碰到橋墩形成下潛的水流進(jìn)而形成馬蹄渦，導(dǎo)致水流流速增大，床面剪切應(yīng)力增加，使得沖刷坑深度迅速增加；在橋墩的兩側(cè)，橋墩對水流的壓縮效應(yīng)導(dǎo)致水流流速和床面剪切應(yīng)力迅速增大，也使得沖刷坑深度迅速增加；在背水側(cè)由于開始時的橋墩阻擋，水流流速降低，床面剪切應(yīng)力減小，而橋墩前側(cè)及兩側(cè)的水流沖刷效應(yīng)使水流中帶有大量的泥沙顆粒，這些影響因素共同導(dǎo)致在沖刷過程的初期,背水側(cè)泥沙出現(xiàn)堆積現(xiàn)象。

隨著時間的推移，迎水側(cè)與兩側(cè)的沖刷坑深度增長速率逐漸變緩；并在60 min左右開始趨于穩(wěn)定；背水側(cè)出現(xiàn)卡門渦街致使床面的剪切應(yīng)力增大，床面的泥沙堆積又變?yōu)榫植繘_刷，最終的沖刷深度大約只有最大沖刷深度的一半左右。

圖5為不同觀測點位置且只有水流流速改變的情況下局部沖刷深度隨時間的變化。在局部沖刷剛開始時，隨著水流的流速增大，橋墩的迎水側(cè)、橋墩最大沖刷深度處和橋墩兩側(cè)的局部沖刷深度增加且沖刷深度的增長速率也增加；隨著時間的進(jìn)行，這3處最終的變化是流速越大，沖刷深度越大。背水側(cè)開始時出現(xiàn)泥沙堆積現(xiàn)象，在4～6 min時泥沙堆積深度出現(xiàn)極值，隨后便出現(xiàn)快速的沖刷現(xiàn)象，但是最終的沖刷深度卻遠(yuǎn)低于同時間其他3處的沖刷深度。對于流速不同的流場，流速增大后導(dǎo)致床面剪切應(yīng)力增大，使泥沙顆粒更容易達(dá)到起動條件，從而造成最終沖刷深度變大。

圖5 不同水流速度下的局部沖刷深度隨時間的變化Fig.5 Variation of local scour depth against time at different flow rates

圖6顯示了不同觀測點位置上隨著時間的發(fā)展河床顆粒粒徑的大小對局部沖刷深度的影響。在各個觀測點上，當(dāng)顆粒中值粒徑從0.145 mm增加到0.385 mm時，沖刷深度有顯著增大；當(dāng)顆粒中值粒徑從0.385 mm增加到0.625 mm時，二者在各個時刻下的沖刷深度相差不大；當(dāng)顆粒中值粒徑從0.625 mm增加到0.815 mm和1.015 mm時，在同一時刻下沖刷深度顯著減小?？梢缘贸鰧τ诤胁煌w粒粒徑的河床來說，當(dāng)粒徑減小時泥沙顆粒更容易起動并最終得到較大的沖刷深度。但是粒徑的大小超過某個數(shù)值時，粒徑的減小會使沖刷深度變小。結(jié)合本文第二節(jié)中的式(10)計算以及趙凱[25]的試驗研究可以得出，這是由于顆粒粒徑過小導(dǎo)致顆粒之間的黏聚力增大，從而在沖刷過程中起到限制沖刷的作用。

圖6 不同河床中值粒徑下的局部沖刷深度隨時間的變化Fig.6 Variation of local scour depth against time with different median sizes of particles

5 結(jié) 論

本文應(yīng)用計算流體力學(xué)方法模擬了橋墩周圍的水流流態(tài)以及河床的形態(tài)變化，研究了水流流速和河床中值粒徑2個變量的依次改變對橋墩沖刷的影響，主要得出以下幾點結(jié)論：

(1)水流的流動由于橋墩的阻擋在橋墩周圍形成了復(fù)雜的流場結(jié)構(gòu)，并改變了橋墩周圍的河床形態(tài)。最大沖刷深度的位置位于迎水側(cè)兩側(cè)約60°處。

(2)由于橋墩的阻擋，橋墩的背水側(cè)在模擬初始時段發(fā)生泥沙堆積現(xiàn)象，但隨著模擬的進(jìn)行仍會發(fā)生泥沙的沖刷現(xiàn)象。最終沖刷深度只能達(dá)到最大沖刷深度的一半左右。

(3)河床泥沙的顆粒中值粒徑減小會導(dǎo)致沖刷深度增大，但是當(dāng)顆粒中值粒徑達(dá)到一定范圍內(nèi)，由于顆粒之間的黏聚力增大反而導(dǎo)致沖刷深度減小。

綜上所述，使用大渦模擬來模擬水流紊流狀態(tài)和泥沙沖刷是比較準(zhǔn)確的方法?；诟呕圻M(jìn)行的模擬試驗結(jié)果可以作為參考方案，為類似工程提供相關(guān)指導(dǎo)。