基于Matlab的下荊江二元岸坡崩塌過程動態(tài)模擬

陳 潔，陶桂蘭，吳俊東

(1.河海大學(xué) 港口海岸與近海工程學(xué)院，南京 210098；2.長江勘測規(guī)劃設(shè)計研究有限責(zé)任公司，武漢 430010)

下荊江河段多為上粘下砂二元結(jié)構(gòu)岸坡，且粘土層較薄，高水期下部砂土層受水流沖刷逐漸掏空、上部粘土層失去支撐后易發(fā)生崩塌破壞。特別是2003年三峽水庫蓄水運(yùn)行后，上游來沙量大幅減少，下荊江河段普遍沖刷，局部失穩(wěn)崩岸屢屢發(fā)生[1]，不僅影響堤防安全，也會給航道條件的穩(wěn)定帶來不利影響，進(jìn)而制約長江中下游“黃金航道”尺度的提升。目前關(guān)于二元結(jié)構(gòu)岸坡崩岸機(jī)理的研究已取得一定進(jìn)展：Thorne和Tovey[2]總結(jié)了二元結(jié)構(gòu)河岸可能發(fā)生剪切、繞軸及拉伸崩塌這3種崩岸形式，針對較為常見的繞軸崩塌給出了相應(yīng)河岸穩(wěn)定性計算方法。Fukuoka[3]通過野外試驗(yàn)觀察研究，提出采用上部粘土層臨界懸空寬度來預(yù)報崩塌發(fā)生。夏軍強(qiáng)[4]結(jié)合下荊江河段崩岸特征，改進(jìn)了粘土層臨界懸空寬度及河岸穩(wěn)定系數(shù)的計算公式。

關(guān)于河岸穩(wěn)定性分析模型，目前應(yīng)用較為廣泛的是從經(jīng)典Osman-Thorne理論[5]發(fā)展而來的BSTEM(bank stability and toe erosion model)模型：國外學(xué)者曾運(yùn)用該模型分析了Osage River[6]、Lost Creek[7]和Barry Fork Creek[8]等河流的河岸穩(wěn)定性問題；國內(nèi)學(xué)者王博[9]、宗全利[10]和張翼等[11]分別應(yīng)用該模型模擬了長江中游二元河岸的崩退過程。BSTEM模型主要驗(yàn)算岸坡懸臂的剪切崩塌，且在計算坡腳沖刷時假定所有崩塌土體立刻被近岸水流沖走，若要考慮崩塌土體的坡腳掩護(hù)，則需對采用該模型得到的崩岸后河岸形態(tài)進(jìn)行人為修正[10-11]，操作較為繁瑣。而根據(jù)實(shí)測資料，下荊江典型二元河岸失穩(wěn)以繞軸崩塌為主[4]，崩塌土體部分堆積坡腳，在一定時間內(nèi)掩護(hù)近岸河床，BSTEM模型有待改進(jìn)，需要針對上部粘土層的繞軸崩塌驗(yàn)算岸坡穩(wěn)定、同時考慮崩塌土體在坡腳堆積的掩護(hù)作用建立分析模型。此外，BSTEM模型基于Excel宏命令運(yùn)行，在時間上是靜態(tài)的，僅能探究特定情形或固定指標(biāo)下的河岸穩(wěn)定性，劉萬利、鄧姍姍[12-13]研究指出，除近岸水流沖刷這一主要沖積因素外，河岸形態(tài)、土體組成及河道水位變化等因素也對崩岸過程有較大影響，建立模型時有必要考慮外部因素(水位升降、土體性質(zhì)、坡腳沖淤、岸坡形態(tài)等)隨時間推移的耦合變化使河岸狀態(tài)發(fā)生的改變。

為此，本文針對下荊江河段概化的典型二元結(jié)構(gòu)岸坡，考慮了水位、土體性質(zhì)、坡腳沖淤、岸坡形態(tài)等因素隨時間的耦合變化，利用Matlab編程模擬了水位升降過程中的下部砂土層淘刷、上部粘土層繞軸崩塌、以及崩塌土體堆積輸移這3個循環(huán)階段[3]，建立了下荊江岸坡崩塌動態(tài)分析模型。在此基礎(chǔ)上，重點(diǎn)探討了不同水位變化模式對河岸崩退過程的影響，為二元結(jié)構(gòu)岸坡穩(wěn)定性及其崩岸規(guī)律研究提供參考。

1 岸坡崩塌模型

本文改進(jìn)現(xiàn)有岸坡崩塌的臨界判斷條件，推導(dǎo)崩塌土體在坡腳堆積寬度公式，建立起二元岸坡崩塌動態(tài)分析模型。該模型以計算水流沖刷導(dǎo)致的岸坡形態(tài)改變?yōu)榛A(chǔ)，針對上部粘土層出露懸臂的繞軸崩塌驗(yàn)算岸坡穩(wěn)定性，同時考慮崩塌土體在坡腳的堆積掩護(hù)與輸移，模擬了岸坡崩退過程中循環(huán)發(fā)生的3個主要階段。

1.1 岸坡橫向沖刷距離的計算

水流對河岸土體的橫向沖刷距離ΔB根據(jù)剩余切應(yīng)力法計算[13]

ΔB=K·Δt·(τ0-τc)

(1)

(2)

τ0=γwRJ

(3)

式中：K為土體侵蝕系數(shù)，m3/N·s，與土體自身性質(zhì)有關(guān)，本文采用Hanson等[14]提出同時適用于粘性土和砂性土的經(jīng)驗(yàn)公式計算，見式(2)；Δt為沖刷時間，s；τc為河岸土體起動切應(yīng)力，N/m2，一般通過土體沖刷試驗(yàn)或根據(jù)粒徑大小按照經(jīng)驗(yàn)公式計算得到；τ0為水流切應(yīng)力，N/m2，與水深成正比；γw為水體容重，N/m3；R為水力半徑，m；J為水面比降，m/m。

對于上粘下砂二元結(jié)構(gòu)岸坡，當(dāng)水流切應(yīng)力超過下部砂土層的起動切應(yīng)力后，該土層將逐漸被近岸水流掏空，使得上部粘土層出露懸空，本文根據(jù)每個分析時段沖刷后的岸坡斷面形態(tài)確定粘土層實(shí)際懸空寬度。

1.2 岸坡崩塌臨界條件的判斷

圖1 上部粘性土層懸空土塊的受力分析Fig.1 Stress analysis of upper cohesive cantilever of composite riverbank

上部粘土層的懸空寬度隨水流沖刷逐漸增大，當(dāng)超過某一臨界值時，自身重力矩大于土體抵抗力矩，使其繞某一中性軸產(chǎn)生向河槽方向的旋轉(zhuǎn)運(yùn)動倒入水中，即繞軸崩塌，如圖1所示。

根據(jù)懸臂梁平衡原理，當(dāng)粘性懸空土塊處于臨界狀態(tài)時，其自重W引起的外力矩與土體潛在斷裂面上產(chǎn)生的抵抗力矩相平衡，即[4]

(4)

式中：Bc為粘土層臨界懸空寬度，m；H1、Ht分別為粘土層厚度和其頂部張拉裂隙深度，m；a為粘性土的抗拉、抗壓應(yīng)力之比，即a=σt/σc；σt、σc分別為土體的抗拉、抗壓強(qiáng)度，kPa；W為單位長度懸空土塊自重，kN。

(5)

(6)

式中：h1、h2分別為計算水位以上及以下的粘土層厚度，m，滿足h1+h2=H1；γ′為水位以下的粘性土容重，取其浮容重；γ為水位以上的粘性土容重，考慮到粘性土的低滲透性，落水期土體內(nèi)部水分來不及排出，取其飽和容重；Ht為粘性土頂部張拉裂隙深度，m，本文采用Terzaghi提出的經(jīng)驗(yàn)公式[15]估算，見式(6)；c、φ分別為粘性土的粘聚力，kPa，摩擦角(°)；其他符號含義同上。根據(jù)上部粘土層懸空寬度實(shí)際值B與臨界值Bc之間的相對關(guān)系，可以判斷其是否發(fā)生繞軸崩塌：若B>Bc，則表示粘性懸空土塊崩塌。

1.3 崩塌土體在坡腳堆積寬度的計算

圖2 崩塌土體在坡腳堆積掩護(hù)示意圖Fig.2 Sketch of covering of collapsed soil at bank toe

粘性懸空土塊崩塌后將在坡腳堆積掩護(hù)，隨后在水流作用下逐漸分解并輸移，如圖2所示，本文僅考慮計算斷面崩塌土體在坡腳堆積的掩護(hù)作用。

根據(jù)下荊江二元河岸崩塌過程的概化水槽試驗(yàn)結(jié)果[16]，崩塌土體在坡腳呈三角形堆積，堆積體積占崩塌體積一定比例，堆積坡度近似等于河岸水下穩(wěn)定坡比或泥沙水下休止角。本文據(jù)此確定堆積體積V和堆積坡度α，則崩塌土體在坡腳堆積寬度Lc和堆積高度hc滿足關(guān)系式

Lc=(cotα-cotβ)hc

(7)

(8)

V=KbBcH1

(9)

式中：Lc、hc分別為崩塌土體堆積的寬度和高度，m；α和β分別為坡腳堆積土體坡度和岸坡坡度，°；V為堆積土體體積，m3；Kb為堆積系數(shù)，即堆積土體占崩塌土體的體積比例，可根據(jù)水流條件確定：流速越大，水流對土體沖刷作用越強(qiáng)，Kb值越?。籅c、H1分別為崩塌粘土塊的寬度和厚度，m。根據(jù)式(8)和式(9)求得堆積高度hc，將hc代入式(7)得到崩塌土體在坡腳堆積寬度Lc的計算公式如下

(10)

由于下部砂土層抗沖性較差、易起動，而粘性土受顆粒間粘結(jié)力影響具有較大的起動切應(yīng)力，因此粘性崩塌土體一定時間內(nèi)的掩護(hù)作用使得坡腳抗沖性增強(qiáng)、側(cè)蝕速率減小[16]。本文分析時首先判斷上部粘土層崩塌情況，若岸坡發(fā)生崩塌，根據(jù)式(10)計算崩塌土體堆積寬度Lc，下一時段水流首先沖刷坡腳堆積的粘性土，沖刷距離根據(jù)式(1)～式(3)代入粘性土參數(shù)求得，待粘性土覆蓋層全部被水流沖散輸移，水流繼續(xù)沖刷坡腳砂性土。

2 基于Matlab軟件的岸坡崩塌過程模擬

2.1 基于Matlab的岸坡崩塌動態(tài)模擬程序

圖3 岸坡崩塌模型計算流程Fig.3 Calculation process of bank slope collapse model

Matlab軟件是美國MathWorks公司出品的商業(yè)數(shù)學(xué)軟件，擁有強(qiáng)大的算法開發(fā)、數(shù)值計算、數(shù)據(jù)分析及數(shù)據(jù)可視化功能，本文采用該軟件編制二元岸坡崩退過程的動態(tài)模擬程序，計算流程見圖3。

由此建立的岸坡崩塌模型具有以下特點(diǎn)：(1)高效的數(shù)值計算功能。將計算總時段細(xì)分為若干分析時段，利用循環(huán)語句模擬水位、土體性質(zhì)、坡腳沖淤、岸坡形態(tài)等因素隨時間的耦合變化；(2)動態(tài)的圖形處理功能。利用動態(tài)繪圖實(shí)現(xiàn)岸坡崩退過程中的斷面形態(tài)變化展示；(3)便捷的輸入輸出界面。輸入不同初始條件即可應(yīng)用該程序模擬不同岸坡的崩塌過程，計算結(jié)果便于匯總輸出，實(shí)現(xiàn)與有限元等其他軟件的數(shù)據(jù)交換。

2.2 二元岸坡崩塌模型驗(yàn)證

為了驗(yàn)證岸坡崩塌模型的計算可靠性，本文采用所編寫的程序模擬了下荊江荊98斷面岸坡[11]在2006～2007水文年內(nèi)的崩退過程，相關(guān)輸入?yún)?shù)與現(xiàn)場實(shí)測數(shù)據(jù)見文獻(xiàn)[11]，天然河岸水下穩(wěn)定坡比取為0.323[17]，模擬所得崩后斷面與BSTEM模型計算結(jié)果以及現(xiàn)場實(shí)測數(shù)據(jù)的對比如圖4所示。由圖可見，本文模型計算得到的岸頂崩退總寬度(70 m)與實(shí)際情況(72 m)較為符合，且模擬的崩后河岸形態(tài)與實(shí)際形態(tài)的差別也相對較小。

圖4 荊98斷面岸坡崩后斷面 計算結(jié)果對比Fig.4 Comparison of cross-section after bank collapse at Jing 98

3 下荊江二元岸坡崩岸過程模擬分析

在岸坡崩塌模型驗(yàn)證可靠的基礎(chǔ)上，本文應(yīng)用該模型重點(diǎn)探討了不同水位變化模式對二元河岸崩退過程、崩塌形態(tài)變化及坡腳沖刷距離的影響。

3.1 河段概況及典型斷面概化

根據(jù)相關(guān)地質(zhì)資料[18]，下荊江兩岸主要為河流沖積平原，地勢平緩，地形地貌及地質(zhì)結(jié)構(gòu)相近。因此，本文匯總下荊江(藕池口—城陵磯段)岸坡斷面、土體參數(shù)等數(shù)據(jù)，按各河段長度加權(quán)計算其平均值作為概化岸坡參數(shù)，以使研究對象更具普遍性與代表性，岸坡土體參數(shù)見表1。二元河岸的岸坡形態(tài)與上部粘土層厚度有重要關(guān)系：當(dāng)粘土層較薄時，岸坡坡度上陡下緩，本文將粘土層坡度簡化為垂直，概化岸坡斷面如圖5所示。

圖5 二元結(jié)構(gòu)概化岸坡斷面Fig.5 Generalized cross-section of composite river banks

土體參數(shù)天然容重/kPa含水率/%孔隙比粘聚力/kPa內(nèi)摩擦角/°粘性土18.6320.9221.911砂性土20.2220.58028.8

3.2 計算參數(shù)與模擬工況

影響岸坡沖刷的因素主要有土體侵蝕系數(shù)、起動切應(yīng)力、水力半徑和水面比降等。水力半徑R近似為岸坡水深；水面比降J根據(jù)岸坡斷面上、下游水位差求得，漲水期取0.008、落水期取0.012。土體參數(shù)取值如表2所示[18]。

表2 土體計算參數(shù)[18]Tab.2 Soil calculation parameters[18]

表3 水位變化模式Tab.3 Modes of water level variation

影響粘土層崩塌的因素主要有粘性土的粘聚力、摩擦角、容重和抗拉強(qiáng)度，前三者取值見表1。粘性土抗拉強(qiáng)度σt可根據(jù)其抗剪強(qiáng)度估算得到，并隨含水率增加呈先增加后減小的規(guī)律[4]，高、低水位下分別將其取為4 kPa和7 kPa。根據(jù)下荊江河岸崩塌概化試驗(yàn)[16]，高、低水期分別取崩塌土體堆積系數(shù)Kb為0.27和0.58，取其水下堆積坡比tanα=1/2.5。長江監(jiān)利站2012～2014年的水位資料表明[19]，長江水位最高、最低值分別為34 m和24 m左右，漲水速率和退水速率相近，約為0.1 m/d。據(jù)此，本文選取計算工況如表3所示。

3.3 二元結(jié)構(gòu)岸坡崩退過程模擬

兩種水位變化模式下，二元結(jié)構(gòu)岸坡崩退形態(tài)變化如圖6、圖7所示，其中帶五角星符號的斷面線代表兩工況的第一次崩后形態(tài)，其余帶相同符號的斷面線分別對應(yīng)于相近沖刷時長后的岸坡形態(tài)，矩形框內(nèi)為坡腳的局部放大圖。

對于處于水位下降模式的工況1，水位下降初期岸坡受水流沖刷的范圍較廣且沖刷量較大，上部粘土層懸空、不斷發(fā)生繞軸崩塌。隨著水位繼續(xù)下降，岸坡橫向沖刷范圍及沖刷量均隨之減小，高程大的岸坡點(diǎn)逐漸露出水面不受沖刷，上部粘土層保持穩(wěn)定不再崩塌。對于處于水位上升模式的工況2，其岸坡崩退形態(tài)的變化規(guī)律與工況1恰好相反：隨著水位逐步上升，岸坡橫向沖刷范圍及沖刷量均由小變大，高程大的岸坡點(diǎn)逐漸開始受水流沖刷，上部粘土層進(jìn)入懸空狀態(tài)發(fā)生繞軸崩塌。兩工況下岸坡的具體崩塌情況匯總?cè)绫?所示。

圖6 水位下降期岸坡崩退過程Fig.6 Collapse process of bank slope under water level decline圖7 水位上升期岸坡崩退過程Fig.7 Collapse process of bank slope under water level rise

表4 兩種水位變化模式下岸坡崩塌情況Tab.4 Slope collapse under two modes of water level variation

結(jié)合表4和圖6、圖7中兩工況的最終崩后斷面對比可見，水位下降過程中岸坡的崩塌總次數(shù)、崩退總寬度均大于水位上升過程，其主要原因和岸坡沖刷幅度與土體性質(zhì)有關(guān)：一方面，水位下降前高水位對河岸土體的浸泡降低了其抗拉強(qiáng)度、土體斷裂面產(chǎn)生的抵抗力矩較小；另一方面，相同水位下，水位下降期較大的水面比降[18]使得其岸坡沖刷速率略大于水位上升期，由此導(dǎo)致該時期的粘土層實(shí)際懸空寬度較大，同時由于粘性土的低滲透性，土體內(nèi)部水分來不及排出，飽和懸空土塊產(chǎn)生的自重力矩較大。綜上，較小的有利抵抗力矩與較大的不利自重力矩，使得水位下降模式下岸坡上部粘土層更易發(fā)生繞軸崩塌。

根據(jù)式(3)，岸坡沖刷速率與水位成正比，由此水位下降初期粘土層的實(shí)際懸空寬度累加最快，但高水位時的水體浮托力有利于其穩(wěn)定性；隨著水位下降，粘土層臨界懸空寬度隨浮托力的逐漸消失而減小，較大的岸坡沖刷速率使得其實(shí)際懸空寬度仍較快累加，最終導(dǎo)致水位變化至兩土層交界面附近時(33～31m)最易發(fā)生繞軸崩塌。同理對于漲水期，中高水位時(31～33 m)粘性土懸空寬度較大的實(shí)際值與較小的臨界值使其成為崩岸的最危險時期。

雖然水位下降期的坡腳沖刷速率稍大于水位上升期，但兩工況下的坡腳最終累計沖刷距離相近，分別達(dá)5.44 m和5.56 m，這與崩塌土體的掩護(hù)作用有關(guān)：工況1在水位下降初期(33.4 m)就發(fā)生上部粘土層繞軸崩塌，崩塌土體堆積掩護(hù)，減緩坡腳沖刷速率。隨后不斷發(fā)生岸坡崩塌、崩塌土體堆積，水流持續(xù)以較小速率沖刷坡腳粘性土。粘性土掩護(hù)層在水位下降末期(24.7 m)被完全沖散，水流重新以較大速率沖刷坡腳砂性土；而工況2在水位上升中后期(31.6 m)才開始發(fā)生岸坡崩塌、受崩塌土體掩護(hù)的影響減緩坡腳沖刷速率，由此導(dǎo)致其最終累計沖刷距離接近工況1。與未考慮坡腳掩護(hù)相比，落、漲水期的坡腳累計沖刷距離分別減少約45%和17%，說明崩塌土體堆積坡腳的掩護(hù)作用一定程度上影響了岸坡崩退過程，崩塌次數(shù)越多、影響程度越大。

4 結(jié)論

本文結(jié)合下部砂土層沖刷計算、上部粘土層穩(wěn)定性分析與推導(dǎo)的崩塌土體坡腳堆積公式，利用Matlab編程建立了二元岸坡崩塌分析模型。在此基礎(chǔ)上，針對下荊江河段典型概化斷面，重點(diǎn)探討了不同水位變化模式對二元河岸崩退過程的影響，主要結(jié)論如下：

(1)應(yīng)用本文所建立的岸坡崩塌模型動態(tài)模擬下荊江二元河岸崩退過程，與現(xiàn)場實(shí)測符合良好，且具有計算高效、模擬多因素耦合變化、動態(tài)展示斷面形態(tài)、輸入輸出便捷等優(yōu)勢。

(2)水位下降期較之水位上升期上部粘土層更易發(fā)生繞軸崩塌：一方面，水位下降時較大的岸坡沖刷速率導(dǎo)致粘土層懸空寬度累加較快，同時由于土體內(nèi)部水分來不及排出，飽和懸空土塊產(chǎn)生的自重力矩較大；另一方面，落水期土體較小的抗拉強(qiáng)度使得潛在斷裂面產(chǎn)生的抵抗力矩較小。

(3)考慮河道水體浮托力隨水位變化后，該浮托力對粘性懸空土塊穩(wěn)定性的有利影響主要體現(xiàn)在高水期，并隨水位下降消失，水位變化至兩土層交界面附近時是岸坡崩塌的最危險時期。

(4)與未考慮坡腳掩護(hù)相比，考慮該掩護(hù)作用后落、漲水期的坡腳累計沖刷距離分別減少約45%和17%，岸坡崩塌次數(shù)越多、崩塌土體坡腳掩護(hù)對岸坡崩退過程的影響越大。