基于CFD的防淘墻河床沖刷模擬方法

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

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

0 引 言

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

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

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

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

1.1 CFD模型

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

連續(xù)性方程:

動量方程:

k方程:

ε方程:

1.2 沖刷坑的形成及控制

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

式中:u為沖坑底部水流時均流速；uc為床面顆粒臨界起動流速。

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

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

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

式(7)便可作為對沖刷坑發(fā)展過程中控制沖坑底部邊界運動的條件。

定義:

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

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

床面顆粒的水下休止角則通過金臘華[10]等提出的散粒體均勻沙的水下休止角公式來計算:

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

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

1.3 邊界條件及初始條件

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

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

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

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

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

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

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

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

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

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

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

采用丁壩局部沖刷實驗結(jié)果[14]驗證采用CFD模擬沖刷坑的可行性。實驗所用水槽長10 m，寬0.4 m，模型沙粒徑d50=0.94 mm，全水槽鋪設(shè)6 cm厚，模型丁壩長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)，實驗要有足夠長的水流沖刷時間。通過對水槽中部丁壩周邊河床進(jìn)行測量，圖3給出了河床沖刷坑順?biāo)飨駻-A(丁壩橫向)剖面、垂直河流向B-B(丁壩軸向)剖面形狀與數(shù)值模擬計算剖面形狀對比圖。

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

圖3 沖刷坑剖面形狀對比圖

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

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

將該法應(yīng)用于某電站導(dǎo)流洞下游防淘墻附近河床的沖刷中?；举Y料如下:導(dǎo)流洞出口下游左岸布置一道混凝土擋墻。在擋墻基礎(chǔ)中設(shè)一道混凝土防淘墻(見圖4)，長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，河道中央覆蓋層較厚(見圖5)。

圖4 防淘墻平面布置圖

圖5 防淘墻橫剖面圖

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

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

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

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

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

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

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

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

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

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

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

圖11給出了沿 AB線(見圖4)、CD線(見圖4)沖刷坑剖面的最終形態(tài)。在圖11(a)中，沿 AB線剖面距防淘墻1.5 m，原點處即為防淘墻起始端，水流對防淘墻上游部位河床產(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 剖面位于樁號防0+0.004.0處，可以看出，距離防淘墻12 m范圍內(nèi)的河床發(fā)生沖刷，沖刷最深處距防淘墻1.5 m。防淘墻附近的河床沖刷僅限于覆蓋層，并未對河床基巖造成沖刷。

3 結(jié) 論

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

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

圖11 沖刷坑剖面

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

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

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

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

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

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

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

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

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

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

[9]王興奎，陳稚聰，張 仁.長江寸灘站卵石推移質(zhì)的運動規(guī)律[J].水利學(xué)報，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ù)模直墻式建筑物前的波浪破碎過程[J].自然科學(xué)進(jì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.