假冬冬,楊 俊,陳長英,張幸農,應 強

(南京水利科學研究院水文水資源與水利工程科學國家重點實驗室,江蘇 南京 210029)

岸灘側蝕崩塌現象廣泛分布于世界各大江河之中[1-3],是一種危害較大的自然現象,與水沙動力等因素共同引起復雜河床沖淤調整[4]。對于中國東北地區(qū)的季節(jié)性冰凍河流(季凍區(qū)河流)而言,河道岸灘、島嶼汊道不穩(wěn)定,凍融作用強、崩岸頻發(fā)[5-6],直接引起農田土地流失。此外,季凍區(qū)河流受水動力-凍脹/凍融耦合作用,其岸灘崩塌與一般河流存在明顯差異,機理更為復雜。因此,開展季凍區(qū)河流岸灘崩塌機理及其數值模擬研究,可為江河治理與保護提供科技支撐,具有重要的理論和實踐意義。

岸灘崩塌是指在近岸水沙與河床邊界的相互作用下,河岸受到各種因素影響而發(fā)生部分土體崩塌的現象[7]。隨著對崩岸現象的深入認識,國內外學者在崩岸機理與力學模式方面開展了具有重要價值的研究工作,建立了不同類型土質河岸的岸灘崩塌描述模式[8-11],為河道崩岸的數值模擬提供了力學基礎?；诖?不少研究者將岸灘崩塌力學模式與傳統(tǒng)水沙輸移模型相耦合,開展了岸灘侵蝕過程的二維數值模擬和三維數值模擬研究,有效復演了岸灘崩塌及河道演變過程[12-15],并在長江、黃河等河流岸灘失穩(wěn)過程中得到成功應用。此外,傳統(tǒng)水沙數學模型在松花江等季凍區(qū)河流河道治理中也得到廣泛應用,例如:宗原等[16]對黑龍江雪水溫河段護岸的防護效果進行了研究,護岸的修建可有效控制河勢穩(wěn)定，防止因河岸崩塌引起的國土流失；陸永軍[17]采用二維泥沙數學模型,模擬了松花江干流三姓淺灘航道整治一期工程丁壩群作用下的回流邊線、汊道分流比以及推移質運動引起的河床變形,但此類模型并未考慮季凍區(qū)河流岸灘崩塌機理。從當前岸灘崩塌機理及其數值模擬的研究進展來看,主要以水動力沖刷下岸灘崩塌機理及其與河床沖淤交互作用為主,尚未深入考慮水動力-凍脹/凍融耦合作用,但對于季凍區(qū)河流而言,受氣候影響則存在顯著凍脹/凍融特點,且對岸灘力學強度和穩(wěn)定性具有重要影響[18-20]。除此之外,對于松花江等季凍區(qū)河流而言,凍融作用強,岸灘崩塌現象常見,但受模擬技術的限制,其河道演變模擬尚未考慮岸灘崩塌過程,亦未涉及凍脹/凍融作用對岸灘崩塌的影響。

本文針對季凍區(qū)河流岸灘崩塌特點,以松花江干流大頂子山航電樞紐下游近壩段為例,在三維水沙數學模型研究成果[14-15]的基礎上,通過建立水動力-凍脹/凍融耦合作用岸灘崩塌力學描述模式,提出季凍區(qū)河流岸灘崩塌及河床變形的數值模擬方法,采用實測數據對數學模型進行驗證,在此基礎上分析岸灘崩塌對河床沖淤的影響。研究成果可為季凍區(qū)河流岸灘崩塌及河道演變的深入研究提供技術支撐。

1 水動力-凍脹/凍融耦合作用岸灘崩塌機理

對于季凍區(qū)河流而言,岸灘侵蝕及其演變過程除了受水沙條件影響外,因氣候影響還存在顯著凍脹/凍融特點,且對岸灘穩(wěn)定性造成重要影響,機理更為復雜,傳統(tǒng)河流動力學理論已不能完全適用。每年封凍初期,河岸土體內的水體在負溫度條件下結冰,土在凍結時易發(fā)生體積膨脹,即發(fā)生凍脹現象。凍脹作用破壞了河岸土層的結構,土壤結構變松,經常使得岸灘土層出現裂縫,甚至坍塌。而到次年春天,當上升到正溫度時,因凍結土層的冰融化,作為土體臨時結構黏結力的冰晶體消失,土體結構變松,力學強度降低,對于處于極限平衡狀態(tài)的岸坡而言,凍融作用會進一步使河岸邊坡處于不穩(wěn)定狀態(tài),最終導致河岸坍塌。

基于上述分析可知,凍脹作用通過凍脹力會使岸坡開裂,而凍融過程則會使岸灘土體力學強度降低,兩者均會影響岸坡的穩(wěn)定。因此,對于季凍區(qū)河流岸灘崩塌而言,既要考慮水動力對岸灘的沖刷侵蝕,又要考慮凍脹/凍融對岸灘穩(wěn)定的影響機理。因此,需建立坡腳水力沖刷計算模式以及凍脹/凍融耦合作用下岸灘失穩(wěn)模式。

1.1 坡腳水力沖刷計算模式

水流淘刷作用下,坡腳側蝕沖刷過程可由筆者提出的岸灘側蝕速率公式計算[11]:

(1)

式中:ωb為岸灘側蝕速率；λ為側蝕系數,由實測資料進行率定；γ為水體重度；γb為岸坡土體重度；u為近岸處水流流速；uc為泥沙起動流速。

1.2 岸灘失穩(wěn)模式

1.2.1 上部土體裂縫深度計算

在水力沖刷與凍脹作用下,坡頂開裂,會形成一定深度的裂縫,并影響岸坡穩(wěn)定。根據土力學朗肯理論,考慮凍脹作用時,岸灘坡頂裂縫深度為

(2)

式中:Ht為考慮凍脹作用的坡頂裂縫深度；pd為凍脹力；c為土體黏聚力；Kz為朗肯主動土壓力系數。

凍脹力與土體性質(黏聚力)密切相關,設ξ=pd/(2c),則式(2)可簡化為

(3)

分析認為,凍脹作用會使岸灘坡頂開裂深度進一步增加,從而降低岸灘的穩(wěn)定性。

1.2.2 岸灘崩塌力學模式

水動力-凍脹/凍融耦合作用下岸灘崩塌過程,同樣基于岸灘受力分析,可推導獲得岸灘崩塌力學模式。從三維受力分析角度,考慮沿縱向相鄰土體對其穩(wěn)定性的影響,筆者借鑒Osman和Thorne[21]的思路,推導建立了岸灘崩塌力學模式[15]。與此類似,考慮凍脹對土體上部開裂深度影響以及凍融對土體力學強度影響,得到水動力-凍脹/凍融耦合作用下岸灘崩塌力學模式。

(1)初次崩塌。在已知季凍區(qū)河流初始岸灘高度(H0)、初始岸坡(i)的情況下,根據水沙數學模型可計算得到坡腳橫向沖刷值為ΔB,河床沖深為ΔZ，并根據垂向、橫向變形計算結果,修正岸灘形態(tài)。圖1為發(fā)生沖刷后初次崩塌的示意圖。初始岸灘地形由ABCD連線表示，經水流沖刷及崩岸后的岸灘地形變?yōu)锳′B′C′D。FS為土體滑動力,FR為土體抗滑力,W為岸灘單位寬度坍塌土體的重量,WB為坍塌土體的寬度,β為發(fā)生初次崩塌時,滑動面與水平面的夾角,可由式(4)計算:

圖1 岸灘初次崩塌Fig.1 Sketch of the initial bank failure

(4)

式中:H為沖刷后的岸灘高度；H′為岸坡線轉折點以上的岸灘高度；k=Ht/H；i為初始岸坡；φ為岸灘土體的內摩擦角。

對于縱向長度為Δl的岸灘土體,假設FS=FR時,岸灘處于崩塌的臨界狀態(tài),由此可求得岸灘崩塌的臨界條件為:

(5)

式中:Δl為岸灘崩塌的縱向長度,pi-1、pi+1分別為上、下游側面所受黏聚力系數。

在計算得到臨界值(H/H′)cr和實際值(H/H′)m的情況下,按如下方法可判斷岸灘是否會崩塌:當(H/H′)m<(H/H′)cr時,岸灘穩(wěn)定,不是發(fā)生岸灘崩塌的臨界條件,進入下一個時段的水沙計算；當(H/H′)m≥(H/H′)cr時,岸灘失穩(wěn),利用幾何形態(tài)關系,可計算坍塌土體尺寸。

(2)二次崩塌。岸灘發(fā)生前述初次崩塌后,并假定之后的岸灘崩塌方式為平行后退,即岸灘崩塌的破壞角仍為β,崩塌示意圖見圖2。

圖2 岸灘二次崩塌Fig.2 Sketch of the subsequent bank failure

此時,假設FS=FR時,可獲得岸灘二次崩塌臨界條件為

(6)

與初次岸灘崩塌判別類似,在計算得到臨界值(H/H′)cr和實際值(H/H′)m的情況下,即可判斷岸灘是否發(fā)生崩塌。

(3)凍融作用引起的岸灘土體力學性能變化。土體凍融過程中,因凍結土層的冰融化,作為土體臨時結構黏結力的冰晶體消失,土體結構變松,強度降低。其中，抗剪強度是對岸灘穩(wěn)定性具有重要影響的力學性質指標,通常采用黏聚力和內摩擦角來表示。根據已有研究表明,凍融循環(huán)作用下土體黏聚力顯著減小,但內摩擦角的變化相對不明顯[18-19]。根據文獻[19]的研究,多根樣品曲線試驗結果表明:黏聚力隨著凍融循環(huán)次數的增加而逐漸減小,單次凍融后黏聚力的減幅范圍約為15%～40%；而內摩擦角無規(guī)律可循呈上下波動的狀態(tài)，單次凍融后內摩擦角有增有減,變幅不大。根據文獻[20]的研究,單根樣品曲線試驗結果表明:單次凍融后黏聚力的減幅約為19%；單次凍融后內摩擦角的增幅約為14%。因此,凍融作用對岸灘土體抗剪強度和穩(wěn)定具有重要作用,本次模擬研究中將考慮凍融影響,對抗剪強度進行相應調整。

松花江干流屬于典型的季節(jié)性冰凍河流,每年4月至11月一般為暢流期,造床作用相對較大,12月至次年3月一般為封凍期,河床沖淤相對較小。針對這一特點,本研究主要對暢流期進行模擬,計算過程中對其進行概化處理:4月考慮凍融影響,即抗剪強度相應調整；11月考慮凍脹作用,即式(3)中考慮凍脹作用,坡頂開裂深度增加,從而降低岸灘的穩(wěn)定性；其余時段模擬則考慮常規(guī)岸灘穩(wěn)定計算。

2 岸灘崩塌與河床沖淤數值模擬技術

岸灘崩塌與河床沖淤的模擬,一方面需要對傳統(tǒng)水沙輸移和河床垂向沖淤過程進行計算,另一方面還需對岸灘崩塌過程進行模擬,并實現兩者的有機聯系。本節(jié)將對季凍區(qū)河流岸灘崩塌與河床沖淤的模擬技術進行介紹。

2.1 局部網格可動技術

與常規(guī)河流模擬類似,在應用數學模型對季凍區(qū)河流岸灘崩塌及河床沖淤過程進行模擬時,崩岸的發(fā)生會改變數值模型的地形和邊界,而岸灘崩塌寬度不一定與崩岸處的模型網格寬度一致,計算網格很難準確跟蹤崩岸后的地形和邊界。本文采用一種基于非正交網格的局部網格可動技術對岸灘崩塌過程進行岸坡邊界的跟蹤,在整個大計算域內生成計算網格。在模擬過程中,僅對岸灘崩塌附近的網格進行移動,實時準確地跟蹤新的岸坡邊界,其余計算網格位置保持不變。這樣做既可較準確地實時跟蹤岸灘崩塌后的岸坡邊界,又無需重新生成整個計算域內的模型網格,可彌補傳統(tǒng)定網格和動網格在這方面存在的不足之處。具體過程可詳見文獻[14]。

2.2 數值模擬流程

本文模擬過程中采用的三維水流泥沙模型為非均勻非平衡泥沙輸移模型,模型基于非正交曲線網格系統(tǒng),采用控制體積法進行方程離散與求解[14-15],該模型已在實際天然河流中得到了較好的驗證和應用。岸灘崩塌模擬過程中,首先給定各變量初值,采用三維水沙模型計算河道內流場、水位等水動力信息；依據水動力模擬結果,計算懸移質含沙量及推移質輸沙率；隨后計算河床沖淤變形,并更新地形；在進行凍融/凍脹作用判斷與計算的基礎上,根據岸灘崩塌機理與力學模式,計算和判斷岸灘穩(wěn)定狀態(tài),若岸灘失穩(wěn),則修改岸灘及近岸地形以及河床級配信息；依次重復上述計算過程,直至完成設置的計算時間。

3 典型季凍區(qū)河流岸灘崩塌數值模擬

以松花江干流大頂子山航電樞紐下游近壩段為例,進行典型季凍區(qū)河流岸灘崩塌數值模擬分析。松花江干流屬于典型的平原沖積河流,全長約940 km,區(qū)間集水面積18.64萬km2,兩岸除部分區(qū)段為丘陵崗地外,基本為沖積平原,向東流至同江附近匯入黑龍江。大頂子山航電樞紐工程位于松花江干流中游、哈爾濱市區(qū)下游約63 km處,是松花江航道梯級開發(fā)總體規(guī)劃中8座樞紐之一,工程于2004年9月開工建設、2008年末交工使用。樞紐的運行以及上游來沙減小,使得樞紐下游近壩段出現較為明顯的沖刷,深槽擺動,洲灘崩退現象時有發(fā)生。本節(jié)將依據前文模擬技術,建立大頂子山航電樞紐下游近壩段河道演變動力學模型,并采用實測資料對模型進行驗證,分析岸灘崩塌對河道演變計算結果的影響。

3.1 計算范圍及計算條件

計算河段上起大頂子山航電樞紐,下至富江島匯流口下游順直段,模擬河段全長約22 km?；?008年10月和2014年10月實測河床地形構建數學模型,模型采用貼體曲線非正交網格與復雜岸線貼合良好,模型平面計算網格為880×280,平均網格尺度沿水流方向約為25 m、沿斷面方向約為10 m,并在局部區(qū)域對計算網格進行適當加密,使網格能夠較好地反映復雜河道邊界,模型垂向網格分為13層。模型平面計算網格劃分及河勢分別見圖3、圖4。模擬時段內模型進口流量過程如圖5所示,出口水位由率定的水位—流量關系確定(圖6)。松花江干流屬于少沙河流,主要由推移質造床為主,且計算河段位于大頂子山航電樞紐下游近壩段,樞紐運行初期模型進口不考慮上游懸沙。

圖3 計算網格劃分Fig.3 Computational meshes

圖4 研究河段河勢Fig.4 River regime of the study area

圖5 模型進口流量過程Fig.5 Flow discharge processes of model inlet

圖6 計算水位—流量關系與實測水位—流量關系對比Fig.6 Calibration of relationship between water level and discharge

3.2 模型計算分析與討論

(1)水位—流量關系驗證。為使模型能夠準確復演天然河道水流運動狀態(tài),首先采用模擬河段進口(大頂子山航電樞紐壩址)及下游(模型出口處)的水位—流量關系進行計算,并與天然水位—流量關系進行比較。各級流量下模型計算水位—流量關系與天然水位—流量關系對比情況見圖6,由圖可見,模型計算水位—流量關系與天然情況基本吻合,各級水情下水位差值最大一般不超過0.10 m,說明模型能夠較好地復演河道阻力情況,率定的河段綜合糙率一般為0.022～0.035,主槽阻力相對較小,灘地阻力相對較大(綜合糙率一般在0.030～0.050左右)。

(2)典型流場分布特征。本河段典型流量下(多年平均流量Q=1 200 m3/s)流場分布特征如圖7所示。由圖可見,中水期河道內水流流速一般在1 m/s以內,順直段流速較為平順,彎曲段三維流場特征較為明顯,表層、中層流速大于底層流速,表層流速指向彎道凹岸、底層流速指向彎道凸岸,局部最大夾角可達40°左右。江中的富江島將該處河道分為左右兩汊,其中左汊為支汊、右汊為主汊,兩汊分流之比約為1∶4；兩汊河道均較為彎曲,凹岸側流速較凸岸側大,彎頂處尤為明顯,易造成河道凹岸沖刷。從三維流場中截取典型斷面(CS1、CS2)流速分布(由下游往上看)如圖8所示,斷面位置見圖7(b)；從圖中可看出,彎道斷面環(huán)流特征較為明顯,表層流速流向凹岸、底層流速流向凸岸。

圖7 流速分布特征Fig.7 Computed velocity distribution

圖8 典型斷面流速分布特征Fig.8 Computed cross-sectional velocity distribution

(3)典型斷面沖淤驗證分析。根據實測河道地形資料,對研究河段的河床沖淤變化進行驗證。收集了計算河段2008年10月實測地形資料,作為模型計算初始地形,根據2014年10月實測的典型大斷面資料對模型進行驗證。根據現場采樣分析,岸灘土體中值粒徑約為0.09 mm,級配曲線如圖9所示,黏聚力以及重度分別取值為14.0 kN/m2和18.0 kN/m3,內摩擦角取值23°,模擬過程中采用此值作為基準值。凍融期對抗剪強度基準值進行相應調整,本文首先根據文獻[19]的研究進行調整,即單次凍融后黏聚力減幅為30%(取均值),內摩擦角保持不變。

圖9 岸灘土體級配Fig.9 Sediment gradation of riverbanks

大頂子山航電樞紐運行初期,依據收集到的資料,下游近壩段典型斷面沖淤變化驗證見圖10,斷面位置見圖4。由圖可見,總體而言2008年10月至2014年10月各斷面以沖刷為主,近壩處斷面沖刷尤為明顯,斷面平均沖深在2 m左右,最大沖刷在4 m以上,下游彎道處斷面彎道凹岸沖刷明顯,岸線崩退,最大崩退距離可達150 m左右,彎道凸岸則出現淤積。具體而言,從各典型斷面的沖淤變化來看,DM1位于樞紐下游約1.5 km處,樞紐運行初期(2008—2014年)河床在整體沖刷下切的同時,右側深槽相鄰河岸崩退距離約100 m；DM2位于樞紐下游約3.0 km 處,與DM1類似,河床整體沖刷下切約2.0 m,右側河岸后退距離約150 m；DM3位于樞紐下游約6.0 km處,河床表現為左側沖刷,右側淤積,總體上仍呈現為沖刷,平均沖刷1.3 m左右；DM4位于樞紐下游約10.0 km處,位于彎道進口段,受水流取直影響,右側岸坡不斷崩退,崩退距離約200 m,河道展寬明顯；DM5位于樞紐下游約12.0 km處,位于彎道段,受彎道環(huán)流和水流頂沖影響,凹岸側河床沖刷、岸線后退約140 m,凸岸側河床則出現淤積；DM6位于富江島右(主)汊彎道處,凸岸出現淤積,凹岸則沖刷明顯；DM7位于富江島兩汊匯流處下游,河道中部沖刷,兩側有所淤積。由上述各斷面沖淤驗證結果可看出,考慮岸灘崩塌的模擬結果與實測結果沖淤定性基本一致,沖淤幅度基本相當,數學模型可較好地反映河道的沖淤變化、岸灘崩退以及河道展寬過程。

圖10 典型斷面沖淤對比分析Fig.10 Simulated and measured results of deposition and erosion distribution at typical cross-sections

(4)岸灘崩塌對斷面沖淤影響分析與討論。為分析岸灘崩塌對斷面沖淤的影響,對考慮與不考慮岸灘崩塌的模擬結果進行了對比,見圖10。從結果對比來看,兩者存在較明顯的差異,考慮岸灘崩塌的模擬結果與實際情況相對吻合,由此也說明岸灘崩塌對于季凍區(qū)天然河道演變準確模擬的重要性。相對而言,順直段兩者差別較小,河床主要以沖刷下切為主,彎曲段受水流取直和頂沖作用,岸灘沖刷后退現象較明顯,考慮岸灘崩塌后的模擬結果與實際觀測較為一致,而未考慮岸灘崩塌的模擬結果則與實際情況存在較大差異。

為進一步分析凍融期抗剪強度調整對岸灘崩退的影響,除了前述根據文獻[19]的研究進行調整外,即單次凍融后黏聚力減幅為30%(取均值),內摩擦角保持不變；還根據文獻[20]的研究,對岸灘崩退進行模擬,即單次凍融后黏聚力減幅為19%,內摩擦角增幅為14%。前者簡稱工況1,后者簡稱工況2。對于岸灘崩退不明顯的斷面而言,兩工況的模擬結果基本一致；為方便比較兩者模擬結果之間的差異,選取具有典型岸灘崩退的DM4、DM5進行對比分析,見圖11。由圖可見:兩工況均有較明顯的岸灘崩退過程,其差異并不十分明顯,主要是工況1的岸灘崩退距離較工況2略大,最大差異在15 m左右。這可能是因為凍融期,河道流量和流速相對較小,河床沖刷幅度亦相對小,岸灘形態(tài)主要由前期河床變形所主導,只有當前期岸灘已處于或接近于臨界崩塌狀態(tài)時,此時抗剪強度降低才易引發(fā)岸灘崩塌。

圖11 凍融抗剪強度變化對典型斷面岸灘崩塌影響Fig.11 Effects of freeze-thaw on bank erosion at typical cross-sections

同時也應指出,盡管考慮岸灘崩塌后的模擬結果與實測結果相對吻合,但具體定量上的沖淤分布還存在一定差異,一方面是因為季凍區(qū)河流岸灘崩塌機理極為復雜,涉及水、土、氣候等方面因素,不僅存在平行崩岸類型的崩退,也存在淘刷落崩的類型,而本文只考慮了平行崩退模式,且實際土質條件與模擬設置的土體抗沖性存在差異,這會對模擬結果產生影響,其精細描述模式及其模擬方法還有待進一步深入研究；另一方面,在春季融化時期,北方河流普遍存在冰凌漂浮現象,雖然持續(xù)時間不長,但冰凌漂浮對岸灘穩(wěn)定性也存在較明顯的影響,本文模型尚未對此進行考慮；此外,季凍區(qū)河流演變除了暢流期的主要沖淤過程外,封凍期亦存在一定幅度的沖淤變化,且涉及明渠流與有壓流的相互轉換,本文僅考慮了暢流期沖淤過程,并未對封凍期沖淤過程進行模擬,這同樣也會對計算結果造成一定影響,這一復雜過程的精細表達與模擬亦有待進一步深入研究。

4 結 論

基于水動力-凍脹/凍融耦合作用岸灘崩塌機理,以松花江干流大頂子山航電樞紐下游近壩段為例,建立典型季凍區(qū)河流岸灘崩塌三維水沙動力學模型,開展岸灘崩塌與河床沖淤的數值模擬分析,主要研究結論為:

(1)季凍區(qū)河流岸灘崩塌除受水沙條件影響外,因氣候影響還存在顯著凍脹/凍融特點,且對岸灘穩(wěn)定性造成明顯影響。坡腳側蝕主要由水力沖刷所控制,凍脹作用通過凍脹力會使岸灘開裂深度增加；凍融過程則會使岸坡土體力學強度降低,影響岸灘穩(wěn)定；通過力學分析,初步建立水動力-凍脹/凍融耦合作用岸灘坍塌失穩(wěn)描述模式,為季凍區(qū)河流岸灘崩塌模擬提供力學依據。

(2)將傳統(tǒng)水沙輸移和河床垂向沖淤計算與岸灘崩塌模擬相耦合,提出季凍區(qū)河流暢流期河道演變過程的數值模擬方法。岸灘崩塌對季凍區(qū)天然河流演變具有重要影響,考慮岸灘崩塌的模擬結果與實測結果沖淤定性基本一致,沖淤幅度基本相當,數學模型可較好地反映河床沖淤、岸灘崩退以及河道展寬現象。建立的季凍區(qū)河流岸灘崩塌三維數值模擬方法,為深入研究季凍區(qū)河流岸灘崩塌及其對河床演變的影響提供了一種有效技術手段。