南水北調(diào)中線總干渠沿線與許多河流交叉，其中在河北穿越七里河的交叉建筑物型式為渠穿河倒虹吸工程。倒虹吸設(shè)計(jì)長(zhǎng)度初選為700m，設(shè)計(jì)洪水標(biāo)準(zhǔn)為百年一遇，設(shè)計(jì)洪峰流量2410m³/s。渠穿河倒虹吸工程的修建對(duì)該處河道水流及河床演變產(chǎn)生的影響，倒虹吸工程的位置、尺寸、埋置深度是干渠設(shè)計(jì)所關(guān)心的重大問(wèn)題。利用二維水沙數(shù)學(xué)模型，可以較好地模擬反映渠、河交叉工程附近水流與河床變形狀況，主要研究：在該河修建渠穿河倒虹吸后，交叉工程附近河段流場(chǎng)流速、壅水及河床沖淤變形程度。通過(guò)河流模擬，分析、評(píng)價(jià)交叉工程設(shè)計(jì)方案并提出工程修改建議。

2 交叉段河道特性

　　在交叉工程附近七里河分為兩汊，屬寬淺型沙質(zhì)河床：河道北槽較大，主流靠左岸。中泓處的河床質(zhì)表層為粗沙，d₅₀在1.0mm左右，灘地及兩岸為中沙及壤土。

　　該河為沖積性河流，河床形態(tài)、演變規(guī)律與洪水造床作用和常年水沙特性密切相關(guān)。據(jù)河道查勘與河床地形圖分析：1963年洪水前期引起河床較強(qiáng)的沖刷，洪水后河道又普遍回淤；河段灘槽高差一般為2～3m,河槽寬850～1200m。交叉工程附近河床組成沿流程分選明顯，由粗至細(xì)變化較大。1980年以來(lái)河道受人類活動(dòng)影響嚴(yán)重，由于河道內(nèi)多處挖取粗沙，原主槽回淤的泥沙被大量挖走，有的挖深達(dá)3～４m。多年小水作用及人類挖沙已經(jīng)使現(xiàn)河道形成人為不連續(xù)窄深槽，窄槽寬度一般為100～300m。這造成在現(xiàn)狀河道小洪水易于歸槽且側(cè)侵蝕較為明顯，見(jiàn)圖1。

3.1.1 控制方程

　　對(duì)于散粒沙河床，床面穩(wěn)定控制的輔助方程可表示成局部區(qū)域穩(wěn)定控制條件

式中ξ—水位；H—水深；u,v-x,y向水深平均流速；U—合速度；U_*—摩阻流速，R—水力半徑；e,ν—水流渦粘系數(shù)、運(yùn)動(dòng)粘滯系數(shù)；β—對(duì)流項(xiàng)修正系數(shù)，S—含沙量，z—河床高程，J—水力坡度；D，α—泥沙擴(kuò)散系數(shù)、恢復(fù)飽和系數(shù)；ω—泥沙沉速；C—Chezy系數(shù)；S_*—挾沙力，k、m—挾沙力系數(shù)、指數(shù)，γ_S、γ＇—泥沙容重與干容重；τ_*0—床面無(wú)量綱水流切應(yīng)力，τ_c—無(wú)量綱謝爾茲臨界切應(yīng)力；τ₀—床面水流切應(yīng)力，τ_c—謝爾茲臨界切應(yīng)力；d—床面分層粒徑。

3.1.2 數(shù)值計(jì)算格式

　　(1)離散網(wǎng)格及變量分布 模型采用非均勻網(wǎng)格，可以在研究量變化梯度較大的重要局部區(qū)域設(shè)置細(xì)密網(wǎng)格，在研究量變化較平緩的非主要區(qū)域設(shè)置較稀疏的網(wǎng)格。在離散網(wǎng)格上，標(biāo)量(ξ,H，S)被安排在單元中央，矢量(速度u,v)安排在單元的四周，交錯(cuò)網(wǎng)格上物理變量的位置相互錯(cuò)開(kāi)。分別用四個(gè)一維數(shù)組(xu(),yv(),xh(),yh())來(lái)確定水位、各流速分量的坐標(biāo) ^[３]。

　　(2)離散格式 本數(shù)學(xué)模型采用較成熟的“交替方向隱式差分逐行求解”方法，其特點(diǎn)是：將時(shí)間步長(zhǎng)(TIME STEP)分成前后兩個(gè)半步，在前半個(gè)時(shí)間步長(zhǎng)取某一個(gè)方向?yàn)殡[式；為保持對(duì)稱，在后半個(gè)時(shí)間步長(zhǎng)改變隱式方向。每一個(gè)時(shí)間步長(zhǎng)，都這樣交替的改變隱式方向計(jì)算。在前后兩個(gè)時(shí)間半步對(duì)控制方程進(jìn)行離散，為了物理概念上的清晰和格式的穩(wěn)定有效，引進(jìn)控制體同時(shí)在進(jìn)行對(duì)流項(xiàng)離散時(shí)引起迎風(fēng)格式。在前半個(gè)時(shí)間步長(zhǎng)，將連續(xù)方程與Y向動(dòng)量方程聯(lián)立，對(duì)u,ξ進(jìn)行隱式求解；在得到水流條件后隱式求解關(guān)于含沙濃度的傳移輸運(yùn)方程。在后半個(gè)時(shí)間步長(zhǎng)，將連續(xù)方程與Y向動(dòng)量方程聯(lián)立，對(duì)ν,ξ進(jìn)行隱式求解；也在獲得水流條件以后隱式求解泥沙傳移輸運(yùn)方程。

3.1.3 初始條件及邊界條件

　　(1)初始條件 對(duì)于給定的計(jì)算區(qū)域，在時(shí)間t=0時(shí)，令：

�И�(2)開(kāi)邊界條件ξ(x,y,t)=ξ^opb(x,y,t)或Q=Q^opb(t) 以及Sξ(x,y,t)=S^opb(x,y,t)。

　　其中ξ^opb，Q^opb以及S^opb分別為開(kāi)邊界上已知的水位、流量以及含沙量，一般由實(shí)測(cè)水文資料確定。

　　對(duì)于具有寬灘深槽復(fù)雜地形的河道，模型進(jìn)口開(kāi)邊界上宜給定洪水水位過(guò)程或流量過(guò)程，便于處理全斷面流量分配。上游開(kāi)邊界所需相應(yīng)的懸移質(zhì)含沙量及過(guò)程如無(wú)法實(shí)測(cè)資料，則借用附近河流相應(yīng)洪水的含沙量及過(guò)程。下游開(kāi)邊界處的含沙量可按第二類邊界條件確定。陸地邊界即河流岸邊，滿足固壁非穿越及無(wú)滑移條件，邊界上的法向、切向流速為零。

　　(3)動(dòng)邊界處理 在寬淺河流中的邊灘和江心洲隨著水位的變化，其水邊線也不斷改變，形成所謂動(dòng)邊界。模型相應(yīng)采用水邊界全區(qū)自動(dòng)跟蹤的處理方法，設(shè)置一個(gè)跟蹤指標(biāo)數(shù)組Iwet()。先將最大可能的淹沒(méi)區(qū)域包納在計(jì)算域區(qū)，在計(jì)算過(guò)程中根據(jù)計(jì)算單元內(nèi)的水深來(lái)判斷該單元是淹沒(méi)或是露出，即Iwet()應(yīng)該賦0還是賦1(0表示露出，為陸地；1表示淹沒(méi)，為應(yīng)計(jì)算的水域)，凡是陸地單元均不納入計(jì)算范圍。這種動(dòng)邊界智能跟蹤的處理方法特別適合邊灘和江心洲淹沒(méi)與出露頻繁的水域，使流場(chǎng)邊界條件模擬更為真實(shí)。

3.2

3.2.1 數(shù)學(xué)模型范圍及網(wǎng)格剖分

　　根據(jù)數(shù)值模擬的研究任務(wù)和交叉工程附近河段的河道形態(tài)，河床組成及洪痕調(diào)查情況，充分考慮上、下游邊界的水流條件，確定模型上邊界取在距交叉工程上游4.0km處；模型下邊界取在距交叉工程下游4.5km處，沿河道模型總長(zhǎng)約8.5km。側(cè)邊界在考慮洪水的最大可能淹沒(méi)范圍及地形條件，基本沿陡坎近岸高地選取。模型最大寬度為2.5km，計(jì)算區(qū)域總面積為21.25km^２。

　　為了反映交叉工程河段河道地形，適應(yīng)不同區(qū)域流場(chǎng)計(jì)算精度要求，交叉工程附近的重點(diǎn)區(qū)域，采用加密網(wǎng)格為50×50m;對(duì)于一般非重點(diǎn)計(jì)算區(qū)域，網(wǎng)格尺寸為100×50m。模型全部節(jié)點(diǎn)為4500個(gè)。

3.2.2 數(shù)學(xué)模型調(diào)試與驗(yàn)證

　　(1)調(diào)試依據(jù)與要求 模型調(diào)試主要依據(jù)是1963年8月和1996年8月該河發(fā)生的兩場(chǎng)重要洪水的調(diào)查洪痕(以下簡(jiǎn)稱“63.8”洪水、“96.8”洪水)；水文分析所得計(jì)算河段進(jìn)出口附近斷面的水位流量關(guān)系以及河床沖淤變形有關(guān)的查勘資料與河床探測(cè)資料。河道地形采用1966年1/10000地形圖并參考1994、1996兩年汛后河道的縱橫斷面測(cè)量資料�？紤]“96.8”洪水河道側(cè)侵蝕嚴(yán)重的特點(diǎn)，對(duì)河道地形進(jìn)行了適當(dāng)修正。“63.8”洪水屬特大洪水，“96.8”洪水屬一般大洪水，分別選用這兩場(chǎng)洪水驗(yàn)證模型流場(chǎng)，可使模型能正確模擬不同洪水條件，不同阻力特征的流場(chǎng)。

　　(2)調(diào)試結(jié)果 二維水沙數(shù)學(xué)模型反演“63.8”、.96.8”洪水的計(jì)算成果經(jīng)后處理，繪制成流場(chǎng)流速矢量圖。流場(chǎng)矢量圖所反映的流場(chǎng)的主流走向、流態(tài)基本與調(diào)查情況一致。這里僅給出“96.8”洪水流場(chǎng)中主要洪痕點(diǎn)水位調(diào)查值與計(jì)算值的對(duì)比，見(jiàn)圖2。河床糙率上段取0.030～0.035，下段取0.025～0.032。

　　通過(guò)對(duì)洪水流場(chǎng)主流走向、流態(tài)、流勢(shì)、流速分布等方面模擬與洪痕點(diǎn)水位校驗(yàn)對(duì)比，數(shù)學(xué)模型反演的“63.8”和“96.8”洪水與實(shí)際洪水調(diào)查情況基本一致，正確反映了兩類洪水的洪水河勢(shì)及行洪特點(diǎn)。數(shù)學(xué)模型在河道地形處理，糙率選擇，沖淤?gòu)?qiáng)度控制參數(shù)方面的選擇是合理的。

3.3

　　本次計(jì)算根據(jù)交叉工程設(shè)計(jì)方案，對(duì)設(shè)計(jì)洪水進(jìn)行了數(shù)值模擬，然后根據(jù)流場(chǎng)狀況，進(jìn)一步提出修改工程方案并進(jìn)行數(shù)值模擬和分析對(duì)比。

3.3.1 自然狀態(tài)計(jì)算成果

　　在未修工程的現(xiàn)狀條件下，對(duì)設(shè)計(jì)洪水進(jìn)行了自然狀態(tài)河段的數(shù)值模擬。由洪峰期計(jì)算河段平面流速分布(圖3)可以看出，在設(shè)計(jì)洪水時(shí)，交叉工程斷面附近水流收縮集中，河槽沖刷嚴(yán)重。受江心洲影響，河道分成兩汊又漸次合一(北汊較大)，南岸高地下游有部分回流區(qū)。

3.3.3 修改方案計(jì)算成果

　　根據(jù)自然狀態(tài)和設(shè)計(jì)工程狀態(tài)下流場(chǎng)流態(tài)、壅水、沖刷深度及部位的計(jì)算成果分析，得知南岸灘地水流流速不高，近期調(diào)查南槽僅略有增大趨勢(shì)，在此適當(dāng)約束水流不會(huì)產(chǎn)生很大影響。因此可以考慮進(jìn)行修改、優(yōu)化工程布置，這里提出了分別在南、北岸倒虹吸工程進(jìn)口、出處將倒虹吸段縮窄50m、和將倒虹吸工程全部安排在北槽的三種修改方案，相應(yīng)倒虹吸工程段長(zhǎng)度均為650m。對(duì)不同修改方案在設(shè)計(jì)洪水條件下進(jìn)行了數(shù)值模擬計(jì)算，不同方案時(shí)在交叉工程附近引起的壅水程度、沖刷程度及發(fā)展規(guī)律各不相同。交叉工程斷面處的水位、流速、河床變形可見(jiàn)表1、圖6。

Characteristic values of cross-section in different modified alternatives

3.3.4 修改方案成果分析

　　由于受河道地形的影響，在設(shè)計(jì)方案和修改方案中，洪水主流略偏靠北側(cè)深槽。受交叉建筑物擠壓，南北兩汊流速均較大，北汊最大單寬流量可達(dá)15m²/s以上，應(yīng)注意邊岸工程的防護(hù)。修改方案中壓縮段基本屬于近灘范圍，故這對(duì)倒虹吸工程處河道實(shí)際主流過(guò)水影響不大。交叉工程斷面附近壅水最大高度(修改方案1)為0.40m，壅水最大范圍約1000m。由于流速較大，所以交叉工程斷面在不同流量級(jí)洪水時(shí)，均有主槽全斷面均有沖刷，而兩側(cè)灘地則有沖有淤。北汊河槽主要是河床下切，落峰時(shí)略有回淤；南槽除河床下切外還有側(cè)向侵蝕，河床被不斷沖深、拓寬。洪水過(guò)程中河床一直沖刷，但洪峰后期沖刷速度轉(zhuǎn)緩。

　　由于河床變形劇烈，洪水演進(jìn)中一般Z_max出現(xiàn)較早，而此時(shí)流量Q未達(dá)到最大值，沖刷強(qiáng)度還繼續(xù)增加，而后次第出現(xiàn)Q_max、ν_max，此時(shí)沖刷強(qiáng)度達(dá)到最大，此后沖刷雖繼續(xù)進(jìn)行但發(fā)展速率逐漸減小。在洪水后期，全斷面便開(kāi)始由沖轉(zhuǎn)變?yōu)槲⒂佟?/p>

　　比較不同方案的計(jì)算結(jié)果，可以看出修改方案1既可改善北側(cè)水流條件又適當(dāng)考慮保留南槽過(guò)流能力，能更好適應(yīng)上游河勢(shì)改變可能引起的主流調(diào)整。

4 結(jié)語(yǔ)

　　1.采用平面二維水沙數(shù)學(xué)模型研究七里河交叉工程河段洪水及河床變形問(wèn)題是可以滿足工程要求的。

　　2.通過(guò)歷史洪水的反演與洪痕校核，進(jìn)行了模型參數(shù)的調(diào)試和糙率調(diào)整，計(jì)算表明這對(duì)于正確模擬和預(yù)測(cè)不同條件洪水與河床變形是必須的。

　　3.根據(jù)設(shè)計(jì)洪水條件，對(duì)不同工程方案進(jìn)行了模擬計(jì)算，所提出的有關(guān)流場(chǎng)流態(tài)、主流走向、壅水、沖刷狀況，交叉工程斷面處壅水沖刷的有關(guān)特征值可做為交叉建筑物初步設(shè)計(jì)的依據(jù)。

　　4.根據(jù)不同工程方案的流場(chǎng)流態(tài)、流速及工程斷面處的壅水、沖刷等狀況的分析，建議采用修改方案1(北側(cè)倒虹吸縮短50m)并適當(dāng)布設(shè)導(dǎo)流堤等邊岸防沖措施。最佳方案的選定還須進(jìn)行全面的技術(shù)經(jīng)濟(jì)比較。

[1] 曹祖德，王運(yùn)洪。水動(dòng)力泥沙數(shù)值模擬。天津大學(xué)出版社，1994.

[2] T.J.Chung. Finite Element Analysis G-H International Book Company,1978.

[3] 孫東坡等。河流模擬理論與實(shí)踐。武漢水利電力大學(xué)出版社，1998，40