齊梅蘭，石粕辰

（1.北京交通大學(xué) 土木建筑工程學(xué)院，北京 100044；2.北京交通大學(xué) 結(jié)構(gòu)風(fēng)工程與城市風(fēng)環(huán)境北京市重點(diǎn)實(shí)驗(yàn)室，北京 100044）

1 研究背景

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

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

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

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

2 數(shù)值模型

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

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

2.1 紊流模型式（2）中紊流的雷諾應(yīng)力項(xiàng)采用湍流渦黏模型，其本構(gòu)關(guān)系為：

式中：δij為克羅內(nèi)克符號(hào)；為紊動(dòng)能；νt為紊動(dòng)渦黏系數(shù)。

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

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

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

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

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

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

2.2.2 泥沙輸運(yùn) 泥沙的推移質(zhì)運(yùn)動(dòng)主要是求解推移質(zhì)輸沙率qb，這里采用Meyer-Peter推移質(zhì)輸沙率公式[22]：

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

懸移質(zhì)泥沙運(yùn)動(dòng)由含沙濃度紊動(dòng)擴(kuò)散方程表示[23]：

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

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

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

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

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

圖1 計(jì)算域設(shè)置

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

3 模型驗(yàn)證

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

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

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

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

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

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

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

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

圖3 流場驗(yàn)證

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

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

圖6 A、B點(diǎn)位置

4 沖刷過程中的懸移質(zhì)輸運(yùn)與推移質(zhì)貢獻(xiàn)

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

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

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

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

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

圖7 清水沖刷輸沙率隨沖刷時(shí)間的變化

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

圖8 清水沖刷懸推比變化

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

圖9 動(dòng)床沖刷輸沙率隨沖刷時(shí)間的變化

圖10 動(dòng)床沖刷懸推比的變化

進(jìn)一步分析可知，動(dòng)床沖刷與清水沖刷兩種條件的來流時(shí)均切應(yīng)力之比為1.89，而最大懸移質(zhì)輸沙率之比則達(dá)約10.0（A、B點(diǎn)基本相同），最大推移質(zhì)輸沙率之比在A點(diǎn)約為2.9，在B點(diǎn)約為2.3。可見懸移質(zhì)和推移質(zhì)輸沙率增長率大于來流強(qiáng)度增長率，而懸移質(zhì)輸沙率增長率為來流強(qiáng)度增長率的數(shù)倍。

5 泥沙運(yùn)動(dòng)狀態(tài)與水流特征量的相關(guān)性

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

圖12 動(dòng)床沖刷過程的墩側(cè)橫斷面渦量（單位：1/s）

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

圖14 動(dòng)床沖刷過程中墩周床面切應(yīng)力分布（單位：Pa）

墩柱周圍具有馬蹄渦流特征，很多文獻(xiàn)針對(duì)墩前對(duì)稱面進(jìn)行了研究，本文著重討論沖刷過程中墩側(cè)沿柱中心橫剖面的渦量、切應(yīng)力與泥沙輸運(yùn)的相關(guān)性。

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

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

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

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

圖15 切應(yīng)力與懸掛比關(guān)系

圖16 渦量與床面切應(yīng)力關(guān)系

5.3 懸移質(zhì)輸運(yùn)貢獻(xiàn)與懸浮指數(shù)懸移質(zhì)輸沙率與懸浮指數(shù)Z有關(guān)：

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

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

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

圖17 懸浮指數(shù)與懸移質(zhì)輸沙貢獻(xiàn)比

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

6 結(jié)論

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

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

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

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

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