彎曲河道是天然河道中常見的一種河流形態(tài)，多年來人們從各個(gè)方面對(duì)彎曲河道特有的水流運(yùn)動(dòng)規(guī)律、河床演變規(guī)律，進(jìn)行了廣泛的研究[1，2].如張紅武[2]從模型試驗(yàn)和理論分析出發(fā)，較為深入地研究了彎道內(nèi)的水面形態(tài)、環(huán)流運(yùn)動(dòng)以及縱向流速的沿程分布規(guī)律等.隨著計(jì)算機(jī)性能的提高，以及數(shù)值計(jì)算方法的發(fā)展，很多研究者開始用數(shù)學(xué)方法模擬彎道內(nèi)的水沙運(yùn)動(dòng)及河床變形.最早以Engelund[3]、Odgaard[4]等人提出的恒定水流模型為代表.近年來，部分研究者提出的三維數(shù)學(xué)模型，既能模擬彎道內(nèi)的水沙運(yùn)動(dòng)，又能考慮床面的沖淤變化，較以前有了更大的進(jìn)步[5，6].盡管三維數(shù)學(xué)模型可以模擬彎道內(nèi)較為復(fù)雜的水流運(yùn)動(dòng)和河床變形問題，不過在實(shí)際工程中用的最廣的還是平面二維模型[7，8].

通常現(xiàn)有這些模型只考慮河床的縱向變形，而很少注意彎道內(nèi)的橫向變形問題，也即河岸沖刷問題.本文針對(duì)當(dāng)前在彎道水流和河床變形模型的研究狀況，以前人的研究結(jié)果為基礎(chǔ)，建立了正交曲線坐標(biāo)下的平面二維水沙數(shù)學(xué)模型，用于模擬彎道內(nèi)的水沙運(yùn)動(dòng)及河床縱向變形；在深入分析河岸沖刷機(jī)理的基礎(chǔ)上，采用力學(xué)模型模擬粘性河岸沖刷、崩塌的過程.最后利用彎道模型的試驗(yàn)資料，驗(yàn)證了本模型計(jì)算的流場(chǎng)、岸邊水位等與實(shí)測(cè)值符合較好.同時(shí)模擬了90°彎道在持續(xù)清水沖刷下的河床變形過程，對(duì)主流線、斷面形態(tài)、主槽比降等計(jì)算結(jié)果隨時(shí)間變化的分析表明，模擬結(jié)果符合彎道變形規(guī)律.

1 數(shù)學(xué)模型

　　本文建立的模型具有如下特點(diǎn)：(1)為適應(yīng)不規(guī)則的岸邊界，建立正交曲線坐標(biāo)下的平面二維水沙模型，并用MADI方法[9]求解水流方程組.(2)改進(jìn)以往對(duì)泥沙恢復(fù)飽和系數(shù)取為經(jīng)驗(yàn)常數(shù)的做法，模型中采用了張紅武[10]提出的不平衡輸沙理論，用于研究懸移質(zhì)泥沙的輸移以及河床變形過程.(3)引入Osman[11]提出的河岸沖刷模型，用于模擬粘性河岸的橫向沖刷過程以及在重力作用下的崩塌過程.

1.1 平面二維水沙數(shù)學(xué)模型

1.1.1水流基本方程

　　式中：ε_ξ、ε_η分別表示ξ、η方向的泥沙擴(kuò)散系數(shù)；S_0k、S_k、S_*k、ω_k分別為第k粒徑組泥沙的側(cè)向輸入項(xiàng)、分組含沙量、分組挾沙力及有效沉速；α_*、f₁、K₁分別為平衡含沙量分布系數(shù)，非飽和系數(shù)以及附加系數(shù).

　　上述3個(gè)參數(shù)以及水流挾沙力的計(jì)算方法，詳見文獻(xiàn)[10].而分組挾沙力級(jí)配以及床沙級(jí)配的計(jì)算方法，見文獻(xiàn)[12].

1.1.3河床縱向變形方程 由懸移質(zhì)泥沙不平衡輸移引起的河床縱向變形方程為：

式中：ρ＇為床沙干密度，N為非均勻泥沙的分組數(shù).

1.2 幾個(gè)關(guān)鍵問題的處理

1.2.1正交曲線網(wǎng)格的生成 設(shè)(x,y)為物理平面上的笛卡爾坐標(biāo)系，(ξ,η)為相應(yīng)計(jì)算平面上的變換坐標(biāo)系，用以下方程可實(shí)現(xiàn)兩個(gè)坐標(biāo)系之間的變換，即：

　　求解上述方程組，即可得到物理平面(x,y)與計(jì)算平面(ξ,η)上對(duì)應(yīng)點(diǎn)的關(guān)系.當(dāng)網(wǎng)格正交時(shí)，β=0.控制函數(shù)的具體表達(dá)式，參見文獻(xiàn)[13].

1.2.2初、邊界條件 對(duì)初始條件，一般給出水位、流速、含沙量等變量的初值.對(duì)進(jìn)口邊界、給定流量、沙量過程線，以及含沙量級(jí)配.對(duì)出口邊界、給定水位過程線，并認(rèn)為S/ξ=0.對(duì)岸邊界，對(duì)水流可采用滑移或無滑移邊界條件；對(duì)泥沙，可取S/=0(為岸邊界法線方向).

　　在本文模擬的彎道中，橫斷面由主槽和灘地組成.假設(shè)灘地不過水，則岸邊界是灘岸，同時(shí)認(rèn)為灘岸可以沖刷.實(shí)際計(jì)算區(qū)域包括灘地和主槽，故采用“凍結(jié)法”技術(shù)[14]處理灘地上那些不過水的節(jié)點(diǎn).一般岸坡坡頂和坡腳間的水平距離較小，可能遠(yuǎn)小于網(wǎng)格尺寸，在程序中不能分辨.為了能模擬河岸的沖刷過程，同時(shí)又不增加計(jì)算量，故本文在計(jì)算中用一數(shù)組記錄岸邊相鄰兩節(jié)點(diǎn)之間的實(shí)際地形.

1.3 河岸沖刷模型 現(xiàn)有的很多泥沙數(shù)學(xué)模型，很少考慮到河岸沖刷問題.即使那些考慮了河岸沖刷的模型，往往對(duì)彎道中河岸沖刷機(jī)理、沖刷速率，做了大量的假設(shè)和簡化.如Ikeda[15]認(rèn)為彎道凹岸附近的縱向垂線平均流速大于彎道中心處的垂線平均流速時(shí)，彎道凹岸沖刷，反之則淤積.Hasegawa[16]認(rèn)為河岸沖刷速率與近岸的剩余流速成正比，根據(jù)泥沙連續(xù)方程，得出了適用于彎曲河道通用的河岸沖刷系數(shù).不過，這類模型中存在一個(gè)共同缺點(diǎn)，即很少考慮到河岸組成物質(zhì)、土體特性、幾何形態(tài)等對(duì)河岸沖刷的影響.眾所周知，粘性河岸和非粘性河岸的沖刷機(jī)理和沖刷速率是完全不同的.

　　為準(zhǔn)確模擬河岸的橫向沖刷過程，必須深入分析河岸沖刷機(jī)理.一般來講，河岸沖刷后退是河岸土體和近岸水流相互作用的結(jié)果.除了岸邊植被生長情況、河道內(nèi)水位升降、滲流、管涌等因素影響河岸后退外，但主要是以下2種情況，是導(dǎo)致河岸后退的主要原因：(1)通過水流直接橫向沖刷河岸導(dǎo)致河岸后退.近岸水流直接作用于河岸，沖動(dòng)河岸邊坡上水面以下的表層土體，并被水流帶走，從而導(dǎo)致河岸后退.(2)通過河岸崩塌導(dǎo)致河岸后退.崩塌是河岸上的一部分土體在重力作用下，沿某一滑動(dòng)面發(fā)生移動(dòng)的過程.一般床面發(fā)生沖刷，導(dǎo)致河岸高度增加，或者水流淘刷河岸坡腳，使河岸坡度變陡，都會(huì)降低河岸的穩(wěn)定性，當(dāng)河岸的穩(wěn)定性降低到一定程度后，河岸便會(huì)發(fā)生崩塌，導(dǎo)致岸邊界后退.

圖1 粘性河岸沖刷的計(jì)算模式

　　本文主要考慮粘性河岸的沖刷后退過程，采用Osman[11]提出的河岸沖刷模型(圖1).該模型首先計(jì)算河岸橫向沖刷距離，然后分析河岸是否會(huì)失穩(wěn)、崩塌.在Δt(sec)時(shí)間內(nèi)，粘性河岸被水流橫向沖刷后退的距離為：

式中：γs河岸土體的容重(kN/m3)；ΔB為Δt時(shí)間內(nèi)河岸因水流橫向沖刷而后退的距離(m)；τ為作用在河岸上的水流切應(yīng)力(N/m2)；τc為河岸土體的起動(dòng)切應(yīng)力(N/m2).Cl為橫向沖刷系數(shù)，取決于河岸土體的物理化學(xué)特性.

　　當(dāng)由式(7)得河槽沖寬ΔB，用平面二維水沙模型算出河床沖深ΔZ后，河岸高度增加，坡度變陡，穩(wěn)定性降低.根據(jù)土力學(xué)中的邊坡穩(wěn)定性關(guān)系，采用若干假定，可得到河岸發(fā)生初次崩塌時(shí)的臨界條件.若河岸已發(fā)生初次崩塌，則假定以后的河岸崩塌方式為平行后退，即崩塌后的邊坡角度恒為β，仍可用土力學(xué)的方法判斷是否會(huì)發(fā)生二次崩塌.

1.4 計(jì)算方法和求解過程 首先，采用MADI法[9]計(jì)算流場(chǎng).這種方法采用非交錯(cuò)網(wǎng)格，將水位、流速等變量均布置在同一網(wǎng)格點(diǎn)上，對(duì)水流連續(xù)方程和動(dòng)量方程均不按方向剖開，由此將基本運(yùn)動(dòng)方程離散而形成新的差分代數(shù)方程組，并建立一種新的解法.具體求解過程如下：對(duì)式(1)進(jìn)行差分離散，時(shí)間導(dǎo)數(shù)項(xiàng)采用前差表示，空間導(dǎo)數(shù)項(xiàng)采用中心差分表示，將時(shí)間步長Δt分成兩部分.在前半時(shí)間步長內(nèi)，利用連續(xù)方程(1)、ξ方向動(dòng)量方程(2)，沿ξ方向采用隱式離散，對(duì)η方向采用顯式離散，求解得Zn+1/2,Un+1/2，再利用η方向動(dòng)量方程(3)，顯式求解得Vn+1/2.在后半時(shí)間步長內(nèi)，利用式(1)、(3)，沿η方向采用隱式離散，對(duì)ξ方向采用顯式離散，求解得Zn+1，Vn+1，再利用式(2)，顯式求解得Un+1.由于要進(jìn)行長時(shí)間的河床變形計(jì)算，計(jì)算量特別大，在實(shí)際計(jì)算中通常用偽恒定方法計(jì)算出恒定流場(chǎng).然后，采用顯隱混合格式計(jì)算懸移質(zhì)含沙量分布.先按不同方向ξ、η對(duì)式(4)進(jìn)行算子分裂，在前半步長內(nèi)，對(duì)ξ方向的分步方程用指數(shù)格式顯式離散求得Sn+1/2，在后半步長內(nèi)，對(duì)η方向的分步方程用Crank-Nicholson型格式隱式離散求得Sn+1.以時(shí)間為迭代參數(shù)，計(jì)算出恒定的濃度場(chǎng).

　　最后，顯式求解式(5)，可得出時(shí)段末的床面高程.在已知河床變形的基礎(chǔ)上，采用河岸沖刷模型，模擬彎曲河道河岸沖刷后退的過程，具體計(jì)算方法見文獻(xiàn)〔12〕.此外，本文在計(jì)算中還采用以下假定：直接從河岸沖刷下來的和河岸崩塌后產(chǎn)生的泥沙，全部轉(zhuǎn)化為懸移質(zhì)泥沙，并作為下一個(gè)時(shí)段的側(cè)向輸沙量.

2模型的驗(yàn)證

　　為了準(zhǔn)確模擬彎道的河床變形過程，本文首先對(duì)水流模型進(jìn)行驗(yàn)證.驗(yàn)證的資料來自羅索夫斯基[1]的180°彎道的水槽試驗(yàn)資料.該水槽由一6m長的順直進(jìn)口段、180°的彎道段以及3m長的出口段組成.水槽寬度為0.80m，過水?dāng)嗝鏋榫匦?彎道外半徑為1.2m，內(nèi)半徑為0.4m.計(jì)算中的水流條件如下：流量為12.3×10-3m3/s，進(jìn)口斷面選在水槽進(jìn)口以下2m處；出口斷面選在距水槽出口0.1m處，斷面水深為0.054m.水槽糙率取0.012，計(jì)算區(qū)域內(nèi)共劃分網(wǎng)格數(shù)為95×11.

圖2 彎道凹岸與凸岸岸邊水位的沿程分布

圖3 縱向垂線平均流速分布

　　圖2給出了彎道中凹岸和凸岸水位沿程變化情況.從圖中可以看出，在彎道進(jìn)口段，凹岸和凸岸的水位基本相同；在彎道段，凹岸水位明顯高于凸岸；在彎道出口段，兩岸水位基本相同.盡管計(jì)算值略小于實(shí)測(cè)值，但兩者的變化趨勢(shì)相當(dāng)符合，當(dāng)水流由順直河段進(jìn)入彎道后，由于受到離心力的作用，使彎道凹岸一側(cè)的水位恒高于凸岸一側(cè)，最大水位差一般出現(xiàn)在彎道頂點(diǎn)附近，而向上下游兩個(gè)方向逐漸減少.由于未能考慮彎道內(nèi)橫向環(huán)流對(duì)水流動(dòng)量方程的影響，從而導(dǎo)致水位計(jì)算值總體偏小.

　　圖3給出了縱向垂線平均流速的沿程分布情況.從圖可知，實(shí)測(cè)值與計(jì)算值符合較好，尤其在彎道的進(jìn)口段附近.從圖中還可以看出，在彎道的進(jìn)口段，最大縱向平均流速的位置(主流線)緊靠凸岸；在彎道段，主流線仍緊靠凸岸；在彎道出口段附近，主流線逐漸向凹岸轉(zhuǎn)移.出彎后的水流，在相當(dāng)長的距離內(nèi)，最大流速仍靠近外側(cè)(凹岸).

3 計(jì)算結(jié)果與分析

　　由于目前缺少由懸移質(zhì)泥沙不平衡輸移的實(shí)測(cè)資料，故本文采用一概化的90°彎道，研究其在持續(xù)清水沖刷、兩岸邊界可動(dòng)的條件下，彎道內(nèi)流速、斷面形態(tài)，以及比降的變化過程.對(duì)模擬結(jié)果的分析表明：本文所建立的平面二維水沙數(shù)學(xué)模型以及河岸沖刷模型，能較好的反映彎曲河道的演變過程.

3.1 河道概況 假定有一概化彎道，包括順直進(jìn)口段(1000m)、90°彎道段(2100m)以及順直出口段(2000m).河床初始縱比降為0.0002，初始糙率為0.02，床沙平均粒徑為0.087mm.河道主槽寬度約450m，兩岸灘地總寬約250m；初始河岸坡度和高度分別為60°、5.0m(見圖5).河岸土體特性為：容重γs=18kN/m3、內(nèi)摩擦角φ=14°、凝聚力C=20kN/m2、τc=1.2N/m2.模擬河段共劃分網(wǎng)格為50×28，每個(gè)計(jì)算時(shí)段恒為2d.假設(shè)該彎道進(jìn)口處的懸移質(zhì)含沙量很小，計(jì)算中取0.10kg/m3，而流量始終為3000.0m3/s；彎道出口斷面初始水位為104.0m，在水流的持續(xù)沖刷下，水位不斷下降，在計(jì)算中假定水位下降幅度為河床沖刷厚度的0.2倍.

3.2 流速分布 圖4為初始時(shí)刻的流速分布圖.在初始時(shí)刻，主流線在進(jìn)口處靠近凹岸，當(dāng)水流進(jìn)入彎道后，由于離心力的作用，主流線逐漸貼近彎道凸岸，在彎道出口處，主流線又逐漸偏離凸岸，而向凹岸靠近.當(dāng)該彎道經(jīng)過30d的持續(xù)清水沖刷后，平均水深3.9m增大4.7m，平均流速由1.6m/s減小1.2m/s.不僅彎道內(nèi)的水流條件發(fā)生變化，而且斷面形態(tài)、平面形態(tài)也發(fā)生變化，從而使彎道內(nèi)的主流線位置也發(fā)生了變化.從圖5中可以看出，在彎道進(jìn)口段，主流線基本上貼近凹岸；在彎道段，主流線仍貼近凸岸；在彎道出口段，主流線逐漸向河道中心位置轉(zhuǎn)移.彎道內(nèi)水深增加，流速減小，從而導(dǎo)致主流頂沖凹岸的位置逐漸上移.主流線頂沖凹岸的部位，符合小水上提，大水下挫的彎道水流運(yùn)動(dòng)規(guī)律.此外，從圖5中還可看出，由于河岸沖刷、崩塌后退，原先的那些不過水的灘地位置，逐漸開始過水.

　　由此可見，彎道內(nèi)垂線平均流速的分布，不僅與進(jìn)口處的流速分布、出口斷面的水位有關(guān)，而且還與彎道的斷面形態(tài)、平面形態(tài)等因素有關(guān).

圖4 初始流速分布

圖5 30d后的流速分布

圖6 斷面形態(tài)變化過程

圖7 主槽平均高程變化過程

3.3斷面形態(tài)變化 圖5給出了彎道段第21號(hào)斷面的斷面形態(tài)隨時(shí)間變化過程.從變化圖中可以看出，由于兩岸灘地可以沖刷，主槽在不斷沖刷下切的同時(shí)，河床逐漸展寬.由于在彎道段，主流線的位置一直靠近凸岸附近，導(dǎo)致凸岸附近的流速大于凹岸附近流速，因而凸岸的灘地沖刷、崩塌后退的距離(28.6m)略大于凹岸灘地后退距離(28.4m).而且隨著河床逐漸展寬，灘地沖刷、后退的速率逐漸減小.另外，30d內(nèi)河床縱向總沖深為1.37m，每隔10d的沖刷深度分別為0.61m、0.43m、0.33m.由此可見，在進(jìn)口水沙條件不變的條件下，河床的縱向沖刷速率隨時(shí)間增加而逐漸減小.

　　彎道進(jìn)口段和出口段的斷面形態(tài)變化，也有類似的變化規(guī)律.不過，應(yīng)當(dāng)指出，這里未出現(xiàn)一般彎道演變中的凹岸沖刷，凸岸淤積的現(xiàn)象，主要由以下兩方面的原因：一是河床受清水沖刷，凸岸附近的含沙量(平均為0.27kg/m3)遠(yuǎn)小于當(dāng)?shù)氐乃�流挾沙�?平均為3.2kg/m3)，因此發(fā)生沖刷.二是由于水流中挾帶的懸移質(zhì)含沙量本身就很少，由橫向環(huán)流從凹岸帶到凸岸的懸沙數(shù)量也少，加上凸岸沒有高程較低的邊灘，因而凸岸不發(fā)生淤積.

3.4主槽比降的變化 圖6給出了河岸可沖刷情況下，主槽平均高程隨時(shí)間的變化過程.在清水的持續(xù)沖刷下，整個(gè)河床都有不同程度的沖深，30d內(nèi)，進(jìn)口處沖深較多(1.3m)，出口處沖深相對(duì)較小(1.0m).主河槽縱比降初始值為2×10-4，河床沖刷10d、20d、30d后的比降變?yōu)?.82×10-4、1.63×10-4、1.59×10-4.隨著河床的逐漸沖深和橫向展寬，比降變小的速度降低.本文還計(jì)算了河岸不可沖刷情況下，主槽平均高程隨時(shí)間的變化過程.在河岸不可沖刷時(shí)，河床沖深明顯加大.在30d內(nèi)進(jìn)口處沖深約1.5m，出口處約1.2m.河床縱比降由原來的2×10-4下降為1.4×10-4，下降幅度為30%.

　　從上述分析來看，在河岸可沖刷的情況下，由兩個(gè)重要原因可以減緩河床沖刷：一是河槽寬度增大，導(dǎo)致過水面積增大，流速較小，水流挾沙力減小，從而減少對(duì)床面的沖刷.二是從河岸沖刷和崩塌下來的泥沙數(shù)量部分滿足了懸移質(zhì)水流挾沙力的要求，因而可減少從床面沖起的泥沙數(shù)量.此外，還有河床沖刷過程中產(chǎn)生的床沙粗化現(xiàn)象，可引起河道阻力增加，也可減緩對(duì)河道的沖刷.

4 結(jié)論

　　本文在分析前人對(duì)彎道內(nèi)水流運(yùn)動(dòng)和河床變形的基礎(chǔ)上，針對(duì)當(dāng)前很多模型中較少考慮河岸沖刷問題的狀況，提出了正交曲線坐標(biāo)系下的平面二維水沙數(shù)學(xué)模型，結(jié)合Osman提出的粘性河岸沖刷模型，用于模擬彎道縱向及橫向的河床變形.通過對(duì)一水槽彎道模型和一概化彎道的模擬，本文得到以下結(jié)論：

(1)采用曲線坐標(biāo)下的水流控制方程，能有效的解決不規(guī)則岸邊界問題；用MADI法求解方程，不僅能減少計(jì)算時(shí)的存儲(chǔ)量，而且減少計(jì)算量.通過對(duì)水槽模型試驗(yàn)資料的對(duì)比，本文提出的水流模型的計(jì)算結(jié)果具有一定的精度和穩(wěn)定性.

(2)采用非均勻懸移質(zhì)泥沙輸移方程，引入張紅武提出的非平衡輸沙模式，克服了以往對(duì)方程中的泥沙恢復(fù)飽和系數(shù)任意調(diào)整的局面.

(3)平面二維水沙數(shù)學(xué)模型與河岸沖刷模型相結(jié)合，能更好的模擬河道的河床變形過程.通過對(duì)一概化彎道的模擬，表明彎道內(nèi)主流線變化不僅與進(jìn)出口水沙條件有關(guān)，而且還與斷面形態(tài)有關(guān).而斷面形態(tài)的變化與河岸抗沖性能有關(guān).對(duì)于彎道在持續(xù)清水沖刷下，當(dāng)河岸可沖刷時(shí)的主槽比降變化比河岸不可沖刷時(shí)要慢.