基于有限差分法的跌坎沖刷一維數(shù)值模擬

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

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

1 研究背景

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

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

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

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

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

2.1 控制方程

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

水流連續(xù)方程：

(1)

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

(2)

河床變形方程：

(3)

輸沙率公式：

(4)

2.2 數(shù)值方法

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

(5)

(6)

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

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

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

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

(7)

(8)

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

2.3 跨臨界流處理方法

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

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

采用水槽試驗(yàn)數(shù)據(jù)進(jìn)行模擬計(jì)算驗(yàn)證[28]，計(jì)算條件如下：計(jì)算長(zhǎng)度為6 m，時(shí)間步長(zhǎng)為0.005 s，空間步長(zhǎng)為0.01 m，總共601個(gè)斷面，分別為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ā)生明顯波動(dòng)。263#斷面水位開(kāi)始出現(xiàn)波動(dòng)，263#～270#斷面水位上下波動(dòng)范圍在0.02 m以內(nèi)，270#～282#斷面水位波動(dòng)很大，水位誤差超過(guò)0.35 m。282#～286#斷面波動(dòng)范圍逐漸縮小，水位降到0.119 m左右。跌坎處280#、282#斷面水位出現(xiàn)負(fù)值，程序中斷無(wú)法計(jì)算。出現(xiàn)數(shù)值振蕩的原因在于進(jìn)行水流模擬時(shí)邊界條件設(shè)置錯(cuò)誤，由于跌坎處水流出現(xiàn)急流，而給定的邊界條件只適用于緩流，急流時(shí)應(yīng)給定兩個(gè)上游邊界條件。

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

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

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

=0

(9)

=0

(10)

如果Frj≥1(j=0，1，2，…，j)，即各斷面弗勞德數(shù)均為急流，則需要給定兩個(gè)上游邊界條件，即流量過(guò)程線Q=Q(t)和水位過(guò)程線z=z(t)，通過(guò)迭代法進(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為微波的相對(duì)速度，m/s。

(17)

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

(18)

(19)

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

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

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

(20)

(21)

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

2.4 邊界條件及求解方法

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

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

3 模型試驗(yàn)驗(yàn)證

3.1 基本資料

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

表2 模型試驗(yàn)跌坎下游水深控制變化過(guò)程

(22)

3.2 結(jié)果分析

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

圖3 水槽試驗(yàn)120 s時(shí)各斷面Fr分布模擬計(jì)算結(jié)果

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

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

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

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

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

圖4 不同時(shí)刻跌坎沖刷模擬計(jì)算與試驗(yàn)實(shí)測(cè)結(jié)果

表3 跌坎前后局部水深及河床沖淤厚度計(jì)算值與試驗(yàn)實(shí)測(cè)值對(duì)比(t=120 s)

4 結(jié) 論

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

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

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

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