官壩河泥石流含沙量變化過程數(shù)值模擬

李俊， 陳寧生， 趙苑迪

(1.四川理工學(xué)院土木工程學(xué)院， 四川自貢643000；2.中國科學(xué)院大學(xué)水利部成都山地災(zāi)害與環(huán)境研究所， 成都610041)

官壩河泥石流含沙量變化過程數(shù)值模擬

李俊1， 陳寧生2， 趙苑迪1

(1.四川理工學(xué)院土木工程學(xué)院， 四川自貢643000；2.中國科學(xué)院大學(xué)水利部成都山地災(zāi)害與環(huán)境研究所， 成都610041)

泥石流輸沙規(guī)律是科學(xué)布置泥石流防治工程的基礎(chǔ)。基于Euler模型模擬了1998年7月6日官壩河泥石流含沙量變化過程。研究結(jié)果有兩點(diǎn)：(1)泥石流輸沙運(yùn)動(dòng)中含沙量大小和溝道寬度成正比關(guān)系。(2)泥石流在溝道地形較陡的地區(qū)攜帶泥沙能力強(qiáng)，在地形較緩的地區(qū)攜帶泥沙能力較弱。研究成果有助于科學(xué)布置官壩河泥石流防治工程，也為解決邛海泥沙淤積問題提供更多的科學(xué)依據(jù)。

西昌市官壩河；泥石流；輸沙過程；Euler方程

引言

對(duì)于大部分泥石流溝，已發(fā)生的泥石流相關(guān)參數(shù)是確定該溝泥石流治理工程設(shè)計(jì)值的重要依據(jù)[1-3]。比如已發(fā)生的泥石流容重、流速、流量和輸沙量是泥石流攔沙壩和排導(dǎo)槽的重要設(shè)計(jì)依據(jù)。所以在泥石流工程防治中，泥石流工作者的首要任務(wù)是調(diào)查某泥石流溝的歷史泥石流發(fā)生時(shí)間，規(guī)模大小，泥石流性質(zhì)和輸沙規(guī)律等[4-9]。本文以邛海官壩河近100年來最大規(guī)模的泥石流為例(即1998年7月6日官壩河泥石流)，模擬該次泥石流輸沙過程，總結(jié)該次泥石流輸沙規(guī)律，為官壩河泥石流工程布置提供科學(xué)依據(jù)，以期通過官壩河泥石流治理工程減少邛海泥沙淤積量。

四川省西昌市的邛海為第四紀(jì)構(gòu)造斷陷淡水湖盆，呈北北西向不規(guī)則長方形展布，湖盆面積27 km2。其中流域面積為137.24 km2的官壩河為匯入邛海水系中匯水面積最大的支流。據(jù)實(shí)地調(diào)查，1998年7月6日官壩河發(fā)生特大型泥石流，泥石流輸移的泥沙一次性把邛海海岸線推進(jìn)了172 m，入海淤積量為68.98×104m3，推進(jìn)淤積面積為0.089 km2，為近百年邛海流域內(nèi)最大規(guī)模的泥石流[10-11]。該次泥石流導(dǎo)致人們賴于發(fā)展的邛海明珠“光澤”逐漸暗淡[12-13]。為減少官壩河泥石流攜帶的泥沙進(jìn)入邛海，本文通過Euler數(shù)值模擬方法研究1998年7月6日官壩河泥石流的輸沙過程，為官壩河泥石流工程布置提供科學(xué)依據(jù)，也為徹底解決邛海的泥沙淤積問題奠定科學(xué)基礎(chǔ)。

1 模擬原理與方法

根據(jù)錢寧團(tuán)隊(duì)的研究成果[14-16]，泥石流可視為固液兩相流。因而Euler方法適合于模擬1998年7月6日官壩河泥石流輸沙過程。通過實(shí)地調(diào)查，獲取主溝(巴拉拉窩段)上游300 m處洪痕斷面，并得出該處洪痕斷面處泥石流洪峰流量。采用Euler模型模擬泥石流含沙量的變化過程，得出該段的泥沙體積分?jǐn)?shù)等場量，根據(jù)實(shí)際溝道泥石流堆積情況驗(yàn)證模擬結(jié)果的正確性，最后分析該段溝道泥石流輸沙過程。在確定泥石流水面線方面，Euler方法的計(jì)算精度低。為了解決這一問題，本文嘗試用計(jì)算自由水面線精度較高的VOF方法確定泥石流漿體的自由水面線，再根據(jù)Euler方法模擬泥石流含沙量的變化過程。利用RNGk-w紊流模型模擬泥石流的運(yùn)動(dòng)規(guī)律，并將泥石流視為具有一定粘度的漿體，即二相流中的液相。這樣即可用VOF得到泥石流漿體與空氣間的自由水面線。該水面線計(jì)算精度高且與實(shí)際調(diào)查結(jié)果相符。在此基礎(chǔ)上，把自由水面線加在已建模型里，形成封閉的有壓流，再用Euler方法確定泥石流含沙量的變化過程。技術(shù)路線如圖1所示。

圖1研究技術(shù)路線

在模擬軟件中建立官壩河局部河流段數(shù)值模型，采用VOF跟蹤計(jì)算官壩河局部河流段泥石流漿體的自由水面，將自由水面線加入已建立的官壩河局部河流段數(shù)值模型，從而建立固液兩相流的Euler數(shù)值模型，最后根據(jù)拉格朗日計(jì)算方法，在多重參考系坐標(biāo)下，選用標(biāo)準(zhǔn)k-ε模型，對(duì)官壩河局部段泥石流輸沙過程進(jìn)行數(shù)值模擬。

2 模擬方程

2.1VOF氣液兩相流模型

連續(xù)方程：

(1)

動(dòng)量方程：

(2)

k方程：

Gk-ρε

(3)

ε方程：

(4)

式中：ρ為液體體積平均的密度，μ為分子粘性系數(shù)，μt為紊流粘性系數(shù)，P為修正的壓力，Gk為由平均速度梯度引起的紊動(dòng)能產(chǎn)生項(xiàng)。

2.2Euler多相流模型

歐拉多相流模型是將顆粒相作為擬流體，認(rèn)為顆粒與流體是共同存在且相互滲透的連續(xù)介質(zhì)，并一同在Euler坐標(biāo)系下處理，為連續(xù)流體模型。其模型包括3個(gè)重要部分：場方程、構(gòu)造方程和界面條件的推導(dǎo)。場方程遵循動(dòng)量、質(zhì)量、能量守恒原理。描述固液兩相流的連續(xù)方程為：

(5)

式中：α1為液體體積分?jǐn)?shù)，v1為液相流速，ρ1為液相的物理密度，ms1為固相與液相間的質(zhì)量傳遞分?jǐn)?shù)。

歐拉多相流模型規(guī)定，體積分?jǐn)?shù)代表了每相所占據(jù)的空間，并且每相獨(dú)自的滿足質(zhì)量和動(dòng)量守恒定律，各項(xiàng)體積分?jǐn)?shù)之和為1。

3 模型建立及條件給定

3.1 主溝(巴拉拉窩段)上游300 m溝道模型建立

數(shù)值模擬的基礎(chǔ)資料為巴拉拉窩上游300 m溝道的地形圖(圖2)，該位置的經(jīng)緯度為N27°48′59.37" E102°26′9.93″，建立該段的模型網(wǎng)格圖，并根據(jù)VOF方法模擬得出該段的自由水面圖(圖3)。

圖2巴拉拉窩上游300 m溝道地形圖

圖3巴拉拉窩上游300 m附近模型

該模型總長度為80 m，泥石流自由水面變化幅度為1.5 m～4 m。根據(jù)實(shí)地調(diào)查，該段泥石流容重為1.77 g/m3。泥石流漿體含有大量的水沙混合物，平均粒徑為10 mm，細(xì)顆粒含量為5%，整體泥石流流動(dòng)狀態(tài)為連續(xù)流。

為保證數(shù)值模擬的來流條件與實(shí)際情況相吻合，模型的網(wǎng)格間距為1 m。在兼顧計(jì)算機(jī)的計(jì)算能力下，計(jì)算區(qū)域網(wǎng)格以結(jié)構(gòu)化網(wǎng)格為主，總的計(jì)算單元約為10.5萬。模型進(jìn)水口泥石流最高水位高程為2340 m，泥石流流量為74.18 m3/s。

3.2 主溝(巴拉拉窩段)上游300 m溝道初始條件和邊界條件

VOF模型中泥石流漿體流量為2700 m3/s，氣相采用天然空氣，液相和氣相體積分?jǐn)?shù)之和為1。計(jì)算泥石流粘滯系數(shù)η：

(6)

其中，rc為泥石流容重，為1.77 g/cm3。計(jì)算得泥石流漿體的粘滯系數(shù)為0.12 kg/(m·s)。

VOF模型中，氣相作為主相，液相作為次相。在巴拉拉窩上游300 m溝道模型進(jìn)口處采用流速控制模擬精度，液相流速為6 m/s，粘滯系數(shù)為0.12 kg/(m·s)。出口處采用大氣壓力控制模擬精度。溝床和溝坡采用無滑移邊界條件。模型上方與空氣接觸的表面設(shè)為開敞邊界。對(duì)于無滑移邊界，各點(diǎn)輸沙率均取為0。需要指出的是該模型未考慮泥石流側(cè)蝕和下蝕作用。

Euler模型中，固相為大于2 mm的泥石流漿體。液相為小于2 mm的泥石流漿體，即為泥沙，其流量為2700 m3/s，粘滯系數(shù)為0.12 kg/(m·s)，固相作為主相，液相作為次相，液相和固相體積分?jǐn)?shù)之和為1。在模型進(jìn)口處采用流速控制模擬精度，液相流速為5.5 m/s。出口處采用大氣壓力控制模擬精度。溝床和溝坡采用固定邊界。溝頂采用大氣壓力控制。對(duì)于固定邊界，各點(diǎn)輸沙率均取為0。

4 模擬結(jié)果驗(yàn)證

根據(jù)Euler方法計(jì)算得出X=5 m和X=75 m處泥沙體積分?jǐn)?shù)(圖4和圖5)，圖4中，X=5.5 m，Z=2234.5 m處的泥沙體積分?jǐn)?shù)為0.6。而在表1中，泥石流在相應(yīng)位置的顆粒百分含量為0.54，與計(jì)算相比小了0.06，這是因?yàn)檎{(diào)查時(shí)間(2010年)距離1998年相差了12年，在這些年中，雨水不斷下滲，帶走了部分泥沙顆粒。所以實(shí)地調(diào)查的泥沙比計(jì)算得出的泥沙體積分?jǐn)?shù)少。在相應(yīng)的深度范圍內(nèi)，實(shí)地調(diào)查的泥沙比Euler方法計(jì)算得出的泥沙體積分?jǐn)?shù)相差較小，說明Euler計(jì)算結(jié)果是與實(shí)際情況較吻合。圖5中，在Z=2.5 m深度范圍內(nèi)，實(shí)地調(diào)查的泥沙比Euler方法計(jì)算得出的泥沙體積分?jǐn)?shù)相差較小，這說明了Euler公式計(jì)算得出泥沙體積分?jǐn)?shù)較為準(zhǔn)確，并且表明這種方法計(jì)算泥石流運(yùn)動(dòng)過程是可行的。通過以上分析，可以采用Euler方法計(jì)算分析官壩河局部段泥石流輸沙過程的變化規(guī)律，為研究官壩河泥石流輸沙機(jī)制奠定基礎(chǔ)。

圖4X=5.5 m時(shí)泥石流漿體與泥沙體積分?jǐn)?shù)分布圖

圖5X=75 m時(shí)泥石流漿體與泥沙體積分?jǐn)?shù)分布圖

取樣位置顆粒百分含量沿深度分布(單位：m)005115225X=5m001004017023054/X=75m001001009015020055

5 主溝(巴拉拉窩段)上游300 m處泥石流輸沙過程分析

泥石流漿體和泥沙體積分?jǐn)?shù)的橫向變化如圖6、圖7所示。

圖6X=15 m時(shí)泥石流漿體與泥沙體積分?jǐn)?shù)分布圖

圖7X=55m時(shí)泥石流漿體與泥沙體積分?jǐn)?shù)分布圖

由圖6、圖7可知，泥沙體積分?jǐn)?shù)在溝底均為0.6%。說明泥石流中泥沙運(yùn)動(dòng)主要集中于泥石流流體中下部，有利于泥沙沉積。另外泥沙在地形較緩的地方體積分?jǐn)?shù)主要聚集在下部，如圖6所示，聚集深度為0.5 m～1.2 m之間。泥沙在地形較陡的地方體積分?jǐn)?shù)主要聚集在中下部，如圖7所示，聚集深度為1.5 m～2.3 m之間，這說明泥石流在溝道地形較陡的地區(qū)攜帶泥沙能力強(qiáng)，在地形較緩的地區(qū)攜帶泥沙能力較弱，此現(xiàn)象符合泥石流的侵蝕規(guī)律。

通過tecplot軟件得出各個(gè)斷面的輸沙量百分?jǐn)?shù)，繪制該段泥石流輸沙量變化圖(圖8)。

圖8官壩河泥石流輸沙量變化過程

在不考慮山坡兩側(cè)泥石流側(cè)蝕和溝底泥石流侵蝕作用所搬運(yùn)的泥沙情況下，泥石流含沙量從斷面1-6的大小先增加后減小，含沙量在斷面3處達(dá)到最大值，為76%。而斷面1-6的溝道寬度分別為18 m，25 m，28 m，31 m，20 m，19 m。溝道寬度也是先增加后減小的曲線。因而含沙量的大小與溝道寬度成正比關(guān)系。

6 結(jié)束語

研究結(jié)果表明，Euler數(shù)值模擬方法能夠較好地模擬官壩河局部溝道泥石流的泥沙分布范圍和體積分?jǐn)?shù)等水力參數(shù)。通過Euler法獲得的泥沙體積分?jǐn)?shù)和泥沙分布范圍，結(jié)合實(shí)際溝道的泥沙分布范圍和泥石流動(dòng)力學(xué)參數(shù)，發(fā)現(xiàn)該溝的泥沙運(yùn)動(dòng)規(guī)律有兩點(diǎn)。(1)官壩河局部溝道泥石流輸沙運(yùn)動(dòng)中含沙量大小和溝道寬度成正比關(guān)系。(2)泥石流在溝道地形較陡的地區(qū)攜帶泥沙能力強(qiáng)，在地形較緩的地區(qū)攜帶泥沙能力較弱。值得說明的是，因該溝道模型較小，可在小流域內(nèi)直接建立Euler計(jì)算模型且精度較高。在中大流域里，由于軟件限制和計(jì)算機(jī)性能限制，不能建立整個(gè)溝道的Euler計(jì)算模型，這時(shí)可采用分段建立模型的方法，分別計(jì)算泥石流輸沙過程中泥沙體積分?jǐn)?shù)和泥沙分布范圍，得出單個(gè)模型的輸沙量，最后把各個(gè)模型的輸沙量累加得出整個(gè)流域的輸沙量。

[1] 齊憲陽,陳寧生,李俊,等.新疆阿合奇縣牙爾巴西溝泥石流特征與防治[J].人民長江,2016,47(2):23-26.

[2] 崔鵬,莊建琦,陳興長,等.汶川地震區(qū)震后泥石流活動(dòng)特征與防治對(duì)策[J].四川大學(xué)學(xué)報(bào):工程科學(xué)版,2010,42(5):10-19.

[3] 胡凱衡,葛永剛,崔鵬,等.對(duì)甘肅舟曲特大泥石流災(zāi)害的初步認(rèn)識(shí)[J].山地學(xué)報(bào),2010,28(5):628-634.

[4] 陳寧生,崔鵬,劉中港,等.基于黏土顆粒含量的泥石流容重計(jì)算[J].中國科學(xué),2003,33(S1):164-174.

[5] 胡桂勝,陳寧生,趙春瑤,等.長江上游梯級(jí)電站開發(fā)區(qū)泥石流的治理工程效果評(píng)估與減災(zāi)策略—以金沙江白鶴灘水電站為例[J].水土保持通報(bào),2017,37(1):241-247.

[6] 李俊,陳寧生,歐陽朝軍,等.扎木弄溝滑坡型泥石流物源及堵河潰壩可能性分析[J].災(zāi)害學(xué),2017,32(1):80-84.

[7] 王怡璇,朱利東,楊文光,等.汶川地震與蘆山地震次生地質(zhì)災(zāi)害特征對(duì)比[J].四川理工學(xué)院學(xué)報(bào):自然科學(xué)版,2015,28(2):63-68.

[8] 舒志樂,史寶寧,張德宇.色多溝泥石流動(dòng)力特征及危險(xiǎn)性評(píng)估研究[J].四川理工學(xué)院學(xué)報(bào):自然科學(xué)版,2016,29(6):70-74.

[9] 吳佳曄.日本東北大地震對(duì)我們的警示[J].四川理工學(xué)院學(xué)報(bào):自然科學(xué)版,2011,24(2):125-128.

[10] CHEN N,CHEN M,LI J,et al.Effects of human activity on erosion,sedimentation and debris flow activity-A case study of the Qionghai Lake watershed,southeastern Tibetan Plateau,China[J].Holocene,2015,25(6):973-988.

[11] WEI X L,CHEN N S,CHENG Q G,et al.Long-term Activity of Earthquake-induced Landslides:A Case Study from Qionghai Lake basin,Southwest of China[J].Journal of Mountain Science,2014,11(3):607-624.

[12] 李俊,陳寧生.官壩河泥石流動(dòng)力學(xué)特征參數(shù)[J].水土保持研究,2013,20(1):264-268.

[13] 余斌,王士革,章書成,等.鵝掌河泥石流對(duì)四川邛海影響的初步研究[J].湖泊科學(xué),2006,18(1):57-62.

[14] 王兆印.河流動(dòng)力學(xué)與河流綜合管理[M].北京:清華大學(xué)出版社,2014.

[15] 錢寧,王兆印.泥石流運(yùn)動(dòng)機(jī)理的初步探討[J].地理學(xué)報(bào),1984,39(1):33-43.

[16] WANG Z Y,QI L J,WANG X Z.A prototype experiment of debris flow control with energy dissipation structures[J].Natural Hazards,2012,60(3):971-989.

Simulation of the Change Process of Sediment Concentration of Debris Flow in the Guanbariver

LIJun1,CHENNingsheng2,ZHAOYuandi1

(1.School of Civil Engineering,Sichuan University of Science & Engineering, Zigong 643000,China; 2.Institute of Mountain Hazards and Environment,CAS,Chengdu 610041,China)

The change process of sediment concentration is the basis for arranging debris flow prevention and control project scientifically.The sediment concentration cannot be obtained from field investigation.Therefore,using simulation method of Euler model,the process of sediment concentration of debris flow in Guanbariver on July 6th1998 is simulated.The resultsare shown as follows.Firstly,sediment concentration of debris flow is proportional to valley width in Guangbariver.Secondly,the debris flow has strong carrying capacity of sediment in areas with steep topography,and has less carrying capacity of sediment in areas with slower topography.The research results are helpful to the scientific arrangement of debris flow prevention project in Guanbariver,and also solve the problem of sediment deposition in Qionghai Lake.

Guanbariver; debris flow; sediment concentration; Euler model

1673-1549(2017)06-0071-05

10.11863/j.suse.2017.06.13

2017-09-12

國家自然科學(xué)基金(41661134012)；四川省安全生產(chǎn)科技項(xiàng)目(aj20170601105926)

李 俊(1989-)，男，四川樂山人，講師，博士，主要從事山地災(zāi)害形成機(jī)理及其防治技術(shù)方面的研究，(E-mail)lijunxiaoyouxiang@163.com

X43

A