淹沒堆積體附近水流特性研究

李馨楠，向 波，張 婧，王 東，趙紅兮

（1.西華大學(xué) 流體及動(dòng)力機(jī)械教育部重點(diǎn)實(shí)驗(yàn)室，四川 成都 610039；2.西華大學(xué) 流體機(jī)械及工程四川省重點(diǎn)實(shí)驗(yàn)室，四川 成都 610039；3.四川省公路規(guī)劃勘察設(shè)計(jì)研究院有限公司，四川 成都 610041）

有關(guān)研究團(tuán)隊(duì)已經(jīng)發(fā)表一系列非淹沒條件下堆積體的研究成果［3－5］，但是對(duì)淹沒堆積體的相關(guān)水沙研究尚屬空白。繞流結(jié)構(gòu)體附近水流紊動(dòng)比較強(qiáng)烈，尤其是回流區(qū)旋渦、氣泡較多，給試驗(yàn)的精密測量帶來困難，因此許多學(xué)者采用數(shù)值模擬技術(shù)對(duì)淹沒丁壩的水流特性展開研究。堆積體和丁壩均屬于繞流結(jié)構(gòu)體中的一種，只是形態(tài)有所不同。與丁壩相似，堆積體存在于河道中，會(huì)影響其附近的水流特性，堆積體的存在會(huì)使上游的壅水現(xiàn)象更加明顯，沿程水深分布將受到流量、堆積體形態(tài)、底坡坡度等因素的影響，呈現(xiàn)較為明顯的規(guī)律。蔡喆偉等［6］為研究新型結(jié)構(gòu)淹沒丁壩的水流特性，基于有限體積法與自由液面捕捉法，通過流體計(jì)算軟件Fluent建立三維數(shù)學(xué)模型對(duì)不同流速、水深與壩長條件下新型結(jié)構(gòu)齒形丁壩的水流特性進(jìn)行研究。寧健等［7］通過建立數(shù)學(xué)模型，模擬在水流作用下的泥沙運(yùn)動(dòng)，預(yù)測了泥沙侵蝕、對(duì)流擴(kuò)散以及沉積等過程，同時(shí)考慮河床坡度對(duì)泥沙輸移以及臨界希爾茲數(shù)的影響，得到丁壩附近流場、床面變形以及沖刷坑形態(tài)等，與試驗(yàn)結(jié)果吻合較好。陳宇豪等［8］為研究勾頭丁壩的水流特性，選用Flow－3D中的RNGk－ε紊流模型建立了勾頭丁壩水流計(jì)算三維數(shù)學(xué)模型，并采用淹沒式丁壩的物理模型試驗(yàn)驗(yàn)證了所建模型能較好地模擬過流丁壩水流。王小明等［9］基于OpenFOAM開源程序包模擬了不同透空率潛壩的三維流場，研究分析了透水潛壩周圍三維水流特性受透空率和流量的影響規(guī)律。

堆積體的存在會(huì)對(duì)下游水電站的安全運(yùn)行造成威脅，而其工程處置的首要問題應(yīng)是明晰堆積體附近的水流特性，淹沒丁壩的研究思路可為本研究提供參考。筆者采用數(shù)值模擬技術(shù)對(duì)淹沒堆積體附近水流的水力特性進(jìn)行研究，以期為后續(xù)關(guān)于堆積體泥沙起動(dòng)及泥沙輸移的試驗(yàn)研究打下基礎(chǔ)，并為堆積體影響下峽谷地區(qū)河道治理方案的優(yōu)化提供理論支持。

1 數(shù)值模型及其驗(yàn)證

1.1 控制方程

本研究采用標(biāo)準(zhǔn)k－ε湍流模型，其水流連續(xù)性方程為

動(dòng)量守恒方程為

湍流模型控制方程為

式中：ρ為水流密度，kg／m3；t為時(shí)間，s；xi和xj分別代表笛卡兒坐標(biāo)系的x、y軸；ui和uj分別為i、j方向上的流速，m／s；μ為液體動(dòng)力黏滯系數(shù)；Si為動(dòng)量守恒方程源項(xiàng)；P為壓力，N；k為湍動(dòng)能；ε為湍動(dòng)能耗散率；μt為湍動(dòng)黏度；Gk為由平均速度梯度得到的湍動(dòng)能的產(chǎn)生項(xiàng)；Gb為浮力引起的湍動(dòng)能的產(chǎn)生項(xiàng)；YM為可壓湍流中脈動(dòng)擴(kuò)張的作用；G1ε、G2ε、G3ε為經(jīng)驗(yàn)常數(shù)；σk、σε為與k和ε對(duì)應(yīng)的Prandtl數(shù)；Sk、Sε為用戶定義的源項(xiàng)。

1.2 網(wǎng)格劃分及邊界條件

根據(jù)試驗(yàn)水槽尺寸，在ICEM中建模和劃分網(wǎng)格，圖1為模型的三維視圖。堆積體由上、下游對(duì)稱面和側(cè)面構(gòu)成，上、下游面均采用半圓錐體結(jié)構(gòu)，其平面投影為扇形，堆積體模型側(cè)向坡度S＝45°，沿水流方向上的最大長度l＝50 cm，橫向?qū)挾萣＝25 cm，高度Hcon＝10 cm，進(jìn)占比為b／B＝0.5，采用非結(jié)構(gòu)性四面體網(wǎng)格進(jìn)行網(wǎng)格劃分，將堆積體模型周圍局部加密，網(wǎng)格總數(shù)為85萬，網(wǎng)格劃分結(jié)果如圖2、圖3所示。本文采用速度和壓力進(jìn)口邊界，出口邊界選擇壓力出口邊界，采用VOF（Volume of Fluid）法計(jì)算自由水面。對(duì)于水氣兩相流，網(wǎng)格內(nèi)水流所占比例采用體積分?jǐn)?shù)α表示，在每個(gè)單元體中，水和空氣的體積分?jǐn)?shù)之和為1。速度進(jìn)口由進(jìn)口面的水相進(jìn)口表示，壓力進(jìn)口由進(jìn)口面的氣相進(jìn)口表示（設(shè)為恒壓）。本文采用SIMPLE算法進(jìn)行求解，即水的體積分?jǐn)?shù)α可由耦合格式計(jì)算。

9月25日下午，中歐水資源交流平臺(tái)地下水管理專題會(huì)議召開。水利部副部長胡四一、丹麥環(huán)境大臣奧肯出席開幕式并致辭。

圖1 三維模型

圖2 計(jì)算區(qū)域網(wǎng)格劃分結(jié)果

圖3 堆積體附近局部網(wǎng)格加密情況

1.3 模型驗(yàn)證

驗(yàn)證試驗(yàn)在西華大學(xué)流體及動(dòng)力機(jī)械教育部重點(diǎn)實(shí)驗(yàn)室進(jìn)行。試驗(yàn)的供水系統(tǒng)為循環(huán)式供水系統(tǒng)，矩形水槽長L＝18 m，寬B＝0.5 m，高H＝0.42 m，水槽底部坡降為0.15%，水深用LH－1自動(dòng)水位儀測量，流速采用螺旋測速儀測量。

試驗(yàn)流量Q＝35.8 L／s，測量布置如圖4所示，x、y、z分別為橫向、縱向、垂向坐標(biāo)。研究采用以下無量綱化的處理方式：x?＝x／l和y?＝y(tǒng)／l分別表示斷面和測線的相對(duì)位置，以z?＝z／h0代表z方向的相對(duì)位置，u?＝u／u0代表相對(duì)流速。

圖4 水槽試驗(yàn)測量布置

選取y?＝0.3、y?＝0.5、y?＝0.7三條測線水深的橫向變化進(jìn)行對(duì)比，實(shí)測值與模擬值如圖5所示，可以看出兩者較為吻合，平均誤差為2.88%，堆積體下游回流區(qū)附近的誤差相對(duì)較大。選取了斷面x?＝－0.5、x?＝0及x?＝0.5的中心測線（y?＝0. 5）的縱向流速的垂向變化進(jìn)行對(duì)比分析，結(jié)果如圖6所示，平均誤差為1.13%，同樣，堆積體下游回流區(qū)附近的誤差相對(duì)較大。其原因是在回流區(qū)伴隨大量氣泡、渦旋的產(chǎn)生，導(dǎo)致水位和流速的測量結(jié)果誤差相對(duì)偏大?？傮w來說，模擬值與實(shí)測值整體較為吻合，水深的橫向變化和縱向流速的沿程變化規(guī)律與水槽試驗(yàn)結(jié)果一致。本文淹沒堆積體數(shù)值模擬研究的變量包括流量、進(jìn)占比和坡度，計(jì)算工況見表1。

表1 模擬計(jì)算工況

圖5 縱向水面線驗(yàn)證結(jié)果

圖6 縱向流速驗(yàn)證結(jié)果

參考前人關(guān)于淹沒丁壩的研究，筆者引入δ（淹沒系數(shù)）作為堆積體處于淹沒狀態(tài)時(shí)衡量其淹沒程度的指標(biāo)，定義如下：

式中：h為對(duì)應(yīng)工況下堆積體頂部中軸線x?＝0斷面的水深，m；h0為對(duì)應(yīng)工況下的均勻流水深，m。

2 水面線分析

從3種淹沒系數(shù)（小、中、大）中選擇具有代表性的工況進(jìn)行分析，其中流量因素以D1、D2、D3為例，進(jìn)占比因素以M1、M2、D3為例，側(cè)向坡度因素以N1、N2、N3為例，以所有測線的平均水深代表各橫斷面的水深，計(jì)算不同工況下沿程相對(duì)水深變化情況，如圖7所示。

圖7 斷面平均水深沿程分布

整體上看，3種不同因素影響下的水面線沿程變化趨勢較為一致。在堆積體上游段，水面線因受到堆積體的阻水作用而壅高，水面線較為平穩(wěn)。進(jìn)入堆積體漸縮段，水面線下降較為迅速，水流的慣性作用使水深在x?＝0.3處最小。水流經(jīng)過堆積體中軸線后，過流斷面逐漸恢復(fù)，水面線開始緩慢回升。堆積體對(duì)水流的影響越往下游越小，水深逐漸恢復(fù)到均勻流水深。由此可以看出，處于淹沒狀態(tài)時(shí)，堆積體附近水面線沿程變化的一般規(guī)律與非淹沒狀態(tài)下無明顯差異［10］。比較分析不同因素影響下的水面線可以發(fā)現(xiàn)，上游水面線的壅水高度與流量成反比，與進(jìn)占比和坡度成正比。從水面線變化幅度來看，堆積體進(jìn)占比、側(cè)向坡度、流量3種因素的影響程度依次加強(qiáng)。與非淹沒情況相比，堆積體側(cè)向坡度因素的阻流影響較強(qiáng)，而進(jìn)占比的影響較弱。

相對(duì)壅水水深h′以x?＝－4.5～－0.5各斷面水深h與均勻流水深h0的比值表示，即h′＝h／h0，h′越大，則上游壅水越明顯。相對(duì)最小水深h′min以斷面x?＝0.3的平均水深與均勻流水深的比值表示，h′min越小，則堆積體段的水深下降幅度越大。將各工況條件下的淹沒系數(shù)與相對(duì)壅水水深及相對(duì)最小水深進(jìn)行數(shù)據(jù)擬合，結(jié)果表明淹沒系數(shù)與相對(duì)壅水水深及相對(duì)最小水深具有良好的指數(shù)關(guān)系：

由圖8、圖9可以看出，淹沒系數(shù)越大，相對(duì)壅水水深越小，相對(duì)最小水深越大，且淹沒系數(shù)對(duì)相對(duì)壅水水深的影響隨著淹沒系數(shù)的增大而減弱。由此可知，堆積體處于淹沒狀態(tài)時(shí)，對(duì)堆積體頂部以下的水流仍然起到調(diào)節(jié)作用，其調(diào)節(jié)作用隨著淹沒系數(shù)的增大而減小，這與丁壩處于淹沒狀態(tài)時(shí)的調(diào)節(jié)機(jī)理相似［11］。當(dāng)淹沒系數(shù)足夠大時(shí)，可以認(rèn)為水面線幾乎不受堆積體影響。

圖8 淹沒系數(shù)與相對(duì)壅水水深關(guān)系曲線

圖9 淹沒系數(shù)與相對(duì)最小水深關(guān)系曲線

3 流場特性分析

3.1 紊動(dòng)特性分析

水流的紊動(dòng)強(qiáng)度對(duì)河道水流挾沙能力的影響較大，本文用紊動(dòng)能來表征水流整體的紊動(dòng)情況，其計(jì)算公式為

式中：TKE為紊動(dòng)能；u、v、w分別為橫向、縱向、豎向的脈動(dòng)流速。

以工況D2為例，分析不同流層的紊動(dòng)能分布情況，選取底層（z?＝0.02）、中層（z?＝0.51）和上層（z?＝0.75）為代表，分析不同流層的紊動(dòng)能分布規(guī)律，如圖10所示。整體上看，堆積體上游壅水區(qū)的紊動(dòng)能數(shù)值均較小，此區(qū)域容易造成泥沙淤積。底部過流斷面進(jìn)占程度較高，水流在近底層受到明顯壓縮，水流擴(kuò)散導(dǎo)致的高紊動(dòng)區(qū)出現(xiàn)在底部，且隨流層升高，水流的擴(kuò)散作用減弱，紊動(dòng)強(qiáng)度也隨之減弱。下游出現(xiàn)的高紊動(dòng)區(qū)主要是主流擴(kuò)散導(dǎo)致的，剪切流所引起的高紊動(dòng)雖然存在，但和水流擴(kuò)散的高紊動(dòng)連成一片，之后隨著流層升高，回流現(xiàn)象開始出現(xiàn)［12］，剪切流所引起的高紊動(dòng)才更加清晰地被觀察到。隨著回流現(xiàn)象的出現(xiàn)，高紊動(dòng)區(qū)域的分布和強(qiáng)度開始受到回流現(xiàn)象的主導(dǎo)，同時(shí)，堆積體束窄作用減弱，水流擴(kuò)散導(dǎo)致的高紊動(dòng)已不再明顯，高紊動(dòng)區(qū)已逐漸轉(zhuǎn)移到回流區(qū)附近，隨著回流影響范圍的增大，其紊動(dòng)能也逐漸增大，越往下游紊動(dòng)強(qiáng)度越弱。由此可知，高紊動(dòng)區(qū)受到回流的影響會(huì)發(fā)生改變，流層越往上，受回流的影響越大。泥沙運(yùn)動(dòng)與水流紊動(dòng)強(qiáng)度有密切的關(guān)系，泥沙發(fā)生運(yùn)動(dòng)，河床剪切力并不是唯一的影響因素，而紊動(dòng)引起的上舉力會(huì)促使水流進(jìn)一步挾帶和運(yùn)移河床表面的泥沙顆粒［13］，因此在回流區(qū)更易發(fā)生泥沙輸移，且水流挾帶泥沙的能力更強(qiáng)。

圖10 不同水深紊動(dòng)能分布

3.2 渦量分析

渦系的存在影響了水流中物質(zhì)的輸移特性且影響水流的能量耗散過程，研究堆積體附近的渦量場使得對(duì)繞流結(jié)構(gòu)體作用下的沖淤特性具有更加清晰的認(rèn)識(shí)。根據(jù)文獻(xiàn)［14］推導(dǎo)的渦量表達(dá)式，在豎直方向上渦量采用如下公式計(jì)算：

式中：ux、uy分別為橫向、縱向的時(shí)均流速。

圖11為工況D1下z?＝0.15（底層）、z?＝0.45（中層）和z?＝0.75（上層）3個(gè)流層水平剖面的渦量分布。整體上可以看出，渦量較大的區(qū)域均產(chǎn)生于堆積體表面和回流區(qū)與主流區(qū)的交界段，這與淹沒丁壩壩后水流結(jié)構(gòu)PIV試驗(yàn)研究［15］得出的結(jié)論一致。渦量大是形成沖坑的主要原因［16］，因此在堆積體表面和回流區(qū)與主流區(qū)的交界段容易造成較為嚴(yán)重的沖刷。堆積體的束窄作用使繞流處存在較大的流速梯度，因此堆積體表面的正渦值較大。堆積體“下大上小”的形態(tài)特征導(dǎo)致水流在近底層（見圖11（a））受到明顯壓縮，流速逐漸增大，主流與下游低流速區(qū)的流速梯度引起剪切流，主流區(qū)水流帶動(dòng)回流區(qū)水流向下游運(yùn)動(dòng)，產(chǎn)生了向下游擴(kuò)散的負(fù)渦。由于底層尚未出現(xiàn)回流，因此下游沿水流方向渦量相對(duì)較小。如圖11（b）、（c）所示，隨著流層的升高，受到主流區(qū)和回流區(qū)水流之間強(qiáng)烈的剪切作用，回流現(xiàn)象逐漸明顯［12］，因此該區(qū)域形成了正渦值較大的渦量帶。同時(shí)，堆積體束窄作用減弱，堆積體表面渦量影響范圍也隨之減小。兩岸邊壁的正負(fù)渦量絕對(duì)值雖不及上述兩區(qū)域大，但也會(huì)對(duì)岸坡造成沖刷影響。文獻(xiàn)［14］指出不可壓縮流體中能量耗散率與流體中渦量絕對(duì)值的大小直接相關(guān)，它與渦量的平方成正比，由此可知水流在經(jīng)過堆積體與繞流的交界處和回流區(qū)與主流區(qū)的交界段后，其能量會(huì)迅速耗散。

圖11 不同流層水平渦量分布

4 結(jié)論與展望

堆積體的存在對(duì)其附近水流結(jié)構(gòu)產(chǎn)生明顯影響，通過水槽試驗(yàn)數(shù)據(jù)驗(yàn)證數(shù)值模擬結(jié)果，在其驗(yàn)證結(jié)果良好的基礎(chǔ)上，展開對(duì)淹沒堆積體附近水流特性的研究，主要結(jié)論如下：

（1）通過分析發(fā)現(xiàn)淹沒堆積體附近水面線沿程變化的一般規(guī)律與非淹沒狀態(tài)下無明顯差異，但堆積體側(cè)向坡度因素在阻流上的影響增強(qiáng)，而進(jìn)占比的影響減弱。通過數(shù)據(jù)擬合發(fā)現(xiàn)，淹沒系數(shù)與相對(duì)壅水水深及相對(duì)最小水深具有良好的指數(shù)關(guān)系。

（2）堆積體上游壅水區(qū)的紊動(dòng)能數(shù)值較小，此區(qū)域容易造成泥沙淤積。水?dāng)U散導(dǎo)致的高紊動(dòng)區(qū)出現(xiàn)在底部，且隨流層上升而逐漸減弱。隨著流層上升，高紊動(dòng)區(qū)的分布和強(qiáng)度開始受到回流現(xiàn)象的主導(dǎo)。

（3）渦量較大的區(qū)域產(chǎn)生于堆積體表面和回流區(qū)與主流區(qū)的交界段。隨著流層升高，堆積體表面渦量影響范圍隨之減小，而回流區(qū)與主流區(qū)交界段的渦量影響范圍隨之增大，在此交界段容易造成較為嚴(yán)重的沖刷。水流經(jīng)這兩個(gè)區(qū)域后，其能量會(huì)迅速耗散。

在模擬研究中，堆積體被設(shè)定為不能透水的剛體，但堆積體在天然狀態(tài)下是由散體砂石組成的，具有透水性，其形態(tài)會(huì)因泥沙起動(dòng)而出現(xiàn)不同程度的改變。在實(shí)際工程問題的模擬研究中，應(yīng)耦合堆積體因泥沙輸移出現(xiàn)的形態(tài)變化過程，對(duì)流場做出相應(yīng)調(diào)整。