李紹武，楊 航

(天津大學(xué) 水利工程仿真與安全國家重點實驗室，天津 300072)

墩柱繞流現(xiàn)象所帶來的墩柱周邊沖刷問題是實際工程所不容忽視的一大問題，國內(nèi)外有許多研究成果，采用數(shù)值模擬方法的研究工作起源于墩柱周圍的流場模擬。1992年，Kobayashi[1]通過將二維離散渦模型進(jìn)行擴展，提出了三維渦段模型，嘗試了在振蕩流的情形下模擬圓柱周圍的三維流場分布，模擬出了尾渦的不穩(wěn)定特征，并和實測的流速場分布吻合較好。1993年，Olsen 和 Melaaen[2]首次將數(shù)值模擬運用于計算三維圓柱周圍局部沖刷問題，他們通過引入泥沙輸運方程，模擬了清水沖刷條件下圓柱周圍無粘性沙局部沖刷過程的初期發(fā)展階段，計算結(jié)果與實驗沖刷型態(tài)能夠較好地吻合。1998年，Olsen和Kjellesvig[3]在1993年模型的基礎(chǔ)上進(jìn)一步開展研究，模擬的沖刷過程不再僅局限于沖刷的初始階段，而是擴展到了沖刷的全過程，最終計算所得的清水條件下圓柱周圍的沖刷深度與經(jīng)驗公式計算結(jié)果吻合良好。2003年，Catalano和Wang[4]在高雷諾數(shù)的條件下運用大渦模擬研究了墩柱周圍的繞流問題，通過與RANS紊流模型所得計算結(jié)果進(jìn)行比對后發(fā)現(xiàn)，大渦模擬所得的計算結(jié)果能夠更好地反應(yīng)高雷諾數(shù)情況下拖曳系數(shù)減小和邊界層分離推后的現(xiàn)象，但是由于大渦模擬對計算資源的要求較高，并未得到廣泛運用。2005年，Roulund和Sumer[5]基于SST k-ω紊流模型模擬了圓柱周圍的三維沖刷問題，較好地模擬出了床面沙紋的產(chǎn)生以及沖刷坑的形成過程，模擬結(jié)果中的平衡沖刷深度與實驗值吻合良好。

由上述國內(nèi)外研究進(jìn)展可以看出，前人在水流作用下圓柱周圍繞流流場和局部沖刷上已取得了顯著進(jìn)展，但仍存在著一些問題。一方面，相對于其他紊流模型而言， LES可以更好地模擬圓柱繞流流場，但是，由于LES模擬的流場結(jié)果依賴于網(wǎng)格尺度，其對于局部沖刷結(jié)果產(chǎn)生的影響尚不明確。另一方面，前人的研究成果大多局限于小尺度的圓柱，不同尺度圓柱周圍的沖刷特性存在何種差異也需要進(jìn)一步探討。本文運用FLOW-3D軟件基于大渦模擬模型，對不同直徑的圓柱周邊的流態(tài)及沖刷過程進(jìn)行數(shù)值模擬探討，研究圓柱尺度變化及網(wǎng)格大小對沖刷坑深度及形態(tài)的影響。

1 計算模型及其數(shù)值理論

1.1 基本控制方程及模型理論

LES過濾后的連續(xù)性方程

(1)

LES過濾后的N-S方程

(2)

(3)

渦旋粘性模型

(4)

(5)

Smagorinsky模型

(6)

LS=CS(δxδyδz)1/3

(7)

(8)

在沖刷模型采用基于Mastbergen和Van den Berg[6]的經(jīng)驗輸沙模型。首先計算無量綱數(shù)d*,i

(9)

式中：di為第i種泥沙的顆粒直徑；ρf為流體的密度；ρi為第i種泥沙的密度；||g||為重力加速度的范數(shù)；μf為流體的動力粘滯系數(shù)。

根據(jù)Soulsby-Whitehouse公式[7]計算θcr,i

(10)

根據(jù)泥沙休止角對臨界謝爾茲數(shù)進(jìn)行坡度修正[7]

(11)

式中：ψ為水流方向與逆坡方向的夾角；β為床面的坡度；φi為第i種泥沙的休止角。

本地謝爾茲數(shù)通過底部切應(yīng)力τ求解

(12)

泥沙顆粒的挾帶抬升速度ulift,i由下式計算[6]

(13)

式中：αi為挾帶系數(shù)，本文中取0.018[6]；ns為床面的外法線方向。

本文使用的泥沙沉降速度usettling,i計算公式來自Soulsby[7]

(14)

式中：vf是流體的運動粘滯系數(shù)。

推移質(zhì)輸沙率采用Meyer-Peter的計算模型[8]，公式為

(15)

(16)

式中：Φi為無量綱的推移質(zhì)輸送率；βi為推移質(zhì)輸沙系數(shù)，取為8[9]；cb,i是第i種泥沙在床面泥沙中所占的體積分?jǐn)?shù)。

懸移質(zhì)泥沙濃度采用輸運方程來求解，方程為

(17)

式中：Cs,i是第i種泥沙的懸移質(zhì)質(zhì)量濃度，即單位體積的水沙混合物中泥沙的質(zhì)量；us,i為懸移質(zhì)速度；D為擴散系數(shù)。

懸移質(zhì)速度us,i即可按下式計算

(18)

1.2 計算域網(wǎng)格劃分

由于模型計算過程頗為耗時，因此在計算域的選擇上，在不影響正常計算結(jié)果的前提下，盡量選取較小的計算區(qū)域。在正式計算前通過一段40 m長的數(shù)值水槽獲得充分發(fā)展的紊動水流，并通過熱起動將充分發(fā)展的紊流條件作為計算域的入流邊界條件。沖刷模型中，圓柱中心位于距入流邊界8D處，距離出流邊界10D，左右邊界和圓柱中心的距離均為5D。計算區(qū)域設(shè)置及網(wǎng)格劃分如圖1所示。

1-a 平面圖 1-b X-Z剖面圖圖1 計算域設(shè)置及網(wǎng)格劃分Fig.1 Setup of computational domain and mesh discretization

1.3 模型參數(shù)及邊界條件設(shè)置

分別選用直徑為0.1 m，0.2 m，0.3 m及0.5 m的圓柱進(jìn)行模擬，模型底部鋪設(shè)均勻床砂，粒徑d為0.26 mm，密度ρs為2 650 kg/m3，休止角α為32°，床面粗糙度取2.5倍中值粒徑。入口處和出口處均布置一個與床砂等高的擋水護(hù)板，為保證其與床砂的交界處平穩(wěn)過渡，設(shè)置其表面粗糙度為0.65 mm。床砂區(qū)上層為水，水體溫度取20°C，水深為0.4 m，斷面平均流速為0.46 m/s。重力加速度取9.81 m/s2。

為在保證精度的情況下減少網(wǎng)格總數(shù)以減少計算量，模型采用嵌套網(wǎng)格，共分為三層，相鄰層網(wǎng)格的尺寸比例為1∶2。計算域的底部邊界采用壁面邊界條件，法向流速為零；側(cè)邊界和頂部邊界采用對稱邊界條件，此邊界處流量為零，流體剪切應(yīng)力為零；入流邊界采用網(wǎng)格重疊邊界條件，將充分發(fā)展的斷面紊動流速作為入流條件；出流邊界采用壓力邊界，控制出流邊界的水深。

2 模型驗證

2.1 水流模型驗證

參照Roulund[5]的物理模型實驗，建立光滑床面下的水流數(shù)值模型。模型中，圓柱直徑D為0.536 m，水深h為0.54 m，平均入流速度V取0.326 m/s，并經(jīng)一段40 m長的無墩柱水槽的模擬，獲得滿足對數(shù)流速率的入流邊界水流流速分布。

圖2和圖3分別展示了墩柱周圍水流的x向流速及z向流速模擬結(jié)果與實測數(shù)據(jù)的對比，選取的水深位置分別為距床面0.005 m，0.01 m，0.02 m，0.05 m，0.1 m及0.2 m處，可以看出在不同的水深上，模擬的x方向流速和z方向流速與實測數(shù)據(jù)均吻合較好。

2-a z=0.005 m 2-b z=0.01 m

2-c z=0.02 m 2-d z=0.05 m

2-e z=0.1 m 2-f z=0.2 m圖2 墩柱周圍x向流速模擬值與實驗值對比Fig.2 Comparison of modeled velocity in x direction around the pile with measured results

3-a z=0.005 m 3-b z=0.01 m

3-c z=0.02 m 3-d z=0.05 m

3-e z=0.1 m 3-f z=0.2 m圖3 墩柱周圍z向流速模擬值與實驗值對比圖Fig.3 Comparison of modeled velocity in z direction around the pile with measured results

2.2 泥沙沖刷模型驗證

參照Roulund[5]的物理模型實驗，建立泥沙沖刷數(shù)值模型。模型中，圓柱直徑D為0.1 m，水深h為0. 4 m，平均入流速度V取0.46 m/s，模擬時長為60 min。

圖4-a是不同時刻圓柱迎水側(cè)沖刷深度模擬值與實測值的對比。由圖可以看出，在沖刷的初始階段，沖刷坑深度的模擬值與實驗值相比略低，但隨著時間的推移，在沖刷過程的中段直至接近平衡，模擬值和實驗值基本保持一致。

圖4-b是不同時刻圓柱背水側(cè)沖刷深度模擬值與實測值的對比。根據(jù)圖中結(jié)果，背水側(cè)和迎水側(cè)類似，在沖刷初始階段模擬值低于實驗值，但在沖刷過程的中后期，二者大體吻合。

4-a 迎水側(cè) 4-b 背水側(cè)圖4 墩柱迎水側(cè)和背水側(cè)沖刷深度模擬值與實驗值對比圖Fig.4 Comparison of modeled scour depth on the upstream side and rear side of the pile with measured results

圖5 沖刷坑縱向剖面模擬值與實驗值對比圖Fig.5 Comparison between modeled results of longitudinal scour profile with experimental results

圖5是沿圓柱對稱面的沖刷坑縱剖面圖。從圖中沖刷坑剖面的輪廓線中可以看出：在圓柱的迎水側(cè)，沖刷坑的最大深度以及其總體輪廓都能與實驗數(shù)據(jù)良好契合；而在圓柱的背水側(cè)，沖刷坑的最大深度和實驗數(shù)據(jù)吻合良好，但離開圓柱后，模擬所得的沖刷坑較實驗測得的沖刷坑略深一些。

綜上所述，在時間上，模擬值在沖刷過程的初始階段與實驗值相比較低，而在沖刷的中后期可以與實驗值良好契合；在空間上，模擬的沖刷坑在圓柱迎水側(cè)與沖刷坑實測結(jié)果基本吻合，而在背水側(cè)圓柱根部吻合較好，離開圓柱模擬值偏大?？傮w而言，所建立的泥沙沖刷模型能夠較好地與實際情況相吻合。

3 沖刷模擬結(jié)果及分析

保持模型的水深h為0.4 m，平均入流速度V為0.46 m/s，均質(zhì)床砂粒徑d為0.26 mm，砂密度ρs為2 650 kg/m3，休止角α為32°，床面粗糙度為0.65 mm不變，圓柱表面設(shè)為光滑，直徑分別取為0.1 m，0.2 m,0.3 m及0.5 m，計算域的大小及沙床厚度進(jìn)行相應(yīng)調(diào)整。圖6是各組模型在t=1h時沖刷坑沿水流方向?qū)ΨQ面的剖面圖。

6-a D=0.1 m 6-b D=0.2 m

6-c D=0.3 m 6-d D=0.5 m圖6 t=1 h時沖刷坑縱向剖面圖Fig.6 Longitudinal scour profile after 1 hour scour

由于沖刷模型的計算量巨大，尤其是對于大墩柱沖刷模型，因此D=0.1 m,0.2 m及0.3 m時，模型的模擬時長為1.5 h，D=0.5 m時，模擬時長取為1 h。之后，統(tǒng)一運用指數(shù)函數(shù)曲線進(jìn)行非線性擬合并推算至t=5 h處，再進(jìn)行比對分析。

在沖刷的初始階段，無論是迎水側(cè)或是背水側(cè)，沖刷坑的發(fā)展過程都比較迅速，而隨著時間的推移，迎水側(cè)和背水側(cè)的沖刷速率都逐漸趨于平緩，逐漸達(dá)到平衡。在馬蹄渦的影響下，圓柱的迎水側(cè)產(chǎn)生持續(xù)且大范圍的泥沙起動，進(jìn)而產(chǎn)生沖刷坑并不斷擴大。從圓柱迎水側(cè)帶走的泥沙起初堆積在圓柱的后側(cè)，但隨著沖刷過程的持續(xù)進(jìn)行，圓柱迎水側(cè)的沖刷坑深度不斷增加，與此同時，沖刷坑形態(tài)圍繞圓柱向著背水側(cè)不斷發(fā)展，并使得原本堆積于圓柱后側(cè)的泥沙向著下游繼續(xù)推移。

對于直徑為0.5 m的圓柱，當(dāng)x和y方向網(wǎng)格尺寸分別為0.04、0.03和0.02 m時，計算所得的迎水側(cè)和背水側(cè)的沖刷深度出現(xiàn)了一定程度上的差異，如圖7所示。

7-a 迎水側(cè) 7-b 背水側(cè)圖7 網(wǎng)格尺寸造成的沖刷深度差異Fig.7 Discrepancies of scour depth caused by grid size

圖7結(jié)果顯示網(wǎng)格尺寸對計算結(jié)果有一定的影響，但隨著網(wǎng)格的細(xì)化，計算結(jié)果的差異呈收斂趨勢?？紤]到再細(xì)化網(wǎng)格模型計算量過大(采用0.02 m的網(wǎng)格尺寸時，全計算域網(wǎng)格數(shù)達(dá)300余萬，32核工作站模擬1 h花費35 d)，以0.02 m網(wǎng)格尺寸所得計算結(jié)果近似代表實際值。

圖8和圖9分別是各組直徑的墩柱迎水側(cè)和背水側(cè)處的沖刷深度隨時間的變化趨勢。

圖8 各組直徑的圓柱迎水側(cè)沖刷深度圖Fig.8 Scour depths at the upstream side of the pile for different pile diameter圖9 各組直徑的圓柱背水側(cè)沖刷深度圖Fig.9 Scour depths at the rear side of the pile for different pile diameter

從圖中可以看出，在沖刷速率隨時間的變化特性上，迎水側(cè)和背水側(cè)大體一致：在沖刷的初始階段，圓柱的兩側(cè)均有較大的沖刷速率，隨著沖刷過程的進(jìn)行，兩側(cè)的沖刷速率均逐漸減緩并最終趨于零，此時圓柱周圍的沖刷坑趨于平衡狀態(tài)，從圖中可以直觀看出，各組模型在經(jīng)過不同沖刷時間后，均呈現(xiàn)沖刷速率趨緩，達(dá)到平衡的趨勢。

對于同一直徑的圓柱而言，在沖刷深度上，整個沖刷過程中圓柱迎水側(cè)的沖刷深度始終大于背水側(cè)；在平衡時間上，圓柱上下游兩側(cè)沖刷到達(dá)平衡狀態(tài)所需的時間大體一致。對于不同直徑的圓柱而言，隨著圓柱直徑的增大，圓柱周圍的沖刷坑達(dá)到平衡狀態(tài)的時間不斷延長，相應(yīng)地，圓柱迎水側(cè)和背水側(cè)的沖刷深度也隨著圓柱直徑的增加而增加，然而二者之間并非簡單的線性關(guān)系。由圖9可以看出，雖然在最終平衡狀態(tài)下，直徑0.5 m圓柱的背水側(cè)沖刷深度大于直徑0.3 m圓柱，但在沖刷過程的前期階段，前前者的沖刷深度卻始終小后者。造成這種現(xiàn)象的原因在于，在圓柱背水側(cè)發(fā)生沖刷作用的同時，從圓柱迎水側(cè)沖刷下來的泥沙在水流的挾帶作用下，在背水側(cè)產(chǎn)生了堆積，導(dǎo)致較大直徑的圓柱在沖刷過程的初始階段，其背水側(cè)的沖刷深度可能小于較小直徑圓柱。

圖10 各組直徑的圓柱對稱側(cè)沖刷深度圖Fig.10 Scour depths at lateral side of the pile for different pile diameter

圖10是各組直徑的墩柱橫向兩側(cè)的沖刷深度隨時間的變化趨勢?？梢钥闯?，不同直徑的圓柱，其橫向兩側(cè)的沖刷過程與上下游兩側(cè)頗為相似，達(dá)到平衡狀態(tài)時的沖刷深度也大體相同，但沖刷過程與背水側(cè)有所不同。

表1給出了按照各組模型的條件參數(shù)，根據(jù)Breusers[10]公式、中國公路工程水文勘察設(shè)計規(guī)范[11]、Sheppard[12]公式計算所得的最大沖刷深度，并與擬合擬合平衡深度數(shù)據(jù)進(jìn)行比對。

從表中可以看出，擬合平衡沖刷深度結(jié)果與Breusers給定的預(yù)測最大沖刷深度較為接近。

表1 各計算公式所估算的最大沖刷深度Tab. 1 Maximum scour depth estimated using different formula

注：h=0.4 m,U=0.46 m/s

表2記錄了各組模型經(jīng)擬合到達(dá)平衡狀態(tài)時，迎水側(cè)，背水側(cè)及對稱側(cè)的沖刷深度。

從表2中可以看出，對于不同直徑的圓柱，其到達(dá)平衡狀態(tài)時，背水側(cè)沖刷深度與迎水側(cè)沖刷深度之比以及橫向兩側(cè)沖刷深度與迎水側(cè)沖刷深度之比均保持為一常數(shù)，前者約為0.7，而后者在0.95～1.0之間?？梢缘贸鰧τ诓煌膱A柱直徑，圓柱背水側(cè)及橫向兩側(cè)的平衡沖刷深度與迎水側(cè)的平衡沖刷深度始終保持著特定的比例關(guān)系。

表2 各組模型到達(dá)平衡狀態(tài)時的沖刷深度Tab.2 Modeled results of scour depth at equilibrium status

圖11 各組直徑圓柱迎水側(cè)沖刷深度與直徑比值S1/D對比圖Fig.11 Variation of S1/D with diameter D

另，從表2和圖11可以發(fā)現(xiàn)，雖然圓柱迎水側(cè)的沖刷深度隨著圓柱直徑的增大而增大，然其增長率隨著圓柱直徑的增大而逐漸減小，從各組模型迎水側(cè)的平衡沖刷深度與圓柱直徑的比值可以看出，相對于直徑為0.1 m的圓柱，直徑為0.2 m的圓柱的S1/D值下降了24.6%，而直徑為0.3 m和0.5 m的圓柱，其S1/D則分別下降了32.6%和44.7%。

4 結(jié)論

本文利用FLOW-3D三維模擬軟件中的紊流模型和泥沙沖刷模型對不同尺度的圓柱局部沖刷過程進(jìn)行了數(shù)值模擬。計算結(jié)果表明，對于不同直徑的圓柱，其迎水側(cè)、背水側(cè)以及橫向兩側(cè)的平衡沖刷深度始終保持著固定的比例關(guān)系。通過建立不同尺度的圓柱沖刷模型，開展了圓柱局部沖刷數(shù)值進(jìn)行模擬，其中最大的圓柱直徑達(dá)到了0.5 m，計算結(jié)果表明，圓柱周圍的局部沖刷深度隨圓柱直徑的增大而增大，但其與圓柱直徑的比值則隨著圓柱直徑的增大而減小。