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

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

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

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

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

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

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

1.1 特征值

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

水流運動方程

泥沙連續(xù)方程

河床變形方程

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

這些基本方程由于簡化程度不同可組合成不同形式的懸移質(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)組成.

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

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

1.2 特征向量

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

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

1.3 特征關(guān)系

特征關(guān)系可以揭示特征線上各變量之間的相互約束關(guān)系.同時，通過特征關(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)生的擾動對水流運動的特征值及特征關(guān)系會產(chǎn)生影響[12].但從式(15)～(18)所示懸移質(zhì)模型的特征關(guān)系來看，上述4種模型在4根特征線上均能綜合反映水流運動、泥沙輸運及河床變形間相互作用的情況.從式(18)所示的特征關(guān)系不難看出，CM_1河床沖淤變形并未受到水流運動的直接影響.

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

2.1 模型耦合分析

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

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

事實上，當(dāng)河床表面存在沙波運動時[16]，式(19)中ubΔz所表示的通量是由沙波運動推動的，而沙波運動與水流的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é)模型的特征分析

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

其中

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

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

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

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

圖1 不同沖、淤強度下特征值與流速的關(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個特征值綜合了水流運動、泥沙輸運及河床沖淤的影響.式(21)和(24)表明，水流運動明顯受到河床變形及泥沙輸運的影響;式(22)表明，水流運動與泥沙輸運同時會影響河床變形的演變.應(yīng)當(dāng)指出的是，若動量方程(5)中不包括左端第5項(泥沙輸運項)和第6項(河床變形項)，式(23)將退化為u，即不平衡輸沙將不受水流運動與河床變形的影響.

3 結(jié) 論

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

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

c.本文基于奇異攝動理論，通過漸進(jìn)展開方法求得的適當(dāng)簡化后的模型所構(gòu)成的雙曲系統(tǒng)的4個特征值，分別表示了4個以不同波速傳播的波的傳播特點.與多數(shù)懸移質(zhì)不平衡輸沙模型的特征分析結(jié)果相比，耦合模型CM_2計算所得的特征值能夠充分反映和定量描述水流運動、泥沙輸運及河床變形的相互作用機制.河床沖淤強度越大，水流所受影響越大.λ′1與λ′2均有一個過渡區(qū)u=hF-1±O(ε1/2).河床變形主要受λ′1，λ′2所表示的波運動的影響，λ′3和λ′4所表示的波運動幾乎不會對河床變形產(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]楊國錄.河流數(shù)學(xué)模型[M].北京:海洋出版社，1993.

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

[4]江恩惠，趙連軍，韋直林.黃河下游洪峰增值機理與驗證[J].水利學(xué)報，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]萬兆惠，宋天成.“揭河底”沖刷現(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]譚維炎.計算淺水動力學(xué)[M].北京:清華大學(xué)出版社，1998.

[13]許協(xié)慶，朱鵬程.河床變形問題的特征線解[J].水利學(xué)報，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]王士強.沙波運動與床沙交換調(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]錢寧，萬兆惠.泥沙運動力學(xué)[M].北京:科學(xué)出版社，2003.

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