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

王 帥，張社榮，戚 藍，王高輝

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

0 引 言

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

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

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

1 數學模型

1.1 CFD模型

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

連續性方程:

動量方程:

k方程:

ε方程:

1.2 沖刷坑的形成及控制

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

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

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

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

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

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

定義:

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

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

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

當d≤10 mm時，φ=36.06+4.66lgd；

當d＞10 mm時，φ=40.06+0.97lgd；

1.3 邊界條件及初始條件

模型的進口邊界通過控制上游水位，采用進口流量邊界條件來實現。

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

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

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

固體邊壁采用壁函數處理:

計算區域水流初始流速為斷面平均流速。

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

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

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

2 基于CFD技術的河床沖刷模擬

2.1 基于CFD技術的沖刷模擬驗證

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

為使沖刷坑達到穩定狀態，實驗要有足夠長的水流沖刷時間。通過對水槽中部丁壩周邊河床進行測量，圖3給出了河床沖刷坑順水流向A-A(丁壩橫向)剖面、垂直河流向B-B(丁壩軸向)剖面形狀與數值模擬計算剖面形狀對比圖。

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

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

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

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

將該法應用于某電站導流洞下游防淘墻附近河床的沖刷中?；举Y料如下:導流洞出口下游左岸布置一道混凝土擋墻。在擋墻基礎中設一道混凝土防淘墻(見圖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 防淘墻橫剖面圖

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

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

圖6 防淘墻區域三維模型

圖7 計算網格剖分

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

圖8 沿C D線剖面關鍵點的M值隨沖坑發展的變化過程

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

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

圖10給出了防淘墻河床在不同深度下的沖坑形態?？梢钥闯?，防淘墻附近河床的沖刷主要集中于防淘墻的上游河床，在導流洞出口、防淘墻起始部位的河床沖刷較為嚴重，形成一個形似椎體的沖刷坑，而下游河床基本未受到沖刷。

圖9 沿AB線剖面關鍵點的M值隨沖坑發展的變化過程

圖10 沖刷坑發展過程

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

3 結 論

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

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

圖11 沖刷坑剖面

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

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

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

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

[3]李愛華，劉沛清，許唯臨.復雜裂隙網絡內沖擊射流脈動壓力傳播過程數值研究[J].四川大學學報，2009，41(5):36-41.

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

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

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

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

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

[9]王興奎，陳稚聰，張 仁.長江寸灘站卵石推移質的運動規律[J].水利學報，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]王永學.VOF方法數模直墻式建筑物前的波浪破碎過程[J].自然科學進展:國家重點實驗室通訊，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]彭 靜，玉井信行，何原能久.丁壩壩頭沖淤的三維數值模擬[J].泥沙研究，2002，(1):25-29.