丁，戴文鴻，鐘德鈺，唐立模，陳洪兵

(1.河海大學(xué)水利水電學(xué)院，江蘇南京 210098;2.清華大學(xué)水沙科學(xué)與水利水電工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100084; 3.上海市水利工程設(shè)計(jì)研究院，上海 200051)

目前，水沙數(shù)學(xué)模型已廣泛應(yīng)用于天然沖積河流的水流運(yùn)動(dòng)、泥沙輸運(yùn)及河床變形模擬計(jì)算.現(xiàn)有常見(jiàn)的水沙數(shù)學(xué)模型多數(shù)采用了2種處理方法:(a)基本控制方程簡(jiǎn)化，例如在水流連續(xù)方程中省略了河床可動(dòng)性影響，運(yùn)動(dòng)方程中忽略泥沙輸運(yùn)及河床變形項(xiàng)等[1-2];(b)求解方法解耦，即在一個(gè)時(shí)間步長(zhǎng)內(nèi)將水流計(jì)算和泥沙計(jì)算分開(kāi)，先求解水流方程得到相應(yīng)流動(dòng)要素后再求解輸沙方程及河床變形方程[2-3].對(duì)于少沙河流而言，其水流含沙量較低，河床沖淤?gòu)?qiáng)度較弱且時(shí)間尺度較長(zhǎng)，這種處理方法具有一定的適用性.但對(duì)于諸如黃河下游這樣的多沙河流而言，水流運(yùn)動(dòng)、泥沙輸運(yùn)及河床變形之間存在強(qiáng)耦合、非線性的復(fù)雜過(guò)程，同時(shí)伴隨一些異常現(xiàn)象的發(fā)生[4-5]，這種處理方法的有效性值得進(jìn)一步研究.研究[6-8]表明，這種處理方式破壞了模型的耦合性，在進(jìn)行多沙河流模擬時(shí)存在較大的局限性.

針對(duì)以上問(wèn)題，國(guó)內(nèi)外開(kāi)展了相關(guān)研究[9-10].本文在前人研究的基礎(chǔ)上，著重從數(shù)學(xué)模型的理論基礎(chǔ)入手進(jìn)行探討.從數(shù)學(xué)性質(zhì)上看，一維水沙數(shù)學(xué)模型構(gòu)成一個(gè)雙曲占優(yōu)的系統(tǒng)，其重要特點(diǎn)是具有特征效應(yīng)，即變量擾動(dòng)沿特征線并以特征值這樣的速度進(jìn)行傳播[11].因此，多沙河流數(shù)學(xué)模型耦合與否的一個(gè)重要標(biāo)志是，水流運(yùn)動(dòng)、泥沙輸運(yùn)及河床變形之間的相互作用能否綜合反映到雙曲系統(tǒng)的特征關(guān)系上.而對(duì)于水沙數(shù)值模擬而言，能否充分反映水流運(yùn)動(dòng)與河床變形間的相互作用及耦合機(jī)制成為有效模擬的關(guān)鍵所在[3].

本文根據(jù)雙曲系統(tǒng)的特征理論及奇異攝動(dòng)理論系統(tǒng)推求了目前數(shù)學(xué)模型的特征值、特征向量及特征關(guān)系，分析了水流運(yùn)動(dòng)、泥沙輸運(yùn)及河床變形的相互作用機(jī)制，并進(jìn)一步研究了一維耦合數(shù)學(xué)模型的特征值.

1 懸移質(zhì)模型的特征值、特征向量和特征關(guān)系

目前對(duì)于推移質(zhì)模型及全沙模型特征分析的結(jié)論揭示了水流運(yùn)動(dòng)與河床變形之間存在的復(fù)雜相互作用機(jī)制[12].由于多數(shù)多沙河流中的泥沙輸運(yùn)主要以懸移質(zhì)運(yùn)動(dòng)為主，因此，懸移質(zhì)不平衡輸沙模型對(duì)于多沙河流中水沙運(yùn)動(dòng)現(xiàn)象的描述和模擬更具有針對(duì)性，且懸移質(zhì)模型目前被廣泛應(yīng)用于多沙河流洪水演進(jìn)與河床沖淤變形計(jì)算.對(duì)于多沙河流水流運(yùn)動(dòng)與河床變形的相互作用機(jī)制的研究而言，懸移質(zhì)模型的特征分析則更具有理論價(jià)值.

1.1 特征值

現(xiàn)有一維水沙數(shù)學(xué)模型常見(jiàn)的基本控制方程為水流連續(xù)方程

水流運(yùn)動(dòng)方程

泥沙連續(xù)方程

河床變形方程

式中:t——時(shí)間;x——縱向坐標(biāo);h——水深;u——流速;zb——河床高程;i0——床面坡降;if——摩阻能坡;g——重力加速度;ρs——泥沙密度;ρw——清水密度;ρm——渾水密度;ρb——床沙飽和密度;p——床沙孔隙率;sv——含沙量;s*v——挾沙力;α——恢復(fù)飽和系數(shù);ω——泥沙沉速.

這些基本方程由于簡(jiǎn)化程度不同可組合成不同形式的懸移質(zhì)輸沙模型.本文用DM(包括DM_1，DM_2和DM_3)表示非耦合模型，CM(包括CM_1和CM_2)表示耦合模型.其中:DM_1[1]由式(1)，(3)，(6)，(8)組成;DM_2[3]由式(1)，(4)，(6)，(8)組成;DM_3[2]由式(1)，(4)，(6)，(9)組成;CM_1[9-10]由式(2)，(5)，(7)，(8)組成.

從計(jì)算結(jié)果可知:DM_1，DM_2，DM_3和CM_1這4種模型的λ4(河床變形特征值)均為0，表示河床沖淤在縱向上不受水流運(yùn)動(dòng)與泥沙輸運(yùn)的影響;DM_1和CM_1的λ1，2(水流特征值)均為，與圣維南方程組特征值相同，λ3(輸沙特征值)均為 u，未受河床沖淤的影響;DM_2和DM_3的 λ1，2均為 u±，雖然不同于DM_1，但均未受水流運(yùn)動(dòng)及河床沖淤的影響，λ3均為u.因此，從雙曲系統(tǒng)特征理論角度看，這4種模型均沒(méi)有達(dá)到特征耦合.

許協(xié)慶等[13]建立了飽和輸沙模型(由式(2)，(5)和(10)組成的耦合模型)，并采用該模型分析了河床變形和水流運(yùn)動(dòng)的影響，其研究結(jié)果揭示了河床變形、泥沙輸運(yùn)與水流運(yùn)動(dòng)的相互關(guān)系.此外，推移質(zhì)模型中河床變形對(duì)水流的影響在其特征關(guān)系上得到了體現(xiàn)[14].

1.2 特征向量

根據(jù)上述特征值可進(jìn)一步導(dǎo)出DM_1，DM_2，DM_3和CM_1分別對(duì)應(yīng)于特征值λ1，λ2，λ3和λ4的右特征列向量矩陣R1，R2，R3和R4，如式(11)～(14)所示.

式(11)～(14)右邊第1列和第4列分別表示水流中的擾動(dòng)向上游和下游傳播的軌跡.與特征值的規(guī)律類(lèi)似，除DM_3外，其他模型的水流運(yùn)動(dòng)未受泥沙輸運(yùn)和河床變形的影響.而DM_3河床沖淤對(duì)水流運(yùn)動(dòng)的影響也未得到體現(xiàn).

1.3 特征關(guān)系

特征關(guān)系可以揭示特征線上各變量之間的相互約束關(guān)系.同時(shí)，通過(guò)特征關(guān)系，可將原有非線性偏微分方程轉(zhuǎn)化成常微分方程，也可使得各變量的相互約束關(guān)系得到更直觀和定量的反映.

DM_1的特征關(guān)系為

DM_2的特征關(guān)系為

DM_3的特征關(guān)系為

CM_1的特征關(guān)系為

從推移質(zhì)模型和全沙模型的特征關(guān)系可以看出，河床沖淤變形產(chǎn)生的擾動(dòng)對(duì)水流運(yùn)動(dòng)的特征值及特征關(guān)系會(huì)產(chǎn)生影響[12].但從式(15)～(18)所示懸移質(zhì)模型的特征關(guān)系來(lái)看，上述4種模型在4根特征線上均能綜合反映水流運(yùn)動(dòng)、泥沙輸運(yùn)及河床變形間相互作用的情況.從式(18)所示的特征關(guān)系不難看出，CM_1河床沖淤變形并未受到水流運(yùn)動(dòng)的直接影響.

2 一維耦合懸移質(zhì)模型的特征值

2.1 模型耦合分析

現(xiàn)有數(shù)學(xué)模型采用的河床變形方程是根據(jù)床面運(yùn)動(dòng)學(xué)邊界條件和某些假定(如平衡輸沙假定)建立的一種床面局部平衡關(guān)系式[2-3]，未能反映河床縱向變化的影響.這是導(dǎo)致現(xiàn)有模型無(wú)法耦合的重要原因.為此，丁等[15]導(dǎo)出了守恒形式的河床沖淤層連續(xù)方程

式中:Δ z——河床沖淤厚度;ub——河床縱向運(yùn)動(dòng)速度;D——沉降通量;E——沖刷通量;ρ′b——床沙干密度.

事實(shí)上，當(dāng)河床表面存在沙波運(yùn)動(dòng)時(shí)[16]，式(19)中ubΔz所表示的通量是由沙波運(yùn)動(dòng)推動(dòng)的，而沙波運(yùn)動(dòng)與水流的Froude數(shù)密切相關(guān)[17]，因此ubΔz可表示成

式中k，m，n是與水流流態(tài)、含沙量等有關(guān)的量.

將式(20)代入式(19)，聯(lián)立式(2)，(5)和(7)，可組成考慮河床縱向沖淤的懸移質(zhì)耦合輸沙模型(記為CM_2).

2.2 一維耦合數(shù)學(xué)模型的特征分析

為便于討論，對(duì)CM_2中的各變量進(jìn)行量綱為1的處理.根據(jù)奇異攝動(dòng)理論[18]，通過(guò)漸進(jìn)展開(kāi)方法求得其4個(gè)特征值，分別為

其中

式中:H，U——恒定流水深及流速;F——恒定流弗勞德數(shù);ε——河床變形對(duì)水流運(yùn)動(dòng)的擾動(dòng)量，河床淤積時(shí)ε＞0(σ=1)，河床沖刷時(shí)ε＜0(σ=-1).

式(21)～(24)中的λ′1表示水流變化產(chǎn)生的負(fù)向波動(dòng)，λ′2表示河床變形波動(dòng)，λ′3表示不平衡輸沙的波動(dòng)，λ′4表示水流變化產(chǎn)生的正向波動(dòng).它們有2個(gè)顯著特點(diǎn):(a)各變量相互影響;(b)λ′1，λ′2有一個(gè)過(guò)渡區(qū)域，即.這表明，水流流態(tài)由緩流向急流變化時(shí)，水流運(yùn)動(dòng)和河床變形所受擾動(dòng)均受到不同程度的影響.

現(xiàn)以λ′1為例說(shuō)明水流運(yùn)動(dòng)與河床變形的關(guān)系.λ′1的結(jié)構(gòu)表明，水流所受擾動(dòng)在負(fù)向上的傳播機(jī)制是不同的.當(dāng)水流處于緩流或急流流態(tài)時(shí)，水流本身的擾動(dòng)占主導(dǎo)地位;當(dāng)水流處于過(guò)渡態(tài)時(shí)，河床變形的擾動(dòng)會(huì)對(duì)水流運(yùn)動(dòng)產(chǎn)生較大影響.

河床不同淤積和沖刷強(qiáng)度對(duì)水流的影響如圖1所示.由圖1可以看出:當(dāng)沖刷強(qiáng)度增大時(shí)，水流負(fù)向、河床變形正向特征值減小，不平衡輸沙特征值變大，水流正向特征值略微變小;而當(dāng)淤積強(qiáng)度增大時(shí)，水流負(fù)向、河床變形正向特征值顯著增大，不平衡輸沙特征值變小，水流正向特征值略微變大.

圖1 不同沖、淤?gòu)?qiáng)度下特征值與流速的關(guān)系Fig.1 Variations of eigenvalues with velocity under different deposition strengths and scouring strengths

與DM_1，DM_3及CM_1的特征分析結(jié)果不同的是，本文所建立的耦合模型所構(gòu)成的雙曲系統(tǒng)的4個(gè)特征值綜合了水流運(yùn)動(dòng)、泥沙輸運(yùn)及河床沖淤的影響.式(21)和(24)表明，水流運(yùn)動(dòng)明顯受到河床變形及泥沙輸運(yùn)的影響;式(22)表明，水流運(yùn)動(dòng)與泥沙輸運(yùn)同時(shí)會(huì)影響河床變形的演變.應(yīng)當(dāng)指出的是，若動(dòng)量方程(5)中不包括左端第5項(xiàng)(泥沙輸運(yùn)項(xiàng))和第6項(xiàng)(河床變形項(xiàng))，式(23)將退化為u，即不平衡輸沙將不受水流運(yùn)動(dòng)與河床變形的影響.

3 結(jié) 論

a.多數(shù)懸移質(zhì)輸沙數(shù)學(xué)模型在特征上是解耦的.

b.懸移質(zhì)輸沙數(shù)學(xué)模型所采用的河床變形方程是基于局部平衡關(guān)系得到的，并不能揭示河床沖淤縱向變化所引起的質(zhì)量和動(dòng)量變化機(jī)制.

c.本文基于奇異攝動(dòng)理論，通過(guò)漸進(jìn)展開(kāi)方法求得的適當(dāng)簡(jiǎn)化后的模型所構(gòu)成的雙曲系統(tǒng)的4個(gè)特征值，分別表示了4個(gè)以不同波速傳播的波的傳播特點(diǎn).與多數(shù)懸移質(zhì)不平衡輸沙模型的特征分析結(jié)果相比，耦合模型CM_2計(jì)算所得的特征值能夠充分反映和定量描述水流運(yùn)動(dòng)、泥沙輸運(yùn)及河床變形的相互作用機(jī)制.河床沖淤?gòu)?qiáng)度越大，水流所受影響越大.λ′1與λ′2均有一個(gè)過(guò)渡區(qū)u=hF-1±O(ε1/2).河床變形主要受λ′1，λ′2所表示的波運(yùn)動(dòng)的影響，λ′3和λ′4所表示的波運(yùn)動(dòng)幾乎不會(huì)對(duì)河床變形產(chǎn)生影響.

[1]CUNGE JA，HOLLY F M，JOINT A A，et al.Practical aspects of computational river hydraulics[M].Boston:Pitman Advanced Publish. Program，1980.

[2]楊國(guó)錄.河流數(shù)學(xué)模型[M].北京:海洋出版社，1993.

[3]謝鑒衡.河流模擬[M].北京:水利電力出版社，1990.

[4]江恩惠，趙連軍，韋直林.黃河下游洪峰增值機(jī)理與驗(yàn)證[J].水利學(xué)報(bào)，2006(12):1454-1459.(JIANG En-hui，ZHAO Lian-jun，WEI Zhi-lin.Mechanism of flood peak increase along the lower Yellow River and its verification[J].Journal of Hydraulic Engineering，2006(12):1454-1459.(in Chinese))

[5]萬(wàn)兆惠，宋天成.“揭河底”沖刷現(xiàn)象的分析[J].泥沙研究，1991(3):20-27.(WAN Zhao-hui，SONG Tian-cheng.On“ripping up the bottom”phenomenon[J].Journal of Sediment Research，1991(3):20-27.(in Chinese))

[6]BUSNELLI M M，STELLING S G，MICHELC L，et al.Numerical morphological modeling of open-check dams[J].Journal of Hydraulic Engineering，2001，127(2):105-114.

[7]SIEBEN J.A theoretical analysis of discontinuous flow with mobile bed[J].Journal of Hydraulic Research，1999，37(2):199-212.

[8]KASSEM A A，CHAUDHRY M H.Comparison of coupled and semicoupled numerical models for alluvial channels[J].Journal of Hydraulic Engineering，1998，124(8):794-802.

[9]CAO Z，EGASHIRA S.Coupled mathematical modelling of alluvial rivers[J].Journal of Hydroscience and Hydraulic Engineering，1999， 17(2):71-85.

[10]WU Wei-ming，WANG S S Y.One dimensional modeling of dam-break flow over movable beds[J].Journal of Hydraulic Engineering，2007，13(1):48-58.

[11]譚維炎.計(jì)算淺水動(dòng)力學(xué)[M].北京:清華大學(xué)出版社，1998.

[13]許協(xié)慶，朱鵬程.河床變形問(wèn)題的特征線解[J].水利學(xué)報(bào)，1964(5):1-19.(XU Xie-qing，ZHU Peng-cheng.Study of riverbed degradation and aggradation by the method of characteristics[J].Journal of Hydraulic Engineering，1964(5):1-19.(in Chinese))

[14]VRIES M.On morphological predictions[C]//A Compilation of the Lectures on River Sedimentation.Beijing:Sedimentation Committee Chinese Society of Hydraulic Engineering，1981:93-158.

[16]王士強(qiáng).沙波運(yùn)動(dòng)與床沙交換調(diào)整[J].泥沙研究，1992(4):14-23.(WANG Shi-qiang.The movement of sandwaves and exchanges of bed material[J].Journal of Sediment Research，1992(4):14-23.(in Chinese)).

[17]錢(qián)寧，萬(wàn)兆惠.泥沙運(yùn)動(dòng)力學(xué)[M].北京:科學(xué)出版社，2003.

[18]HASAN N A.Introduction to perturbation techniques[M].New York:John Wiley&Sons，1981.