張 雪，王新宏，董 梁，郎永彪，倪宇航

(1.西安理工大學(xué) 省部共建西北旱區(qū)生態(tài)水利國家重點實驗室，陜西 西安 710048；2.青海省水利水電勘測設(shè)計研究院，青海 西寧 810001)

1 研究背景

跌坎沖刷現(xiàn)象廣泛存在于水庫溯源沖刷[1]、堤壩潰決[2-3]以及河道過度采砂[4]導(dǎo)致的河床沖刷等物理過程中。為了解跌坎沖刷過程及其機(jī)理，國內(nèi)外學(xué)者對不同尺度跌坎進(jìn)行試驗研究[5]，其中Bennett等[6]進(jìn)行的跌坎高度為5～50 mm的小尺度跌坎沖刷水槽試驗以及Robinson[7]進(jìn)行的0.96～1.55 m的大尺度跌坎沖刷水槽試驗比較有代表性。關(guān)于溯源沖刷問題，我國學(xué)者早在20世紀(jì)60年代就已開始研究[8-9]，近年來也有很多學(xué)者針對實際水流中出現(xiàn)的跌坎沖刷現(xiàn)象進(jìn)行了研究[4,10-11]。由于跌坎沖刷過程中伴隨著復(fù)雜的跨臨界流以及劇烈的河床變形過程，給數(shù)值求解帶來了極大的困難。目前針對跨臨界流模擬主要有兩種方法——激波捕捉法和激波擬合法[12-13]。前者可以自動捕捉間斷，但要求格式具有高性能且滿足三守恒(守恒因變量、守恒型微分方程、守恒型格式)，否則間斷解精度不高或不正確[14]；后者在間斷點引入RH條件及單側(cè)格式，可準(zhǔn)確判斷間斷位置及傳播速度，適合于形狀陡峭的數(shù)值間斷，算法復(fù)雜[15]。目前求解淺水方程的格式有Godunov 格式[16-17]、交替方向隱式格式(ADI)[18]、通量差分裂格式[19-21]和Preissmann格式等，其中Preissmann格式結(jié)構(gòu)最簡單且計算精度高、對步長沒有限制，因而應(yīng)用最為廣泛[22]。但基于Preissmann格式離散的圣維南方程組在不修正的情況下無法適用于跨臨界流等間斷流態(tài)[23-24]。為解決上述問題，茅澤育等[25]針對Preissmann格式處理混合流時的不適進(jìn)行了定性分析，認(rèn)為邊界條件數(shù)目與方程總數(shù)不同導(dǎo)致此格式無法模擬跨臨界流等混合流動。吳曉玲等[26]將Preissmann格式與特征方向分解方法融合，提出了復(fù)雜流態(tài)自適應(yīng)模型，并通過理想算例及實際算例驗證了所建模型可模擬急、緩流交替現(xiàn)象。Sart等[27]根據(jù)特征方向進(jìn)行局部分解的方法，對跨臨界流區(qū)域進(jìn)行流態(tài)劃分，但如果邊界河段流態(tài)復(fù)雜則無法使用該方法。

本文試圖建立適用于跌坎沖刷的一維非恒定推移質(zhì)輸移數(shù)學(xué)模型，模擬動床情況下的跨臨界流問題。對水流方程的求解，考慮到跌坎沖刷現(xiàn)象中存在跨臨界流的問題，依據(jù)Fr大小，將河道劃分為不同流態(tài)的河段，基于Preissmann格式分別離散求解。對河床變形方程的求解，采用向后差分格式離散求解，并利用水槽試驗數(shù)據(jù)對模型的可靠性進(jìn)行驗證。本研究對于水庫溯源沖刷、堤壩潰決等物理過程中出現(xiàn)的跌坎沖刷過程的預(yù)測具有重要的參考意義。

2 數(shù)學(xué)模型構(gòu)建

跌坎沖刷通常發(fā)生在河床高程沿程忽然下降的河段，出現(xiàn)跌坎沖刷的河段可以概化為由3段不同底坡的河段組成，中間為陡坡河段，其上游和下游均為緩坡河段，如圖1(a)所示。通常情況下，上游緩坡河段水流流態(tài)為緩流，中間陡坡河段水流為急流，下游河段水流由急流轉(zhuǎn)變成緩流，有水躍現(xiàn)象發(fā)生。由于陡坡段水流流速突然增大，河床沖刷變形劇烈。在水流沖刷作用下上游變坡點(中間陡坡段與上游緩坡段床面交點)的位置不斷向上游移動，臨界水深位置和沖刷逐步向上游發(fā)展，下游河床產(chǎn)生淤積，最終達(dá)到?jīng)_淤平衡，如圖1(b)。本文假定導(dǎo)致跌坎沖刷、河床沖淤變形的主要原因是推移質(zhì)運動，不考慮懸移質(zhì)運動造成的河床沖刷。在此前提下構(gòu)建合適的推移質(zhì)輸移數(shù)學(xué)模型，對跌坎沖刷過程進(jìn)行模擬，給出不同時刻水面線以及河床高程變化。

圖1 河道跌坎水流流態(tài)及沖刷示意圖

2.1 控制方程

本文建立一維非恒定推移質(zhì)輸移數(shù)學(xué)模型，采用的水流控制方程為水流連續(xù)方程和水流運動方程(統(tǒng)稱為圣維南方程組)，泥沙控制方程為河床變形方程，輸沙率方程采用Meyer-Peter Müller(MPM)推移質(zhì)輸沙率公式，基本控制方程如下。

水流連續(xù)方程：

(1)

水流運動方程：

(2)

河床變形方程：

(3)

輸沙率公式：

(4)

2.2 數(shù)值方法

2.2.1 圣維南方程組離散 采用Preissmann四點隱式差分格式對水流控制方程式(1)、(2)進(jìn)行離散，并對離散后的非線性代數(shù)方程組進(jìn)行進(jìn)一步線性化處理后可以得到線性代數(shù)方程組(5)、(6)。

(5)

(6)

式中：Δzj、Δzj+1分別為第n+1時間層j、j+1斷面的水位增量；ΔQj、ΔQj+1分別為第n+1時間層j、j+1斷面的流量增量；A′1、B′1、…、F′2為線性代數(shù)方程組的系數(shù)，其值只與第n時間層結(jié)點值有關(guān)，通過代入第n時間層各斷面流量及水位即可計算，各系數(shù)表達(dá)式見表1。

表1 線性代數(shù)方程組(5)、(6)中各系數(shù)表達(dá)式

方程組(5)、(6)是針對1個河段的兩個計算斷面j、j+1上寫出的方程，假設(shè)將河段劃分為j個斷面(j-1個河段)，每個斷面含有兩個未知數(shù)，各斷面編號分別為0，1，2，…，j，每個河段可列兩個方程，總共2j-2個方程，補(bǔ)充上、下游兩個邊界條件使方程組封閉，即可聯(lián)立求解得出各斷面不同時刻的流量和水位。

2.2.2 河床變形方程離散 采用向后差分格式對河床變形方程進(jìn)行離散，離散后的方程式(3)可化為方程式(7)、(8)：

(7)

(8)

河床高程Y0為已知，將水流連續(xù)方程求解得出的流量及水位代入方程式(7)，采用迭代法進(jìn)行計算。

2.3 跨臨界流處理方法

2.3.1 數(shù)值振蕩原因分析 由于在跌坎附近水流變?yōu)榧绷?，計算區(qū)域出現(xiàn)急、緩流交替，水流模擬進(jìn)行到6.68 s時數(shù)值發(fā)散，出現(xiàn)數(shù)值振蕩現(xiàn)象，如圖2所示。

圖2 跌坎水流模擬數(shù)值振蕩細(xì)節(jié)圖(模擬至6.68 s時)

采用水槽試驗數(shù)據(jù)進(jìn)行模擬計算驗證[28]，計算條件如下：計算長度為6 m，時間步長為0.005 s，空間步長為0.01 m，總共601個斷面，分別為1#、2#、…、601#斷面；跌坎坡腳位于3 m處，跌坎高度為0.1 m，坡度為30°；河床糙率為0.017，加權(quán)系數(shù)θ取0.65；初始水位為0.13 m，入口流量為5 m3/h，出口控制水深為0.05 m。模擬結(jié)果顯示，260#斷面之前水位維持在0.12 m左右，未發(fā)生明顯波動。263#斷面水位開始出現(xiàn)波動，263#～270#斷面水位上下波動范圍在0.02 m以內(nèi)，270#～282#斷面水位波動很大，水位誤差超過0.35 m。282#～286#斷面波動范圍逐漸縮小，水位降到0.119 m左右。跌坎處280#、282#斷面水位出現(xiàn)負(fù)值，程序中斷無法計算。出現(xiàn)數(shù)值振蕩的原因在于進(jìn)行水流模擬時邊界條件設(shè)置錯誤，由于跌坎處水流出現(xiàn)急流，而給定的邊界條件只適用于緩流，急流時應(yīng)給定兩個上游邊界條件。

根據(jù)茅澤育等[15]對Preissmann格式處理間斷水流的不適定性分析，當(dāng)水流流態(tài)處于單一緩流時，給定一個上邊界條件和一個下邊界條件(通常給定為上游流量過程線和下游水位過程線)，方程數(shù)與未知數(shù)相等，可得唯一解；水流為完全急流時，只在上游給定兩個邊界條件(通常給定為流量過程線和水位過程線)，與緩流時一樣可封閉求解；當(dāng)?shù)沧兤曼c處水流由緩流過渡為急流時，入流點(緩流)給定上邊界條件，方程數(shù)(2j-1)小于未知數(shù)數(shù)目2j；當(dāng)水流由急流過渡到緩流時，入流點(急流)給定上邊界條件，出流點(緩流)給定一個下邊界條件，方程數(shù)(2j+1)大于未知數(shù)數(shù)目2j；上述兩種跨臨界流均不可得到唯一解，因此需要對所建立的水流模型進(jìn)行修正。

2.3.2 圣維南方程組離散求解方法修正 根據(jù)上述討論，Preissmann格式只能直接應(yīng)用于每一個緩流或者急流區(qū)域，對于跨臨界流區(qū)域需要增加附加方程才能求解。首先計算出第n時間層各斷面的弗勞德數(shù)Fr，依據(jù)Fr沿流程變化情況，可將河段按水流流態(tài)劃分為緩流河段、急流河段、水躍河段、水跌河段4種情況，確保各流段為單一流態(tài)的連續(xù)水流，針對每種情況確定不同的方程形式，從而建立求解問題的方程組。

(1)連續(xù)緩流河段或連續(xù)急流河段計算方程。如果Frj<1(j=0，1，2，…，j)，即各斷面弗勞德數(shù)均為緩流，各斷面均可列出(9)、(10)兩個方程組，給定上游邊界條件為流量過程線Q=Q(t)，下游邊界條件為水位過程線z=z(t)，采用經(jīng)典的追趕迭代法即可求解各斷面水位及流量。

=0

(9)

=0

(10)

如果Frj≥1(j=0，1，2，…，j)，即各斷面弗勞德數(shù)均為急流，則需要給定兩個上游邊界條件，即流量過程線Q=Q(t)和水位過程線z=z(t)，通過迭代法進(jìn)行求解。

(11)

(12)

其中：

(13)

(14)

a(1)=v-c

(15)

a(2)=v+c

(16)

式中：a(1)為微波向上游傳播的速度，m/s；a(2)為微波向下游傳播的速度，m/s；j+1/2為河段；j、j+1為斷面；v為斷面平均流速，v=Q/A，m/s；c為微波的相對速度，m/s。

(17)

α(Δzj+2bj+2-ΔQj+2/2cj+2)+(1-α)(Δzj+1bj+1-ΔQj+1/2cj+1)=0

(18)

(19)

(3)水跌河段計算方程。如果Frj<1、Frj+1≥1，則有臨界點(Fr=1)存在，出現(xiàn)水跌現(xiàn)象，但水流仍處于連續(xù)狀態(tài)，在此河段上仍然可以使用水流連續(xù)方程及水流運動方程。此河段上游邊界為緩流，下游邊界為急流，上游給定流量過程線，由于方程數(shù)為(2j-1)小于未知數(shù)數(shù)目2j，需要補(bǔ)充一個方程作為下游邊界條件即可求解。臨界點處a(1)=0,a(2)=2c，補(bǔ)充方程為：

α(Δzj+1bj+1-ΔQj+1/2cj+1)+

(1-α)(Δzjbj-ΔQj/2cj)=0

(20)

(21)

通過流態(tài)將整個河段劃分，不同區(qū)域采用不同的方程求解，避免了Preissmann格式不能適用于跨臨界流的缺點。

2.4 邊界條件及求解方法

在跨臨界流問題中將河道劃分為不同流態(tài)的河段，依據(jù)不同流態(tài)提出相適應(yīng)的邊界條件。河段為緩流時給定上游流量過程線和下游水位控制線，急流時給定兩個上游邊界條件，即水位和流量。水流為水躍和水跌時需加入附加方程進(jìn)行求解。

求解方法采用非耦合解法，首先考慮水動力因素求解水流控制方程，其次將水流方程中計算得到的水位、流量帶入泥沙方程進(jìn)行迭代求解，此方法大大減少了計算難度。

3 模型試驗驗證

3.1 基本資料

本文對跌坎沖刷問題的模擬采用水槽試驗數(shù)據(jù)進(jìn)行驗證[28]。矩形斷面水槽長度為6 m，水面寬度為1 m，跌坎末端位于水槽中點，跌坎坡度為30°，跌坎高度為0.1 m，泥沙顆粒中值粒徑采用0.6 mm(采用均勻沙)，河床糙率取0.017；時間步長取0.005 s，空間步長取0.01 m(即斷面總數(shù)為601)，河床底部糙率取0.017，加權(quán)系數(shù)θ取0.65；入口流量為5 m3/h，上游水深為0.027 m，跌坎下游水深控制變化過程見表2，其中t為時間，Z為水深。河床底部高程由公式(22)給定，x為距進(jìn)口里程方向的距離。

表2 模型試驗跌坎下游水深控制變化過程

(22)

3.2 結(jié)果分析

3.2.1 水流流態(tài)分析 根據(jù)上述水槽基本資料及各參數(shù)對本文所建模型進(jìn)行了驗證計算，120 s時各斷面Fr分布情況見圖3。

圖3 水槽試驗120 s時各斷面Fr分布模擬計算結(jié)果

由圖3可以看出，在260#斷面(距離進(jìn)口2.6 m處)之前，各斷面Fr均小于0.3，上游斷面水流處于均緩流狀態(tài)，260#～280#斷面Fr迅速增長接近于1，出現(xiàn)這一劇增的原因在于跌坎的出現(xiàn)。281#斷面Fr為0.978，282#斷面Fr為1.823。281#～282#斷面中間經(jīng)歷緩流到急流的過渡，出現(xiàn)跨臨界流，該位置正處于跌坎上游邊坡點附近。283#～292#斷面Fr均為大于1且不斷增大到4.856，293#斷面Fr開始劇烈下降，294#斷面Fr急劇降到0.5左右，水流流態(tài)變?yōu)榫徚?，該位置處于跌坎下游坡腳處，水流從急流轉(zhuǎn)變?yōu)榫徚鞒霈F(xiàn)水躍現(xiàn)象。之后水流一直處于緩流流態(tài)，最終Fr維持在0.05基本不變。整個沖刷過程中水流流態(tài)從緩流到急流再到緩流，期間出現(xiàn)了水躍、水跌等水力現(xiàn)象，表明經(jīng)過修正后的模型可以很好地模擬復(fù)雜流態(tài)的轉(zhuǎn)變。

3.2.2 水面線及河床底部高程分析 圖4為t=0、t=5 s、t=10 s、t=50 s以及t=120 s時的跌坎沖刷模擬計算結(jié)果，其中包括水面線計算值與實測值、河床底部高程計算值與實測值。初始時刻水面線及河床底部高程見t=0時刻(圖4(a))，水面高程保持在0.13 m，河床底部高程由公式(22)給定。跌坎沖刷進(jìn)行到5 s時(圖4(b))，水面線稍有下降，可以看到水面線在上游變坡點處(284#斷面)出現(xiàn)明顯降落，水位開始出現(xiàn)波動，底部河床還未見明顯沖刷。沖刷進(jìn)行到10 s時(圖4(c))，上游變坡點處水位沿河床邊界大幅下降，河床底部也開始出現(xiàn)沖刷，下游水位降到0.08 m左右。沖刷到50 s時(圖4(d))，坡腳處水位已降到0.03 m左右，水位沿程下降，河床大幅沖刷，且沖刷范圍不斷向上游發(fā)展，下游河床發(fā)生淤積，跌坎坡度變緩。模擬進(jìn)行到120 s時(圖4(e))，上游水位降到0.11 m左右，下游水位降到0.01 m，上游沖刷及下游淤積速度明顯減慢，基本達(dá)到?jīng)_刷平衡。本文在水動力模型中考慮河床沖淤變形，修正后的模型在模擬動床跨臨界流也可以適應(yīng)復(fù)雜流態(tài)的轉(zhuǎn)變，所建模型對于跌坎沖刷的模擬符合其實際物理現(xiàn)象。

圖5為不同時刻跌坎前后局部水面線高程及河床底部高程。由圖5(a)可看出，跌坎附近水面線隨著沖刷時間逐漸下降，跌坎上、下游水位高差不斷增大，且有明顯的間斷點。t=80 s時水面線在下游坡腳處(3.0～3.5 m斷面)低于t=120 s時的水面高程，其原因是此時下游坡腳處發(fā)生沖刷，此處河床產(chǎn)生淤積所致。由圖5(b)可看出，沖刷初期河床底部高程變化劇烈，隨著時間的延長沖淤速度逐漸變慢；在沖刷作用下跌坎坡度逐漸變緩，沖刷向上游發(fā)展，下游河床發(fā)生淤積，最終達(dá)到基本穩(wěn)定，表明了跌坎沖刷的發(fā)展過程及特點。

圖5 不同時刻跌坎前后局部水面線高程及河床底部高程

3.2.3 計算精度分析 將t=120 s時跌坎前后局部水深及河床沖淤厚度的計算值與實測值進(jìn)行對比，見表3。由表3可見，試驗水槽內(nèi)水深最深處只有1.5 cm，各斷面水深計算值的絕對誤差為0.06～0.12 cm，相對誤差為-6.7%～40%。河床沖淤厚度絕對誤差最大為-0.90 cm，相對誤差最大為24.3%。由于水深很小，試驗本身測量難度大，實測數(shù)據(jù)也有一定的誤差，導(dǎo)致計算值與實測值的誤差偏大，但是模擬趨勢符合實際情況，從圖4中不同時刻跌坎沖刷模擬結(jié)果可以看到，所建模型對于水槽試驗中跌坎沖刷過程的模擬比較準(zhǔn)確。

圖4 不同時刻跌坎沖刷模擬計算與試驗實測結(jié)果

表3 跌坎前后局部水深及河床沖淤厚度計算值與試驗實測值對比(t=120 s)

4 結(jié) 論

跨臨界流現(xiàn)象經(jīng)常出現(xiàn)在潰壩水流及水庫溯源沖刷等實際工程問題中，其研究價值也得到了國內(nèi)外學(xué)者的關(guān)注。本文采用改進(jìn)的Preissmann格式建立了一維水沙動力學(xué)模型，模擬了跌坎沖刷處水沙運動情況，得到以下結(jié)論：

(1)建立了跨臨界流跌坎沖刷一維非恒定推移質(zhì)輸移數(shù)學(xué)模型。該模型可以計算跌坎沖刷河段跨臨界流水面線、河床縱剖面等隨時間的變化過程。

(2)采用經(jīng)典Preissmann格式模擬跨臨界流會出現(xiàn)數(shù)值振蕩現(xiàn)象，修正后的模型通過弗勞德數(shù)Fr識別間斷點位置，據(jù)此將水流劃分為急流、緩流、水躍以及水跌4種流態(tài)，針對4種流態(tài)河段采用不同的離散方程求解，通過河段邊界水流流態(tài)提出合適的邊界條件，迭代進(jìn)行求解。通過水槽試驗數(shù)據(jù)驗證，所建模型可以很好地模擬復(fù)雜流態(tài)的轉(zhuǎn)變，證明了經(jīng)典Preissmann格式通過修正模擬跨臨界流的可行性。

(3)利用跌坎沖刷水槽試驗數(shù)據(jù)，對所建模型進(jìn)行了驗證計算。計算結(jié)果表明，該模型可以反映跌坎沖刷過程中跨臨界流的各種復(fù)雜水流現(xiàn)象，以及陡坡河段由于泥沙沖淤造成的底坡逐漸變緩、長度逐漸增加等河床變形現(xiàn)象。對比水槽試驗實測值與模型計算值可知，各斷面計算水位和河床高程與實測值較為符合，說明所建模型具有一定的計算精度。