齊梅蘭，石粕辰

（1.北京交通大學 土木建筑工程學院，北京 100044；2.北京交通大學 結構風工程與城市風環(huán)境北京市重點實驗室，北京 100044）

1 研究背景

水下結構物的局部沖刷問題，具有水流結構及泥沙輸運隨三維空間和時間尺度變化的復雜性，數(shù)值模擬研究受到廣泛關注[1-3]，以期深入了解沖刷機理并獲得準確的沖刷預測。根據(jù)物理機制，描述局部沖刷的數(shù)學方程包括水流運動、泥沙運動和床面變形方程，其中泥沙運動狀態(tài)很復雜，包括泥沙臨界起動條件以及推移質運動、懸移質運動模式識別等，如考慮不周則可能影響局部沖刷的數(shù)值模擬結果。

天然河流紊流方程瞬態(tài)數(shù)值解的直接模擬甚至大渦模擬因需占用大量的計算資源[4-5]，工程中多用雷諾平均（RANS）方程的數(shù)值模擬。非恒定的雷諾平均方程（URANS）改進了恒定的RANS方程在捕捉復雜渦流能力方面的缺陷，提高了計算精度[6-7]。RANS方程中的紊流脈動項多作線性渦黏應力模型處理，針對渦黏系數(shù)又派生了包括k-ε在內的多種模型。其中，標準k-ε模型對近壁區(qū)和逆壓梯度區(qū)流動的模擬存在不足，重整化群RNGk-ε模型考慮了黏性渦的尺度，在模擬渦流特性方面有一定優(yōu)勢[8-9]。

受水流作用，床面泥沙可呈推移、懸移或二者兼之的運動狀態(tài)。在局部沖刷的數(shù)值模擬中，Olsen等[10]僅考慮了懸移質輸運，更多的學者則只考慮推移質輸運[11-14]，也有考慮推移質和懸移質即全沙輸運的[15-17]。墩柱局部沖刷分為清水沖刷和動床沖刷兩種條件（分別表示床沙未起動和起動的均勻來流條件），Burkow等[12]對清水沖刷條件采用了推移質泥沙運輸；而DIXEN等和JIA等在動床條件下也僅考慮了推移質，但引入了柱前渦流[13]或脈動流速[14]對水流時均切應力做了增強性的修正；Alemi等[18]即使對于清水沖刷也采用了全沙輸運模型；Baykal[15]根據(jù)其計算則認為如不考慮懸移質輸運，墩前最大沖刷深度減小50%，無論是否進行渦流及脈動修正；Roulund等[11]在其計算條件下認為懸移質不重要；Radice等[19]的試驗盡管是動床沖刷，也稱沖刷過程中懸移質輸運不重要?？梢娔壳熬植繘_刷模擬中，對泥沙輸運物理模式的合理選用未取得共識，作者認為，這一問題需要通過深入探討推移質和懸移質輸運在局部沖刷中的貢獻加以解決。

本文基于水、沙、床面變形相互作用的物理機制，采用瞬態(tài)的URANS方程并輔以RNGk-ε湍流模型封閉方程組、推移質和懸移質輸沙模型以及床面質量守恒方程，分別在清水沖刷和動床沖刷條件下，對水流、泥沙運動及床面變形三者之間的循環(huán)作用求解，模擬了墩柱沖刷過程。通過追蹤沖刷過程中墩周的紊流特性、懸移質和推移質輸沙率的變化，分析懸移質和推移質輸運對局部沖刷發(fā)展的貢獻及其與紊流特性的相關趨勢，以期闡明局部沖刷泥沙輸運模式的選用條件。

2 數(shù)值模型

水流運動由不可壓縮流體的連續(xù)方程和雷諾時均紊流方程（URANS）控制：

式中：xi和ui分別為直角坐標系下的坐標分量和對應的速度分量（i，j，k為坐標分量循環(huán)指標，取值1，2，3，分別表示x，y，z方向）；t為時間；n為水的運動黏度；p為水的平均壓強；fi為流體單位質量力為雷諾應力；Sij為平均應變率張量：

2.1 紊流模型式（2）中紊流的雷諾應力項采用湍流渦黏模型，其本構關系為：

式中：δij為克羅內克符號；為紊動能；νt為紊動渦黏系數(shù)。

式中：k為紊動能；ε為紊動能耗散；Cμ為常數(shù)，通常取值0.09。為了能夠準確捕捉橋墩周圍的漩渦結構，本文采用RNGk-ε模型[8]求解k和ε，該模型在沖刷坑數(shù)值模擬中有較好的效果[20]。

2.2 泥沙運動及河床變形模型在沖刷模擬過程中，需要描述床面泥沙的起動條件，并考慮泥沙的懸移和推移運動兩種可能的運動模式。

2.2.1 泥沙起動 沖刷坑形成及發(fā)展的過程中，局部床面與水流方向非平行，流態(tài)為非均勻流，沖刷坑坡面上泥沙起動的臨界Shields數(shù)與平床時不同，需要對其進行修正。引入量綱一的泥沙粒徑d*：

式中：d為泥沙平均粒徑；ρs為泥沙顆粒密度；g為重力加速度。均勻流床面泥沙顆粒起動的臨界Shields數(shù)θcr如下[21]：

設沖刷坑斜坡面與平床面夾角為β、近床面水流方向與斜坡面上坡方向的夾角為ψ，坡面上泥沙臨界起動Shields數(shù)修正為θ′cr：

式中：φ為泥沙水下休止角，本文取為32°。

2.2.2 泥沙輸運 泥沙的推移質運動主要是求解推移質輸沙率qb，這里采用Meyer-Peter推移質輸沙率公式[22]：

式中：θ為水流的Shields數(shù)；b為經(jīng)驗參數(shù)，其取值與θ有關，本文取8～12。

懸移質泥沙運動由含沙濃度紊動擴散方程表示[23]：

式中：C為懸移質泥沙濃度；βs=1；ωs為泥沙沉降速度，底部床面邊界上的泥沙濃度Ca如下計算[24]：

式中：a=2.5；d為懸移質床面參照高度。懸移質輸沙率由濃度與流速乘積的垂線積分得到。

2.2.3 床面變形方程 采用床面質量守恒方程描述床面變形：

式中：zb為床面高程；n=0.4，為泥沙孔隙率；Ds、Es分別為懸移質泥沙的沉降和卷吸率[23]。沖刷坑演變過程中，還采用沙滑模型修正坡面角度大于泥沙水下休止角的情況[11]。

2.2.4 數(shù)值解及邊界條件 采用有限體積法對上述方程進行空間離散，以求得數(shù)值解。在設置的計算域和坐標系下（如圖1所示），方程求解的邊界條件類型有水流入口、出口邊界、對稱邊界和固壁邊界，圖1中h為水深，D為墩柱直徑。按對稱性處理的邊界包括水表面和兩側面。本文計算時，滿足h/D＞3，屬于窄墩[1]，且水流弗勞德數(shù)Fr較小，此時，自由液面變化對底部床面影響不大[11]，故按對稱邊界處理。兩側面作對稱邊界處理，以消除邊壁對于柱體周圍水流的影響。

圖1 計算域設置

各類邊界條件分別為：（1）在水流入口邊界，給定均勻流x方向速度u和紊動能k垂向分布，其分布由紊流模型在圖1所示計算域無墩的情況下計算得到，其余方向速度為零；（2）在水流的出口邊界上，設各水流速度分量、紊動能和水壓力的法向梯度均為零；（3）在對稱邊界上設置各通量均為零，僅有切向信息；（4）在固壁邊界上，設定無滑移條件，即各速度分量和壓力為零。懸移質的含沙量在底部Ca按式（11）計算，其余邊界上設通量為零。

3 模型驗證

3.1 驗證條件本文數(shù)值模型采用Melville的橋墩沖刷試驗[25]進行驗證。該試驗所用水槽長19 m，寬0.456 m，高0.44 m，直徑D=0.05 m的圓柱墩置于水槽中間，斷面壓縮率B/D=9（其中B為水槽寬度），水槽比降為1×10-4。床面鋪中值粒徑d50=0.385 mm的均勻沙，泥沙休止角為32°。試驗的來流條件為水深h=0.15 m，時均流速U0=0.25 m/s，時均切應力為0.19 Pa，小于泥沙臨界起動切應力0.21 Pa。

針對Melville的試驗，為減少計算量，建立模型時在保證橋墩周圍流態(tài)不受影響的情況下，適當設置了計算域的幾何長度。設坐標系原點位于床面墩中心，墩柱中心距上游入口長度為6D，距下游出口長度為14D，距計算域兩側面均為4.5D，見圖1。床沙厚度設為0.07 m。來流條件與試驗相同。

3.2 網(wǎng)格信息網(wǎng)格劃分質量對于提高墩周渦流及床面變形的模擬精度非常重要。本文采用六面體網(wǎng)格對計算區(qū)域空間進行離散，并在以墩柱軸線為中心，邊長為3D×3D×h的局部范圍內采用加密網(wǎng)格與嵌套網(wǎng)格相結合，如圖2，使墩周流區(qū)的水平向網(wǎng)格（圖2（a））平均尺寸達到0.05D。

本文還特別地參照底部馬蹄渦的形成范圍，對控制床面以上第一層網(wǎng)格以及沙層網(wǎng)格高度的設置進行了較多探索。圓柱繞流中當墩雷諾數(shù)Re＞1×104時，墩柱前馬蹄渦中心的x方向位置為0.64D～0.69D，z方向距離床面為0.04D～0.06D[26-27]，本工況墩雷諾數(shù)為1.27×104，故為了能夠刻畫馬蹄渦，z向第一層網(wǎng)格尺寸設置為0.01D左右，小于渦心距，x-z平面墩周網(wǎng)格的布置見圖2（b）。經(jīng)網(wǎng)格測試計算，最終采用的網(wǎng)格總數(shù)量約為80萬，其中局部加密區(qū)網(wǎng)格數(shù)約為24萬。

圖2 計算網(wǎng)格劃分及嵌套示意圖

為滿足水、沙運動過程瞬態(tài)解的穩(wěn)定和收斂，取時間離散步長Δt始終滿足柯朗數(shù)小于1。

3.3 驗證結果分別針對局部流場和床面變形的結果進行驗證，圖3為t=0時刻z=2 mm的x-y平面墩周流場。對比圖3（a）的流線圖可以發(fā)現(xiàn)，數(shù)值計算和試驗測得的墩前繞流及墩后尾渦尺度都幾乎一致，水流與墩柱面分離點的位置幾近相同。圖3（b）為墩周流速相對于來流速度U0的放大倍數(shù)即U/U0分布，可以看出，數(shù)值計算與試驗測得的分布情況吻合良好，其中最大流速為1.2U0，出現(xiàn)在墩側附近。由于墩柱局部區(qū)域流速增大至大于泥沙臨界起動流速，局部床面開始發(fā)生沖刷，且沖刷坑逐步增大。

圖4為t=30 min時刻墩周沖刷深度等值線的對比，可以看出二者的沖刷形態(tài)很相近，僅在墩后出現(xiàn)差異，這也是很多數(shù)值模擬難以解決的問題[17，28]。該時刻墩周最大沖刷深度出現(xiàn)在墩前，數(shù)值計算給出其大小為4.1 cm，與試驗測量值4.5 cm相對誤差約為8.9%。墩周最大沖刷深度是隨時間變化的變量，記為dt，將30 min時刻的最大沖刷深度記為dt30，則在前30 min內dt/dt30隨時間的變化趨勢如圖5所示?？梢姅?shù)模與試驗結果吻合良好。

圖3 流場驗證

圖4 t=30min沖刷深度等值線圖（單位：cm）

圖5 相對最大沖刷深度隨時間的發(fā)展

圖6 A、B點位置

4 沖刷過程中的懸移質輸運與推移質貢獻

本文數(shù)值模擬考慮了兩種不同的局部沖刷狀態(tài)，即清水沖刷（方案1）和動床沖刷（方案2）。兩種方案的計算域設置均同圖1，床沙平均粒徑均為0.385 mm，方案1的來流與Melville的試驗[25]相同，方案2的來流時均流速U0=0.35 m/s，時均切應力τ0=0.36 Pa，滿足動床沖刷τ0＞τc條件。

在沖刷過程中，紊流場隨著墩柱周圍床面的變化而不斷改變，對泥沙懸移質或推移質的輸運能力也將不斷改變，即泥沙的兩種輸運狀態(tài)對沖刷發(fā)展的貢獻是變化的。為揭示其過程中兩種泥沙輸運的貢獻，在床面和流場變化最強烈的墩側平面位置分別選取A（x=0，y=0.75D）和B（x=0，y=1.00D）點（如圖6所示），并取兩點自水面至床面的垂線，對垂線上的懸移質和推移質輸沙率分別進行計算比較。同一垂線上，推移質輸沙率qb由式（9）計算，推移層以上懸移質輸沙率qs用數(shù)值計算點的含沙量和時均流速之積沿垂線積分獲得，即：

墩周水流和泥沙運動具有強紊動和隨機性，其時均值是采用500個時間步長的瞬態(tài)數(shù)值計算結果進行平均計算而得。為便于分析，本文設分別為無量綱懸移質和推移質輸沙率，并設懸移質輸沙率與推移質輸沙率之比為時均懸推比η，即：

由式（14）可定量得出兩種泥沙輸運對沖刷坑發(fā)展的貢獻。

4.1 清水沖刷條件沖刷坑在初始階段變化很快，通常在發(fā)展至沖刷平衡的前20%時段內沖刷深度約為平衡沖刷深度的80%。現(xiàn)分析沖刷發(fā)展較快的前15 min過程中A、B兩垂線的q*s和q*b變化，如圖7所示。圖7（a）為懸移質輸沙率隨沖刷時間的變化，其中A點因距墩柱近，開始沖刷時輸沙率較大，但隨時間一致性快速減??；B點的輸沙率則呈現(xiàn)先增大后減小的變化。主要是B點在開始沖刷時刻的水流強度小于A點，且沖刷坑尚未發(fā)展至此，但隨著沖刷坑逐漸增大，該處水流強度及輸沙能力增大，當該處床面繼續(xù)沖刷加深，輸沙率則下降。圖7（b）為推移質輸沙率隨沖刷時間的變化，A、B兩點的輸沙率變化趨勢與圖7（a）近乎相同，但量值遠大于后者，且A點的推移質輸沙率減小速率遠小于懸移質輸沙率的減小速率。

圖7 清水沖刷輸沙率隨沖刷時間的變化

清水沖刷的整個過程中，懸移質輸運量占總輸沙量的比例很小。圖8給出了時均懸推比隨沖刷時間的變化，由圖8可見，A點懸推比沖刷開始的瞬間最大，在t=0.5 min時，η≈4%；B點懸推比最大值（約為2.5%）較滯后，出現(xiàn)在t=6 min時。在本文清水沖刷條件下，推移質輸運在沖刷坑發(fā)展過程中始終占控制地位。

圖8 清水沖刷懸推比變化

4.2 動床沖刷條件動床沖刷的沖刷坑較清水沖刷發(fā)展快，t=5min時的沖刷坑的最大深度已大于清水沖刷t=30min時刻的最大深度。仍以A、B兩點為例，揭示沖刷坑發(fā)展過程中的懸移質和推移質輸運的變化，如圖9。與清水沖刷的輸沙率（圖7）相比有以下不同：圖9（a）懸移質輸沙率在A點隨時間單調下降，初期時段的下降率很大，主要是起始時刻輸沙率很高，B點雖也有先增大后減小的趨勢，但增大的時段很短，且初始輸沙率很大；圖9（b）推移質輸沙率隨時間的變化在A點與圖7（b）趨勢相似，但在B點初期隨時間增大的時段很短，輸沙率在沖刷至1.5 min時達到最大，然后則持續(xù)緩慢減小。此工況的懸推比變化過程如圖10所示。圖中A點和B點的懸推比大小和隨時間的變化趨勢幾乎相同，與圖8相比，懸推比增大了1～1.5倍，B點的初始時段變化趨勢與圖8明顯不同。

圖9 動床沖刷輸沙率隨沖刷時間的變化

圖10 動床沖刷懸推比的變化

進一步分析可知，動床沖刷與清水沖刷兩種條件的來流時均切應力之比為1.89，而最大懸移質輸沙率之比則達約10.0（A、B點基本相同），最大推移質輸沙率之比在A點約為2.9，在B點約為2.3?？梢姂乙瀑|和推移質輸沙率增長率大于來流強度增長率，而懸移質輸沙率增長率為來流強度增長率的數(shù)倍。

5 泥沙運動狀態(tài)與水流特征量的相關性

圖11 清水沖刷過程的墩側橫斷面渦量（單位：1/s）

圖12 動床沖刷過程的墩側橫斷面渦量（單位：1/s）

圖13 清水沖刷過程中墩周床面切應力分布（單位：Pa）

圖14 動床沖刷過程中墩周床面切應力分布（單位：Pa）

墩柱周圍具有馬蹄渦流特征，很多文獻針對墩前對稱面進行了研究，本文著重討論沖刷過程中墩側沿柱中心橫剖面的渦量、切應力與泥沙輸運的相關性。

5.1 馬蹄渦作用馬蹄渦的動力特征采用渦量描述，總渦量ω在各方向上的分量ω按如下計算：

在墩側過墩柱中心的橫斷面上，清水沖刷和動床沖刷過程中幾個特征時刻的總渦量ω分布分別如圖11和圖12。圖11（a）和圖12（a）均為t=0（沖刷起始）時刻，清水沖刷和動床沖刷懸移質輸沙率最大的時刻分別為t=6 min和t=1 min。計算垂線平均渦量可知，圖11和圖12中A垂線的平均渦量隨時間衰減，B垂線的平均渦量則由（a）到（b）增大，由（b）到（d）減小，與輸沙率的變化趨勢一致。

5.2 床面切應力沖刷坑發(fā)展過程中水流床面切應力變化見圖13（清水沖刷條件）和圖14（動床沖刷條件），t=0時刻墩柱周圍切應力最大，圖13（a）中τmax=2.4 Pa，圖14（a）中τmax=4.0 Pa，均位于與迎流方向交角約45°處。在隨沖刷坑發(fā)展過程中，最大切應力迅速減小，且其位置移向下游。在清水沖刷的t=0～6 min時段和動床沖刷的t=0～1 min時段床面切應力減小速率大，后續(xù)切應力減小速率變緩，如圖13（b）—（d）和圖14（b）—（d）。分析發(fā)現(xiàn)，墩側面A、B床面切應力與輸沙率變化趨勢一致，與懸推比η具有相關性，見圖15。

墩柱周圍床面切應力較來流增大，與局部馬蹄渦流態(tài)有關。對A、B兩點清水沖刷和動床沖刷兩種條件計算的垂線平均渦量與床面切應力進行分析，得到渦量與切應力關系如圖16所示。

圖15 切應力與懸掛比關系

圖16 渦量與床面切應力關系

5.3 懸移質輸運貢獻與懸浮指數(shù)懸移質輸沙率與懸浮指數(shù)Z有關：

式中：κ=0.4為卡門常數(shù)；u*為摩阻流速。本文根據(jù)沖刷過程中數(shù)值計算的床面切應力，由下式計算相應的摩阻流速：

進而可計算懸浮指數(shù)Z。通過由床面切應力計算的懸浮指數(shù)，可探討墩柱局部沖刷過程中懸移質輸運貢獻與懸浮指數(shù)的關系。設ξ為懸移質輸沙率qs與總輸沙率qt之比，其值表示懸移質輸沙率貢獻，即：

式中qt=qs+qb。取本文兩種計算方案Case1和Case2沖刷過程中A、B點所在垂線計算得到的ξ與Z點繪在圖17中，其中還繪入了Laursen曲線及前人的試驗結果[24]。從圖17可見，本文計算結果與Laursen曲線趨勢很接近，d=0.19 mm的前人試驗數(shù)據(jù)當u*/ωs＞1時與Laursen曲線較符合，見文獻[24]。

圖17 懸浮指數(shù)與懸移質輸沙貢獻比

圖17表明沖刷過程中懸移質輸運的貢獻與懸浮指數(shù)明顯相關。清水沖刷過程中的水流切應力、渦量及摩阻流速均較動床沖刷時小，Z在5.1～6.8之間，數(shù)值相對較大，懸移質輸運貢獻小，ξ＜4%；動床時則因水流動力增大，Z在3.5～6.2之間，數(shù)值相對較小，懸移質輸運貢獻增大，ξ在2%～13%之間。本文計算的條件下，盡管是動床沖刷，其懸移質輸沙率貢獻最大也僅為13%，因為本文的動床條件的來流切應力剛達到泥沙起動臨界值，τ0/τc=1.7，本研究還模擬了條件τ0/τc=11.2局部沖刷，此條件的沖刷坑發(fā)展過程中，懸移質貢獻最大可達45%。因此，無論是清水沖刷還是動床沖刷，懸移質輸運對墩柱局部沖刷坑發(fā)展的貢獻大小與懸浮指數(shù)大小有關。但清水沖刷條件下，懸移質輸運對局部沖刷坑發(fā)展的貢獻很小，數(shù)值模擬時可予以忽略，以節(jié)省計算時間，動床沖刷是否考慮懸移質輸運，可視懸浮指數(shù)決定。

6 結論

本文利用URANS方程、考慮懸移質和推移質輸運的水、沙運動和河床變形方程，以圓柱墩為例進行了結構物局部沖刷數(shù)值模擬，討論了清水沖刷和動床沖刷來流條件下，墩柱局部沖刷坑發(fā)展過程中的局部渦流、床面切應力、懸移質和推移質輸沙率的變化，基于水、沙量變化關系討論了懸移質和推移質輸運在沖刷過程中的貢獻，回答了局部沖刷數(shù)值模擬是否考慮及如何考慮懸移質泥沙輸運的問題，得出以下主要結論：

（1）最大的懸移質和推移質輸沙率都發(fā)生在局部沖刷坑發(fā)展的早期時段，且隨沖刷坑的增大而迅速減小，在沖刷發(fā)展至平衡的后80%時段內，輸沙率緩慢降低。

（2）結構物周圍局部沖刷是強三維和時變問題，沖刷坑的形態(tài)與發(fā)展決定了水流的變化，水流作用決定了床面泥沙運動，各空間點懸移質和推移質輸沙率與水流的床面切應力的變化趨勢存在很強的相關性，墩周局部馬蹄渦導致床面切應力的增大，渦量與床面切應力呈正相關關系。

（3）無論是清水沖刷還是動床沖刷，懸移質輸運對局部沖刷發(fā)展的貢獻大小決定于坑內由床面摩阻流速表示的懸浮指數(shù)，懸移質輸沙率的貢獻，在懸浮指數(shù)大于5的范圍內小于5%，但隨懸浮指數(shù)的減小迅速增大。

（4）清水沖刷條件下，沖刷坑內床面懸浮指數(shù)通常較大，懸移質輸沙率的貢獻率較小，局部沖刷數(shù)值模擬中可不考慮，以節(jié)省計算時間；滿足動床沖刷條件的局部沖刷發(fā)展過程中，懸移質輸沙的貢獻大小與局部床面的懸浮指數(shù)有關。