王 帥，張社榮，戚 藍(lán)，王高輝

(天津大學(xué)水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津300072)

0 引 言

在泄水建筑物及導(dǎo)流洞出口下游設(shè)置防淘墻是保護(hù)岸坡，防止坡腳發(fā)生沖刷的有效措施，但是在水流、波浪等的作用下，防淘墻附近河床仍會(huì)發(fā)生沖刷，河床下切，使得防淘墻暴露于基巖外部，嚴(yán)重時(shí)會(huì)導(dǎo)致防淘墻及上部建筑失穩(wěn)，進(jìn)而造成岸坡坍塌。因此研究河床的沖刷問(wèn)題對(duì)于防淘墻的設(shè)計(jì)具有重要意義。

目前河床沖刷主要集中在兩個(gè)方面:①下游為巖石河床時(shí)，其沖刷過(guò)程大致為:射流作用于巖石河床上，強(qiáng)烈的脈動(dòng)水流進(jìn)入巖石縫隙并沿縫隙傳播，產(chǎn)生的脈動(dòng)壓力促使巖石解體與出穴，從而造成河床沖刷，最有代表性的是劉沛清[1-3]等挑射水流對(duì)巖石河床的沖刷研究；②下游為深覆蓋層河床時(shí)，其沖刷相對(duì)簡(jiǎn)單，多采用模型試驗(yàn)[4]或經(jīng)驗(yàn)公式[5]進(jìn)行研究。然而，考慮到模型試驗(yàn)費(fèi)用高、周期長(zhǎng)的特點(diǎn)，以及針對(duì)防淘墻前河床沖刷經(jīng)驗(yàn)公式的匱乏，使用CFD技術(shù)對(duì)河床沖刷進(jìn)行數(shù)值模擬為本研究提供了一條更加可行的途徑，如鄧軍[6]等通過(guò)建立沖刷坑底部壓力、流速和沖刷深度的平衡關(guān)系模擬了射流對(duì)二維河床的沖刷，并進(jìn)行了實(shí)驗(yàn)驗(yàn)證；陳小莉[7]等以實(shí)際平衡沖坑為邊界條件，通過(guò)CFD技術(shù)模擬了橋臺(tái)局部沖刷，并研究了坑內(nèi)的水流運(yùn)動(dòng)；凌建明[8]等使用CFD模擬了圓柱形橋墩附近三維流場(chǎng)，針對(duì)水流對(duì)河床的沖刷，給出了橋墩附近床面的剪切應(yīng)力分布。以上研究多是集中于橋墩、丁壩等建筑，防淘墻作為一種有效的護(hù)岸工程，其沖刷研究相對(duì)較少。

本文針對(duì)下游河床為深覆蓋層的情況，結(jié)合實(shí)際工程，運(yùn)用CFD技術(shù)，通過(guò)建立導(dǎo)流洞下游河床三維模型，對(duì)防淘墻附近河床的沖刷進(jìn)行數(shù)值模擬，初步探討了導(dǎo)流洞按設(shè)計(jì)洪水導(dǎo)流時(shí)沖刷坑的形成，為防淘墻設(shè)計(jì)提供理論依據(jù)。

1 數(shù)學(xué)模型

1.1 CFD模型

三維水流的控制方程組用不可壓縮流體的連續(xù)方程和動(dòng)量方程來(lái)表示，在笛卡爾坐標(biāo)系下其形式如下:

連續(xù)性方程:

動(dòng)量方程:

k方程:

ε方程:

1.2 沖刷坑的形成及控制

在沖刷過(guò)程中，沖刷坑底部水流流速會(huì)隨著沖刷坑深度的增加而逐漸減小，當(dāng)沖刷坑發(fā)展到穩(wěn)定狀態(tài)時(shí)，意味著沖刷坑底部流速與顆粒臨界起動(dòng)流速之間達(dá)到了一種平衡關(guān)系。其關(guān)系式如下:

式中:u為沖坑底部水流時(shí)均流速；uc為床面顆粒臨界起動(dòng)流速。

根據(jù)王興奎[9]等從床面顆粒受力平衡的觀點(diǎn)推導(dǎo)出起動(dòng)流速的一般表達(dá)形式為:

式中:h為水深；d為顆粒粒徑；A與y分別為待定常數(shù)和指數(shù)，由床面顆粒決定，其中顆粒粒徑采用d50時(shí) ，A=0.146，y=0.586；顆粒粒徑采用d96時(shí)，A=1.006，y=0.167。

將式(6)帶入式(5)得出沖刷坑底部時(shí)均流速與沖刷坑水深之間的平衡關(guān)系式:

式(7)便可作為對(duì)沖刷坑發(fā)展過(guò)程中控制沖坑底部邊界運(yùn)動(dòng)的條件。

定義:

式(8)可用來(lái)作為判斷河床是否遭受沖刷的控制方程。M ＞0說(shuō)明該節(jié)點(diǎn)處水流流速大于顆粒起動(dòng)流速，該處河床發(fā)生沖刷，節(jié)點(diǎn)下移H；反之(M≤0)，該處河床未發(fā)生沖刷，節(jié)點(diǎn)不動(dòng)。當(dāng)所有沖刷坑底部節(jié)點(diǎn)都不滿足下移條件時(shí)，計(jì)算終止，此時(shí)沖刷坑達(dá)到穩(wěn)定狀態(tài)。

節(jié)點(diǎn)下移距離H通過(guò)水下休止角來(lái)確定:H≤Rcot φ，式中:H為節(jié)點(diǎn)下移距離；R為節(jié)點(diǎn)水平間距；φ為床面顆粒水下休止角。

床面顆粒的水下休止角則通過(guò)金臘華[10]等提出的散粒體均勻沙的水下休止角公式來(lái)計(jì)算:

當(dāng)d≤10 mm時(shí)，φ=36.06+4.66lgd；

當(dāng)d＞10 mm時(shí)，φ=40.06+0.97lgd；

1.3 邊界條件及初始條件

模型的進(jìn)口邊界通過(guò)控制上游水位，采用進(jìn)口流量邊界條件來(lái)實(shí)現(xiàn)。

出口邊界采用出流(Outflow)邊界條件，其實(shí)質(zhì)是由orlanski提出[11]的sommerfeld輻射邊界條件，認(rèn)為:

其中:φ為所要輻射的變量，在此為速度U；C為波浪傳播速度；n為輻射邊界的法向向量。

自由表面:水面與空氣接觸表面即為自由表面，采用VOF法[12]來(lái)處理。該方法涉及多相流理論，在每一個(gè)網(wǎng)格用一個(gè)變量F來(lái)標(biāo)志它的狀態(tài)。當(dāng)F=1時(shí)表明網(wǎng)格內(nèi)充滿液體，F(xiàn)=0時(shí)表明網(wǎng)格內(nèi)無(wú)液體存在，當(dāng)0＜F＜1時(shí)表明網(wǎng)格包含自由表面。F函數(shù)可通過(guò)下述方程得到:

固體邊壁采用壁函數(shù)處理:

計(jì)算區(qū)域水流初始流速為斷面平均流速。

1.4 基于CFD的防淘墻河床沖刷分析

FLOW 3D采用有限差分法對(duì)控制方程進(jìn)行離散，壓力速度求解采用GMRES[13]法，時(shí)間差分采用全隱格式。求解過(guò)程為:河床未沖刷時(shí)，水流達(dá)到穩(wěn)定狀態(tài)，獲得計(jì)算區(qū)域各節(jié)點(diǎn)流場(chǎng)參數(shù)，以式(8)判定其是否發(fā)生沖刷，若M＞0，則該節(jié)點(diǎn)下移 H；反之，節(jié)點(diǎn)不動(dòng)。調(diào)整計(jì)算模型，返回第一步重新計(jì)算流場(chǎng)參數(shù)，循環(huán)計(jì)算直到所有節(jié)點(diǎn)滿足河床不沖刷條件，如圖1所示。

圖1 基于CFD的河床沖刷分析流程

2 基于CFD技術(shù)的河床沖刷模擬

2.1 基于CFD技術(shù)的沖刷模擬驗(yàn)證

采用丁壩局部沖刷實(shí)驗(yàn)結(jié)果[14]驗(yàn)證采用CFD模擬沖刷坑的可行性。實(shí)驗(yàn)所用水槽長(zhǎng)10 m，寬0.4 m，模型沙粒徑d50=0.94 mm，全水槽鋪設(shè)6 cm厚，模型丁壩長(zhǎng)7.5 cm，寬 1.8 cm，高 3.0 cm，間距15 cm布設(shè)于水槽右側(cè)邊墻，如圖2所示，水槽流量5.83 L/s，平均水深5.0 cm。

為使沖刷坑達(dá)到穩(wěn)定狀態(tài)，實(shí)驗(yàn)要有足夠長(zhǎng)的水流沖刷時(shí)間。通過(guò)對(duì)水槽中部丁壩周邊河床進(jìn)行測(cè)量，圖3給出了河床沖刷坑順?biāo)飨駻-A(丁壩橫向)剖面、垂直河流向B-B(丁壩軸向)剖面形狀與數(shù)值模擬計(jì)算剖面形狀對(duì)比圖。

圖2 丁壩局部沖刷三維示意圖

圖3 沖刷坑剖面形狀對(duì)比圖

從圖3中可以看出，實(shí)驗(yàn)與數(shù)值模擬得出的河床沖刷坑均位于丁壩端偏向上游。計(jì)算所得沖深最大值與實(shí)驗(yàn)值相比略小，計(jì)算值為9.54 mm，實(shí)驗(yàn)值為11.0 mm，相對(duì)誤差為13.3%。沖刷坑位置與實(shí)驗(yàn)結(jié)果相同，最大沖深計(jì)算值與實(shí)驗(yàn)值相近，沖刷坑性態(tài)相似但存在一定差異，一方面是由于床面顆粒起動(dòng)具有一定的隨機(jī)性；另一方面，在相同流速情況下，顆粒順坡起動(dòng)要較水平床面起動(dòng)容易，逆坡起動(dòng)則相對(duì)較難些；同時(shí)數(shù)值模擬是建立在顆粒粒徑相對(duì)均勻的基礎(chǔ)之上的，與模型試驗(yàn)采用的非均勻沙有一定出入。但就沖坑位置、沖刷深度以及沖坑形成規(guī)律而言，采用CFD技術(shù)模擬深覆蓋層河床沖刷具有一定的可行性。

2.2 基于CFD技術(shù)的防淘墻河床沖刷模擬

將該法應(yīng)用于某電站導(dǎo)流洞下游防淘墻附近河床的沖刷中?；举Y料如下:導(dǎo)流洞出口下游左岸布置一道混凝土擋墻。在擋墻基礎(chǔ)中設(shè)一道混凝土防淘墻(見(jiàn)圖4)，長(zhǎng)140 m，深12m，厚3m，防淘墻底高程1 293.0m，頂高程1305.0 m。該段河床平均高程1 305.0 m，下游河道自然順直，左岸邊坡自然坡度40°～50°。河床覆蓋層由沖積漂卵石組成，顆粒粒徑d96=10 mm，厚度3.0 m～10.0 m，河道中央覆蓋層較厚(見(jiàn)圖5)。

圖4 防淘墻平面布置圖

圖5 防淘墻橫剖面圖

建立防淘墻區(qū)域河床三維模型，如圖6所示，模型選取導(dǎo)流洞出口到防淘墻下游125 m處，模型全長(zhǎng)444.5 m，寬472.0 m，高237.5 m，按全年 P=5%的設(shè)計(jì)洪水標(biāo)準(zhǔn)，相應(yīng)導(dǎo)流流量7 180 m3/s，下游水位1 317.6 m，平均水深12.6 m。

計(jì)算網(wǎng)格為470×410×250，如圖7所示，X方向上網(wǎng)格長(zhǎng)度為0.4 m～3 m，Y方向上網(wǎng)格長(zhǎng)度為0.4 m～3 m，Z方向上網(wǎng)格長(zhǎng)度為0.2 m～1 m。采用重點(diǎn)區(qū)域局部加密，邊緣區(qū)域適當(dāng)放寬的原則。

圖6 防淘墻區(qū)域三維模型

圖7 計(jì)算網(wǎng)格剖分

追蹤床面部分關(guān)鍵點(diǎn)M值隨沖坑發(fā)展的變化過(guò)程，如圖8所示。

圖8 沿C D線剖面關(guān)鍵點(diǎn)的M值隨沖坑發(fā)展的變化過(guò)程

圖8(a)給出了沿CD線(見(jiàn)圖4)防淘墻橫剖面中關(guān)鍵點(diǎn)分布情況，其中c′點(diǎn)為床面c點(diǎn)在沖刷坑底部的投影(a，b，d，e同)，圖8(b)為 c點(diǎn)向c′點(diǎn)移動(dòng)過(guò)程中M值的變化曲線(a，b，d，e同)。由圖可知，河床未沖刷時(shí)，各點(diǎn)的M值都較大，當(dāng)沖深1 m后，M大幅下降，隨著沖坑深度的不斷加深，M不斷減小，并最終等于0，此時(shí)沖刷坑達(dá)到穩(wěn)定狀態(tài)，距離防淘墻1.5 m處的c點(diǎn)沖坑最深。

圖9給出了沿防淘墻縱向AB線(見(jiàn)圖4)剖面中 J(h 、i、k 、l)點(diǎn)向 J′(h′，i′，k′，l′)運(yùn)動(dòng)過(guò)程中 M值的變化曲線。從圖中可以看出，各點(diǎn)的M值在河床未沖刷時(shí)較大，河床沖深1 m后，M值大幅下降，并隨著沖坑深度的不斷增加而減小，最終等于0，沖坑達(dá)到穩(wěn)定狀態(tài)，在防淘墻0+004.0處的j點(diǎn)沖坑最深。

圖10給出了防淘墻河床在不同深度下的沖坑形態(tài)?？梢钥闯?，防淘墻附近河床的沖刷主要集中于防淘墻的上游河床，在導(dǎo)流洞出口、防淘墻起始部位的河床沖刷較為嚴(yán)重，形成一個(gè)形似椎體的沖刷坑，而下游河床基本未受到?jīng)_刷。

圖9 沿AB線剖面關(guān)鍵點(diǎn)的M值隨沖坑發(fā)展的變化過(guò)程

圖10 沖刷坑發(fā)展過(guò)程

圖11給出了沿 AB線(見(jiàn)圖4)、CD線(見(jiàn)圖4)沖刷坑剖面的最終形態(tài)。在圖11(a)中，沿 AB線剖面距防淘墻1.5 m，原點(diǎn)處即為防淘墻起始端，水流對(duì)防淘墻上游部位河床產(chǎn)生沖刷，沖刷部位主要集中于防0+000.0～防0+050.0段，其中防0+000.0～防0+005.0段沖刷較為嚴(yán)重，沖深最大值為4.6 m，位于防0+004.0處。圖11(b)中，沿 CD 剖面位于樁號(hào)防0+0.004.0處，可以看出，距離防淘墻12 m范圍內(nèi)的河床發(fā)生沖刷，沖刷最深處距防淘墻1.5 m。防淘墻附近的河床沖刷僅限于覆蓋層，并未對(duì)河床基巖造成沖刷。

3 結(jié) 論

本文依托工程實(shí)例，基于CFD技術(shù)，從流速角度對(duì)深覆蓋層河床沖刷進(jìn)行了研究，為河床沖刷的數(shù)值模擬提供了一種新的思路。研究結(jié)果表明:

(1)以沖坑底部流速與臨界起動(dòng)流速之間的平衡關(guān)系作為判斷河床發(fā)生沖刷的控制條件，基于CFD技術(shù)，對(duì)深覆蓋層河床進(jìn)行了數(shù)值模擬。與文獻(xiàn)[14]中丁壩局部沖刷實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比，采用CFD技術(shù)模擬深覆蓋層河床沖刷具有一定的可行性。

圖11 沖刷坑剖面

(2)對(duì)某電站導(dǎo)流洞下游防淘墻河床沖刷進(jìn)行了數(shù)值模擬，導(dǎo)流洞按全年P(guān)=5%的設(shè)計(jì)洪水導(dǎo)流時(shí)，防淘墻樁號(hào)防0+000.0～防0+050.0段河床發(fā)生明顯沖刷，其中在防0+000.0～防0+005.0段河床沖刷尤為嚴(yán)重，在防0+004.0、距防淘墻1.5 m處河床形成一個(gè)形似圓錐體的沖坑，沖坑最深處達(dá)4.6 m，該處覆蓋層厚4.9m，水流并未沖刷到河床基巖。

(3)通過(guò)數(shù)值計(jì)算結(jié)果與模型實(shí)驗(yàn)進(jìn)行對(duì)比，可以看出，就目前技術(shù)水平而言，數(shù)值計(jì)算在某些方面可以替代模型試驗(yàn)。首次將數(shù)值模擬手段應(yīng)用于防淘墻設(shè)計(jì)中，為防淘墻埋深提供設(shè)計(jì)依據(jù)，對(duì)降低防淘墻的設(shè)計(jì)周期與成本具有重要意義。

[1]李愛(ài)華，劉沛清.脈動(dòng)壓力在消力池底板縫隙傳播的瞬變流模型和滲流模型統(tǒng)一性探討[J].水利學(xué)報(bào)，2005，36(10):1236-1240.

[2]李愛(ài)華，劉沛清.巖石河床在沖擊水流脈動(dòng)壓力作用下的解體破壞機(jī)理[J].水利學(xué)報(bào)，2007，38(11):1324-1328.

[3]李愛(ài)華，劉沛清，許唯臨.復(fù)雜裂隙網(wǎng)絡(luò)內(nèi)沖擊射流脈動(dòng)壓力傳播過(guò)程數(shù)值研究[J].四川大學(xué)學(xué)報(bào)，2009，41(5):36-41.

[4]沈 波，艾翠玲，田偉平，等.小橋涵自由式出流沖刷試驗(yàn)[J].長(zhǎng)安大學(xué)學(xué)報(bào)，2007，27(3):61-66.

[5]詹義正，王 軍，談廣鳴，等.橋墩局部沖刷的試驗(yàn)研究[J].武漢大學(xué)學(xué)報(bào)，2006 ，39(5):1-4，9.

[6]鄧 軍，許唯臨，楊忠超，等.基巖沖刷的數(shù)值模擬[J].水科學(xué)進(jìn)展 ，2005，16(1):47-51.

[7]陳小莉，馬吉明.橋臺(tái)局部沖刷坑內(nèi)水流運(yùn)動(dòng)的三維數(shù)值模擬[J].水動(dòng)力學(xué)研究與進(jìn)展，2007，22(6):689-695.

[8]凌建明，林小平，趙鴻鐸.圓柱形橋墩附近三維流場(chǎng)及河床局部沖刷分析[J].同濟(jì)大學(xué)學(xué)報(bào)，2007，35(5):582-586.

[9]王興奎，陳稚聰，張 仁.長(zhǎng)江寸灘站卵石推移質(zhì)的運(yùn)動(dòng)規(guī)律[J].水利學(xué)報(bào)，1992，(4):32-38.

[10]金臘華，石秀清.試論模型沙的水下休止角[J].泥沙研究，1990，(3):87-93.

[11]Orlanski I.A simple boundary condition for unbounded hyperbolic flows[J].Journal of Computational Physics，1976，21:251-269.

[12]王永學(xué).VOF方法數(shù)模直墻式建筑物前的波浪破碎過(guò)程[J].自然科學(xué)進(jìn)展:國(guó)家重點(diǎn)實(shí)驗(yàn)室通訊，1993，3(6):553-559.

[13]Saad Y，SchultzM H.GMRES:A generalized minimal residual algorithm for solving nonsymmetric linear systerms[J].SIAM Journal on Scientific and Statistical Computing，1986，7:856-869.

[14]彭 靜，玉井信行，何原能久.丁壩壩頭沖淤的三維數(shù)值模擬[J].泥沙研究，2002，(1):25-29.