潘華利，李炳志,3，鄧其娟，歐國強(qiáng)，孔 玲,3

(1.中國科學(xué)院 水利部 成都山地災(zāi)害與環(huán)境研究所，四川 成都 610041；2.中國科學(xué)院 山地災(zāi)害與地表過程重點實驗室，四川 成都 610041；3.中國科學(xué)院大學(xué)，北京 100049；4.四川建筑職業(yè)技術(shù)學(xué)院，四川 成都 610399)

溝道流體運(yùn)動的動力過程決定著流體容重的變化過程，在流域極端降雨條件下，流體沿溝道流動過程中容重的變化規(guī)律又直接決定著溝道各處流體的類型[1]，從而決定著各溝段是否有泥石流發(fā)生。如果能夠根據(jù)流域所具備的物源條件判斷其在流域極端降雨條件下各溝段是否會暴發(fā)泥石流，對于進(jìn)行流域風(fēng)險評估、防治工程設(shè)計以及潛在性泥石流溝判識等具有重要意義。目前，圍繞泥石流動力過程研究領(lǐng)域，國內(nèi)外學(xué)者開展了大量研究。例如，Cao[2]及Wu[3]等基于淺水流動守恒定律，建立了動床條件下1維潰壩洪水的水沙耦合模型，詳細(xì)研究了水流演進(jìn)過程中水沙與床面的相互作用；Simpson等[4]在Cao等[2]的基礎(chǔ)上引入了泥沙濃度守恒式，獨立處理泥沙的侵蝕及沉積，將水沙耦合模型擴(kuò)展到2維；何思明等[5]基于極限分析的上限定理，研究了黏性泥石流運(yùn)動對溝道土體的侵蝕啟動機(jī)制，得到了泥石流侵蝕啟動速度的表達(dá)式；廖麗萍等[6]從土力學(xué)角度進(jìn)行土體臨界狀態(tài)分析，探討了土體性質(zhì)對泥石流起動的影響；為探索降雨型泥石流對溝床的侵蝕機(jī)理，吳永[7–8]等借助水力學(xué)理論，在科學(xué)表征泥石流流深的前提下，分析了溝床在泥石流流動剪切和滲流水壓耦合作用下的侵蝕機(jī)理；為了對坡面–溝道間土壤侵蝕的時空關(guān)系進(jìn)行合理耦合，高佩玲等[9]對小流域侵蝕動態(tài)過程進(jìn)行分析，建立了流域系統(tǒng)土壤侵蝕動態(tài)過程模擬模型；王鈞[10–11]等采用分布式水文模型得到流域徑流過程，根據(jù)徑流條件計算得到溝道各處固體物源厚度，反演計算了各子流域出口處流體自由表面的變化過程；喬成[12]將流體固液兩相分開考慮，構(gòu)建了考慮顆粒剪脹效應(yīng)和動床侵蝕作用的深度積分泥石流兩相動力學(xué)模型。

上述研究表明，在建立流體運(yùn)動模型時，大多學(xué)者考慮了溝床侵蝕對流體運(yùn)動狀態(tài)的影響，而將溝床侵蝕及坡面匯流綜合考慮，并將其與來流流體合理耦合的研究成果相對較少。本研究從泥石流暴發(fā)所必須的物源條件及水動力條件出發(fā)，在充分分析溝床可移動固體物源侵蝕過程的基礎(chǔ)上，將具有時空變異性的坡面匯流及上空降雨合理地分配于溝道中，將來流流體、被侵蝕的固體物源、坡面匯流以及上空降雨4部分進(jìn)行合理耦合，建立了溝道流體運(yùn)動的水土耦合模型；采用迎風(fēng)差分格式及歐拉公式對模型進(jìn)行數(shù)值離散，得到流體運(yùn)動數(shù)值離散模型；基于MATLAB語言對流體運(yùn)動數(shù)值離散模型進(jìn)行編程求解，可得到各溝段流體容重的動態(tài)變化過程，從而為流域開展防災(zāi)減災(zāi)工作及進(jìn)行潛在性泥石流溝判識等奠定理論基礎(chǔ)。

1 溝道流體運(yùn)動水土耦合模型的建立

1.1 流體運(yùn)動的水土耦合方程

流體沿溝道流動的過程中，除受到沿程坡面匯流及上空降雨兩部分水體補(bǔ)充外，還會侵蝕溝道內(nèi)的可移動固體物源得到物源補(bǔ)充。物源被侵蝕后，會與上游來流、兩側(cè)匯流及上空降雨混合，共同流向下游，致使流體運(yùn)動參數(shù)不斷發(fā)生變化，流體流動示意圖如圖1所示。

圖1 流體沿溝道流動示意圖Fig. 1 Schematic diagram of fluid flow along the channel

假設(shè)匯水坡面為基巖，則坡面匯流只對溝道補(bǔ)充水體，不補(bǔ)充泥沙。同時，假設(shè)可移動固體物源沿溝道方向?qū)挾炔蛔兦椅锢砹W(xué)性質(zhì)均一，則流體中的固相顆粒和液相漿體沿x方向的體積守恒方程可表示為：

式中：方程右邊項是指物源、匯流及上空降雨進(jìn)入流體所引起的固液兩相組分體積的增加；x為笛卡爾坐標(biāo)，m；t為 時間，s；h(x,t)為 流體深度，m；u(x,t)為流體流速，m/s；φ(x,t) 為流體中粗顆粒的體積濃度；i為單位寬度固體物源被侵蝕的速率，m/s；φ1為固體物源中粗顆粒的體積濃度；q(x,t)為坡面匯流及上空降雨對單寬溝道的水體補(bǔ)充量，m/s，其包括坡面?zhèn)攘骷吧峡战涤陜刹糠?，可表示為?/p>

式中：p為降雨強(qiáng)度，m/s；w為溝道寬度，m；qs為坡面?zhèn)攘鲗挝婚L度溝道的水體補(bǔ)充量，m3/s。

固相顆粒和液相漿體動量守恒方程可表示為：

方程（4）和（5）右邊第1項是重力沿斜坡向下的分力，第2項是由于流體自由表面變化所引起的控制體前后兩側(cè)的壓力差，第3項是底床摩擦阻力，第4項是物源進(jìn)入流體所引起的動量變化。

1.2 流域匯水坡面劃分及坡面匯流量計算

流域匯水坡面的水流沿匯水方向從不同點進(jìn)入溝道，使得匯水坡面對溝道的水體補(bǔ)充隨時間和空間而發(fā)生改變。如圖2所示，為了處理具有時空變異性的匯流，本研究根據(jù)匯流區(qū)的坡度和下墊面情況沿溝道方向?qū)R水坡面劃分為多個匯水分區(qū)。采用由美國農(nóng)業(yè)部發(fā)明的用于計算小流域坡面匯流過程的SCS匯流模型分別計算每個匯水區(qū)的江流情況[14–17]。再依據(jù)各匯水分區(qū)與溝道的連接關(guān)系將坡面匯流量分配于溝道中，實現(xiàn)坡面?溝道間的合理耦合。

圖2 小流域匯流區(qū)劃分示意圖Fig. 2 Schematic diagram of the division of small watershed confluence area

采用SCS匯流模型對匯水區(qū)進(jìn)行匯流計算時，首先將一整場降雨劃分為多個凈雨時段，根據(jù)SCS匯流模型，在一次凈雨時段中，峰現(xiàn)時間和洪峰流量分別為：

式（6）～（7）中：tp為 峰現(xiàn)時間，h；qp為洪峰流量，m3/s；Q為一次降雨的凈雨量，mm；A為 匯水區(qū)面積，km2；l為匯流長度，m；y為坡度的百分?jǐn)?shù)；S為流域場次降雨最大滯留量，m；凈雨量Q及場次降雨最大滯留量S可表示為：式（8）～（9）中：I為時段降雨量，mm；CN為反映降雨前流域特征的一個綜合參數(shù)，它與流域前期土壤濕潤程度、坡度、植被等因素有關(guān)。

1.3 溝道可移動固體物源侵蝕速率計算

可移動固體物源被流體侵蝕帶走的過程可看作物源在水動力作用下起動運(yùn)移的過程[18]。物源的穩(wěn)定狀態(tài)與外界降雨情況和地質(zhì)構(gòu)造作用緊密相關(guān)，在地下滲流和地表徑流共同作用下，當(dāng)其所受重力、剪切力等動力大于阻力時，固體物源間會存在一潛在滑動面，形成不穩(wěn)定層[19–21]。此時，物源的重力勢能有向動能和摩擦能轉(zhuǎn)化的趨勢[22]。根據(jù)楊順[23–24]等建立的考慮飽和滲流及表面徑流共同作用下可移動固體物源的臨界力學(xué)判別模型，可將不穩(wěn)定層厚度z表示為：

式中：γm為 流體容重，N/m3；β 為物源內(nèi)摩擦角，(°)；θ為溝道坡角；d50為 物源中值粒徑，m；γs為物源顆粒容重，N/m3；γw為 清水容重，N/m3；n1為 物源孔隙度；C為黏聚力，kPa；γsat為物源飽和容重，N/m3；

由于流體尚未到達(dá)的部分存在不動層的影響，認(rèn)為不穩(wěn)定層z不能瞬時地全體移動，侵蝕達(dá)到z深度需要一定的滯后時間[25]。滯后時間與流體流速、流體容重及堆積體物質(zhì)組成緊密相關(guān)，將滯后時間表示為(αd50) /(uρT) ，則侵蝕速率i為：

式中，α為系數(shù)。

在流體侵蝕作用下，物源厚度不斷減小，溝床地形變化可表示為：

式中，b為溝床地形，m。假設(shè)溝道某點物源初始厚度為d，則經(jīng)過n?t時刻物源剩余厚度dL為：

式中：?t為侵蝕時間步長，s；ii為計算點對應(yīng)時刻物源被侵蝕的速率，m/s。

2 溝道流體運(yùn)動模型的數(shù)值離散

為了對所建的流體運(yùn)動模型進(jìn)行連續(xù)求解，得到溝道各點處各時刻流體運(yùn)動參數(shù)，需對其進(jìn)行數(shù)值離散。數(shù)值離散方法采用有限差分法。其中，空間離散采用一階迎風(fēng)差分格式，時間離散采用歐拉公式。時空離散后，得到溝道流體運(yùn)動數(shù)值離散模型如下：

式（14）～（17）中，上標(biāo)n、n+1指時刻，下標(biāo)j–1、j指控制點。

3 模型初始及邊界條件處理

3.1 初始條件

初始時，流域無降雨，溝道內(nèi)沒有流體流過，初始條件為：

由于h要滿足計算公式

故在初始時刻以一無窮小值 ε代 替h。

3.2 邊界條件

1）匯水坡面與溝道的交匯處為模型內(nèi)邊界，溝道在交界處會受到坡面匯流的水體補(bǔ)充，匯流補(bǔ)充量為qs(x,t)，可根據(jù)所劃分的各匯水分區(qū)采用SCS匯流模型計算得到。

2）溝道起點處流體的運(yùn)動參數(shù)可根據(jù)溝口匯流區(qū)計算得到，起點邊界條件可表示為：

得到模型初始條件及邊界條件后，將其代入溝道流體運(yùn)動數(shù)值離散模型，采用基于MATLAB語言編寫的流體運(yùn)動模型求解程序?qū)ζ溥M(jìn)行求解，便可得到溝道各點處各時刻流體的流深、流速、顆粒濃度及容重值。

4 模型驗證

通過室內(nèi)水槽實驗對流體運(yùn)動水土耦合模型進(jìn)行驗證，水槽實驗在中國科學(xué)院、水利部成都山地災(zāi)害與環(huán)境研究所的山地災(zāi)害與地表過程重點實驗室中進(jìn)行。

4.1 實驗?zāi)Ｐ?/h3>

如圖3所示，實驗裝置主要由供水裝置、溝道模擬裝置及降雨模擬裝置組成。供水玻璃水槽尺寸為200 cm×20 cm×25 cm，用于給模擬溝道供水；溝道模擬水槽尺寸為200 cm×20 cm×45 cm，用于堆放實驗物料，模擬泥石流溝道；溝道模擬水槽起點至100 cm處設(shè)置了不透水混凝土棱柱模型，以防止溝槽起點處物源被淘蝕而形成跌坎，避免過流流體紊動劇烈；降雨模擬裝置位于溝道模擬裝置正上方，其降雨噴頭分布長度為120 cm，用于給溝道模擬水槽提供水體補(bǔ)充，模擬野外降雨及坡面匯流對實際溝道的水體補(bǔ)充。

圖3 實驗裝置圖Fig. 3 Experimental equipment

如圖4所示，實驗物料堆放長度為200 cm，堆放時使得物料上表面與供水水槽底面處于同一平面。為了消除物料與實驗水槽底部明顯的分界面的影響，在實驗水槽出口處設(shè)置了高5 cm的侵蝕基準(zhǔn)。

圖4 實驗物料照片F(xiàn)ig. 4 Experimental material photos

實驗選取來流流量及溝槽坡度為控制變量，對所建的流體運(yùn)動水土耦合模型進(jìn)行驗證。來流流量大小通過動力泵控制，動力泵可控流量范圍為0～3 000 cm3/s，流量大小隨動力泵開關(guān)的開合程度而發(fā)生變化，本實驗根據(jù)動力泵開關(guān)可動范圍設(shè)計6個檔位，測量得到各檔位對應(yīng)流量分別為598、967、1 401、1 745、2 113、2 594 cm3/s。溝道模擬水槽坡度可調(diào)范圍為0～15°，為在滿足泥石流暴發(fā)所需的坡度條件下探索坡度對流體運(yùn)動狀態(tài)的影響，在坡度可控范圍內(nèi)均勻地設(shè)計5組變化值，分別為5°、7°、9°、11°、13°。同時，降雨模擬裝置提供的水體補(bǔ)充量為310 cm3/s。侵蝕參數(shù)α和溝床表面糙度n為未知值，需要對其進(jìn)行反復(fù)調(diào)校，找到最合適的參數(shù)值。采用12組水槽實驗數(shù)據(jù)對模型參數(shù)進(jìn)行調(diào)校，發(fā)現(xiàn)當(dāng)侵蝕參數(shù)α為1.5×107，溝床表面糙度n為0.012 5時，流體運(yùn)動模型的計算值與實驗值有著比較高的吻合度。相關(guān)實驗參數(shù)明細(xì)見表1。

表1 實驗參數(shù)明細(xì)表Tab. 1 Experimental parameters

4.2 實驗土體

實驗過程中，所使用的物料取自于四川省都江堰市龍池鎮(zhèn)干家溝，所用物料最大粒徑為dmax=20 mm，實驗水槽寬度與最大粒徑比值為10，符合寬徑比的要求[26]。實驗前，抽取3組樣品采用篩分法進(jìn)行顆粒級配分析試驗，得到物料級配圖，如圖5所示。

圖5 干溝泥石流原樣顆粒級配曲線Fig. 5 Particle grading curve of debris flow sample in Gangou

實驗物料基礎(chǔ)物理力學(xué)性質(zhì)參數(shù)如表2所示。

表2 實驗物料基礎(chǔ)物理力學(xué)性質(zhì)參數(shù)Tab. 2 Physical and mechanical property parameters of experimental materials

其中，實驗物料顆粒密度采用2 000 mL量筒通過水中置換法進(jìn)行3次測量求取平均值得到；飽和密度通過環(huán)刀取3組樣品，然后向樣品中滴水至飽和狀態(tài)，測量后求取密度平均值得到；內(nèi)摩擦角及黏聚力的測量如圖6所示，通過由直剪試驗在4種不同垂直壓力下的應(yīng)力應(yīng)變曲線所繪制的莫爾圓得到。

圖6 直剪試驗測量內(nèi)摩擦角及黏聚力Fig. 6 Measurement of internal friction angel and cohesion by direct shear test

4.3 實驗量測參數(shù)

選取溝道模擬水槽控制點，測量不同工況條件下各控制點處流體流深、流速和流體容重，將其與理論模型計算值進(jìn)行對比，以驗證所建理論模型的合理性。如圖7所示，本實驗共選取2個控制點，為避免降雨模擬裝置對流深及流速的測量造成影響，控制點應(yīng)選在降雨模擬裝置分布以外的區(qū)域，故將1號控制點設(shè)置于溝道模擬水槽下游140 cm處，2號控制點設(shè)置于溝道模擬水槽出口處。

圖7 控制點分布圖Fig. 7 Distribution map of control points

實驗過程中，采用攝影法結(jié)合激光測距法對流體流深進(jìn)行測量。在實驗水槽側(cè)面安裝索尼HDR–PJ670高速攝影機(jī)記錄流體運(yùn)動過程，通過影像解譯獲取控制點處各時刻流體流深。同時，于控制點上方安裝測量精度為1.5 mm的RLM–S30C激光測距儀，測量流體自由表面高程的變化過程，并與通過攝影得到的流體液面高程變化過程進(jìn)行對比，以保證通過影像解譯得到的流體流深的準(zhǔn)確性。流速的測量采用浮標(biāo)法，于供水水槽上游投擲乒乓球浮標(biāo)，通過水槽側(cè)面的高速攝影機(jī)記錄浮標(biāo)流至控制點處的運(yùn)動過程，通過影像解譯計算得到控制點處各時刻流體的流速。流體容重通過在出口處取樣得到，當(dāng)流體流出溝道出口時，不斷用量筒于出口處取樣，同時通過水槽前方的攝影機(jī)記錄各樣本的取樣時刻，待實驗完成后測量樣本的質(zhì)量及體積，計算得到溝道出口處各時刻流體的容重，并取平均值即得到當(dāng)前工況條件下溝槽出口處流體的容重值。流體運(yùn)動參數(shù)測量的實驗過程照片如圖8所示。

圖8 實驗過程照片F(xiàn)ig. 8 Experimental process photos

4.4 實驗與理論結(jié)果對比

1）流體流深實驗結(jié)果與理論結(jié)果對比分析

圖9、10分別為不同來流流量條件下兩個控制點處流體流深實驗值與計算值的變化及對比。圖9和10表明：不同來流流量條件下，溝槽兩個控制點處流深實驗值與計算值變化趨勢相同，均隨來流流量的增加而不斷增大；但流深的實驗值略大于計算值。其原因是：流體從供水水槽流入物源區(qū)時，會先對接覆土較薄的物源區(qū)進(jìn)行沖刷，致使混凝土棱柱模型與物源堆積體的銜接面下移，流體流經(jīng)銜接面時，過流表面的變化會增加流體的紊動作用；同時，物源堆積體表面的凹凸不平也會促使過流流體紊動作用增強(qiáng)，致使通過高速攝影機(jī)觀測到的流體流深略大于其計算流深。對流深進(jìn)行誤差分析，發(fā)現(xiàn)控制點1和出口處流深實驗值與計算值的相對誤差分別為14.3%和16.7%。

圖9 流體流深變化Fig. 9 Evolution of fluid depth

圖10 流深實驗值與計算值對比分析Fig. 10 Comparative analysis chart of the experimental value and calculated value of fluid depth

2）流體流速實驗結(jié)果與理論結(jié)果對比分析

不同流量條件下流速實驗結(jié)果與理論結(jié)果對比情況見圖11、12。由圖11和12可見，兩個控制點處流速的實驗值與計算值吻合良好，且均隨來流流量的增加而不斷增大。其原因是：來流流量越大，單位時間通過溝道橫截面的流體越多，所對應(yīng)的溝道各點處的流速便越大。誤差分析發(fā)現(xiàn)控制點1和出口處流速實驗值與計算值的相對誤差分別為5.7%和6.1%。

圖11 流體流速變化Fig. 11 Evolution of fluid velocity

圖12 流速實驗值與計算值對比分析Fig. 12 Comparative analysis chart of the experimental value and calculated value of fluid velocity

3）流體容重實驗結(jié)果與理論結(jié)果對比分析

不同流量條件下容重實驗結(jié)果與理論結(jié)果對比情況如圖13、14所示。由圖13和14可見，流體容重實驗結(jié)果與理論結(jié)果吻合度較高，且均隨來流流量的增加而呈逐漸減小的趨勢，這是來流流量的增加對流體產(chǎn)生的稀釋作用所致。對其進(jìn)行誤差分析，發(fā)現(xiàn)出口處容重實驗值與計算值的相對誤差為5.3%。

圖13 出口處流體容重變化Fig. 13 Evolution of fluid bulk density at the outlet

圖14 出口處容重實驗值與計算值對比分析Fig. 14 Comparative analysis chart of experimental value and calculated value of bulk density at the outlet

綜上所述，溝道各點處流體流深、流速和容重實驗結(jié)果與理論結(jié)果的相對誤差均小于20%。其中，流體流速和容重的模擬精度均高于90%，流體流深的模擬精度高于80%，表明溝道流體運(yùn)動的水土耦合模型侵蝕參數(shù)的確定較為合理，采用該模型對溝道流體運(yùn)動過程進(jìn)行模擬基本可行。

5 結(jié) 論

本研究包含以下研究成果：

1）沿著溝道方向，根據(jù)匯水坡面坡度及下墊面情況劃分匯水區(qū)，根據(jù)坡面與溝道間的連接關(guān)系將坡面匯流量分配于溝道之中，實現(xiàn)了坡面?溝道間時空關(guān)系的合理耦合，在此基礎(chǔ)上建立了考慮溝床侵蝕和坡面匯流的溝道流體運(yùn)動的水土耦合模型。采用歐拉公式和迎風(fēng)差分格式對流體運(yùn)動水土耦合模型進(jìn)行時空離散，得到了溝道流體運(yùn)動數(shù)值離散模型。

2）通過分析水槽實驗數(shù)據(jù)，發(fā)現(xiàn)溝道各點處流體流深、流速與來流流量呈正相關(guān)關(guān)系，流體容重與來流流量呈負(fù)相關(guān)關(guān)系。

3）對比分析模型實驗結(jié)果與理論計算結(jié)果，發(fā)現(xiàn)流體流深、流速和容重的實驗值與模擬值的相對誤差均小于20%。流體流速和流體容重的模擬精度高于90%，流體流深的模擬精度高于80%，表明該小流域溝道流體運(yùn)動水土耦合模型及其數(shù)值計算方法基本合理，采用該模型模擬動床作用下溝道流體運(yùn)動過程基本可行。