中國北方季節(jié)性冰凍河流岸灘崩塌數(shù)值模擬
——以松花江為例

假冬冬,楊 俊,陳長英,張幸農(nóng),應(yīng) 強(qiáng)

(南京水利科學(xué)研究院水文水資源與水利工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室,江蘇 南京 210029)

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

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

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

1 水動(dòng)力-凍脹/凍融耦合作用岸灘崩塌機(jī)理

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

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

1.1 坡腳水力沖刷計(jì)算模式

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

(1)

式中:ωb為岸灘側(cè)蝕速率；λ為側(cè)蝕系數(shù),由實(shí)測(cè)資料進(jìn)行率定；γ為水體重度；γb為岸坡土體重度；u為近岸處水流流速；uc為泥沙起動(dòng)流速。

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

1.2.1 上部土體裂縫深度計(jì)算

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

(2)

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

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

(3)

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

1.2.2 岸灘崩塌力學(xué)模式

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

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

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

(4)

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

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

(5)

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

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

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

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

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

(6)

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

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

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

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

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

2.1 局部網(wǎng)格可動(dòng)技術(shù)

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

2.2 數(shù)值模擬流程

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

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

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

3.1 計(jì)算范圍及計(jì)算條件

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

圖3 計(jì)算網(wǎng)格劃分Fig.3 Computational meshes

圖4 研究河段河勢(shì)Fig.4 River regime of the study area

圖5 模型進(jìn)口流量過程Fig.5 Flow discharge processes of model inlet

圖6 計(jì)算水位—流量關(guān)系與實(shí)測(cè)水位—流量關(guān)系對(duì)比Fig.6 Calibration of relationship between water level and discharge

3.2 模型計(jì)算分析與討論

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

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

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

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

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

圖9 岸灘土體級(jí)配Fig.9 Sediment gradation of riverbanks

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

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

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

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

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

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

4 結(jié) 論

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

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

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