非恒流量計算中的非恒流量分析

非恒流量計算中的非恒流量分析

由于河流網(wǎng)絡(luò)的河流垂直和水平方向，水流方向不能確定，因此有可能有門、壩等建筑物。當(dāng)研究領(lǐng)域較大時，利用物理模型難以控制和調(diào)節(jié)由于水和土壤變化而發(fā)生的變化。數(shù)學(xué)模型已成為研究和模擬河流網(wǎng)絡(luò)水和沉積物問題的唯一手段。盡管已經(jīng)簡化為一維問題,但由于在分汊點處要考慮水流、泥沙的銜接,而且復(fù)雜河網(wǎng)中大量的河段和汊點數(shù)目將引起計算速度和精度下降,這些都增加了河網(wǎng)水沙模擬的復(fù)雜性。目前,隱式解法中分級、分組解法以其穩(wěn)定性好、收斂快等優(yōu)點,已成為實際中求解河網(wǎng)非恒定流方程組的主要方法。河網(wǎng)水沙運動模擬相比于單一河道水沙模擬的主要困難表現(xiàn)于,要處理好汊點處的沙量連續(xù),即分沙問題,還要解決汊點數(shù)量多而引起的計算速度和精度問題。多位學(xué)者提出了單一河道分汊口的懸移質(zhì)分沙模式,用于分汊口的分沙特性分析,但能夠用于大型河網(wǎng)泥沙數(shù)值模擬的還不多見。目前的河網(wǎng)水沙模型多是以汊點處的分流比直接確定分沙比,顯得過于粗糙。此外,許多模型采取逐斷面求解離散后的不平衡輸沙方程,類似于水流方程的直接解法,必然是以速度和精度的降低為代價。針對以上問題,本文分析提出了多種分沙模式用于河網(wǎng)水沙模擬,并且實現(xiàn)了泥沙方程組的分級解法以提高計算效率和精度,對河網(wǎng)中經(jīng)常存在的可儲水汊點處的水沙分匯問題也進行了初步的分析研究。1建立非固定流水?dāng)?shù)學(xué)模型1.1基本控制方程式與點連接條件(1)泥沙運動方程及含沙量河網(wǎng)區(qū)一般位于中下游的平原或河口三角洲區(qū),水體所挾帶泥沙以細顆粒懸移質(zhì)為主,河床沖淤也以懸移質(zhì)起主要作用,因此,本文在模型中也只考慮懸移質(zhì)泥沙。通常單一河道懸沙數(shù)學(xué)模型基本控制方程可表述如下:水流連續(xù)方程?Q?x+B?Ζ?t=ql(1)水流運動方程?Q?t+??x(βQ2A)+gA(?Ζ?x)+gn2Q|Q|AR4/3=0(2)懸移質(zhì)泥沙運動方程?AS?t+?QS?x=-αBω(S-S*)(3)河床變形方程?AS?t+?QS?x+ρ′B?Ζ0?t=0(4)水流挾沙力采用張瑞瑾公式S*=Κ(U3gRω)m(5)式中x為流程;Q為流量;ql為旁側(cè)單寬入流量;Z為水位;B為河寬;t為時間;β為動量修正系數(shù);A為過水?dāng)嗝婷娣e;R為水力半徑;Z0為平均河底高程;U、H分別為斷面平均流速及水深;ω為懸移質(zhì)沉速;ρ′為泥沙干密度;α為懸移質(zhì)泥沙恢復(fù)飽和系數(shù);S、S*分別為斷面平均含沙量及挾沙力;K、m等系數(shù)需要根據(jù)實測資料率定。(2)河網(wǎng)點與單一河道的銜接河網(wǎng)中,汊點是兩條或兩條以上的河段交匯點,計算中,邊界點也可以作為汊點處理。河網(wǎng)與單一河道的不同之處在于其中分布著大量汊點,各單一河段通過汊點交換水沙和能量,水流、泥沙必須滿足質(zhì)量和能量的守恒,所以對于河網(wǎng)中的汊點(包括內(nèi)部汊點和邊界汊點),水流、泥沙必須滿足一定的銜接條件。qn+dm的增殖比例關(guān)系進出每一汊點的流量必須與汊點實際水量增減率相平衡:Qn+1m+L(m)∑l=1Qn+1m?l=?Vm?tm=1?2???Μ(6)將式(6)寫成增量形式,并將水量增減率寫成差分形式:Qn+1m+L(m)∑l=1Qnm?l+L(m)∑l=1ΔQm?l=SmΔΖmΔtm=1?2???Μ(7)式中M為河網(wǎng)中的汊點數(shù);L(m)為與汊點m相連接的河段數(shù);Qn+1m?l為與汊點m相接的第l條河段流進(或流出)該汊點的流量,Qn+1m為除與汊點m相接的河段外的其它流進(或流出)該汊點的流量;Vm為汊點m的蓄水量;Sm為汊點m的水面面積;ΔZm為汊點m的水位增量。各河段水位增加壓計算汊點動量守恒條件與是否考慮汊點處各河段端點處流速水頭、阻力損失等有關(guān),一般情況下,可近似地認為汊點處各河段端點水位相同。即Ζm?1=Ζm?2=?=Ζm?L(m)=Ζmm=1?2???Μ(8)同樣可以近似地認為汊點處各河段水位增量也相同。ΔΖm?1=ΔΖm?2=?=ΔΖm?L(m)=ΔΖmm=1?2???Μ(9)入pa和bmmQn+1mSn+1m+L(m)∑l=1Qn+1m?l?Sn+1m?l=ρ′Am?Ζb?m?tm=1?2???Μ(10)式中Sn+1m?l為與流量Qn+1m?l相對應(yīng)的含沙量;Sn+1m為與Qn+1m相應(yīng)的含沙量;?Ζb?m?t為汊點m的河床高程變形情況;ρ′為泥沙干密度;Am為汊點m的平面面積。一般情況下,汊點的沖淤可以忽略,式(10)右邊一項可略去。根據(jù)實際水流方向區(qū)別進入和流出的沙量,式(10)化簡為Qn+1mSn+1m+Lin(m)∑l=1Qn+1m?l?Sn+1m?l=Lout(m)∑l=1Qn+1m?l?Sn+1m?l(11)式中Lin(m)、Lout(m)分別為匯入河段數(shù)和流出河段數(shù);流量按照實際的水流方向取值,均為正值。1.2點方程的轉(zhuǎn)化由式(1)和式(2)組成的方程采用任何一種隱格式離散,結(jié)合水流銜接條件式(7)、式(9)和適當(dāng)?shù)倪吔鐥l件即構(gòu)成封閉的河網(wǎng)方程組,一般采用三級解法求解,其基本過程是通過迭代消元將未知數(shù)集中于汊點上,得到上下端點的流量增量與上下端點水位增量之間的關(guān)系:ΔQ1=E1ΔΖ1+F1+Η1ΔΖΙ(l)(12)ΔQΙ(l)=E′1ΔΖ1+F′1+Η′1ΔΖΙ(l)(13)式中1和I(l)分別為第l條河段上下端點(本文規(guī)定1斷面為下游出口,I(l)斷面為上游進口,下同);式中系數(shù)只和各斷面上一時段水力因素有關(guān)。代入汊點連接方程之后,河網(wǎng)方程組可以轉(zhuǎn)化成汊點方程組:[A]{ΔΖ}={B}(14)式中A為系數(shù)矩陣,各元素僅與系數(shù)E1、H1、E′1、H′1及汊點水面面積S/Δt有關(guān);{ΔZ}為汊點水位增量矢量;{B}為汊點處各河段端點流量Qnm?l,其他入流量Qn+1m及F1、F1′組成的矢量。式(14)的階數(shù)與河網(wǎng)中的汊點數(shù)相同,但若汊點很多時,汊點方程組的系數(shù)矩陣也很大,給求解帶來了困難。針對這個問題,文獻提出了汊點分組解法,基本思路是利用矩陣分塊,根據(jù)河網(wǎng)的實際情況和計算的需要,將河網(wǎng)中的汊點分為多組,能夠有效地降低系數(shù)矩陣的階數(shù)、節(jié)省存儲量,從而提高計算的速度和精度。利用分級解法或汊點分組解法確定出各汊點的水位增量,逐步回代可得河段端點流量增量以及各河段內(nèi)部斷面的水位、流量增量,這樣就可進行下次迭代或下個時段的計算。1.3河流網(wǎng)絡(luò)沙模型的構(gòu)建(1)有限分析法式(3)對數(shù)值解法的精度要求較高,采用一般的差分格式,精度不能保證,含沙量的計算結(jié)果有時可能會出現(xiàn)不合理的數(shù)值,甚至為負值;采用高精度的數(shù)值格式,不僅計算復(fù)雜,而且計算量太大,不適于河網(wǎng)計算。為了避免含沙量計算結(jié)果的不合理現(xiàn)象,同時不致于使計算量猛增,對懸移質(zhì)泥沙運動方程的離散可借鑒有限分析法的思想,相鄰時層之間用差分法求解,在同一時層上求式(3)的分析解,求解關(guān)于S的常微分方程可得:Sn+1i=Sn+1i+1e-[(αωˉq)n+1+1ˉUn+1Δt]Δxi+αωˉUn+1Δt?ˉSn+1*+ˉqn+1ˉSnαωˉUn+1Δt+ˉqn+1{1-e-[(αωˉq)n+1+1ˉUn+1Δt]Δxi}(15)式中ˉU、ˉq、ˉS*及ˉS分別為Δxi河段內(nèi)的平均流速、單寬流量、挾沙力及含沙量。由式(15)可見,在水力要素已知的情況下,沿著水流方向,Δxi河段出口斷面的含沙量決定于進口斷面的含沙量,以及兩斷面上一時層的含沙量。(2)遞推關(guān)系的分析一般的水沙數(shù)學(xué)模型采用非耦合解法,對單一河道而言是先通過非恒定流計算得出各個斷面的水力要素和挾沙力,然后就可利用上游含沙量從上至下依次推求至下游出口斷面。但對于河網(wǎng),逐河段逐斷面進行計算就顯得過于繁瑣。因為對某一汊點而言,與其相連的既有匯入河段,也有流出河段,而根據(jù)沙量的守恒,必須計算出所有匯入河段出口斷面的沙量之后,才能進行各輸出河段的計算。逐河段計算的實現(xiàn)過程就是不斷地枚舉、判斷,直至最終實現(xiàn)所有河段的遍歷,效率不高且過于繁瑣。注意到式(15)中,下游斷面僅對上游相鄰斷面有依賴性,如果能夠仿照水流分級解法的思想,利用遞推關(guān)系消去河段中間各斷面的未知量,使河段進口和出口斷面直接發(fā)生聯(lián)系,無疑將使問題簡化。假定河段實際的水流方向與上文的規(guī)定相同,式(15)可以表示為Sn+1i=Ai?Sn+1i+1+Bii=Ι(l)-1???1(16)其中Ai=e-[(αωˉq)n+1+1ˉUn+1Δt]ΔxiBi=αωˉUn+1ΔtˉSn+1*+ˉqn+1ˉSnαωˉUn+1Δt+ˉqn+1{1-e-[(αωˉq)n+1+1ˉUn+1Δt]Δxi}i=Ι(l)-1???1進行完一個時段的非恒定流計算后,各水力要素已知,式(16)中系數(shù)將為已知值。假如對I(l)以下任意斷面i,存在以下遞推關(guān)系:Sn+1i=Ρi?Sn+1Ι(l)+Ri(17)將式(17)代入式(16)可得:Sn+1i-1=Ρi-1Sn+1Ι(l)+Ri-1(18)式中Pi-1=Ai-1Pi,Ri-1=Ai-1Ri+Bi-1。如此可以遞推至出口斷面,得Sn+11=Ρ1Sn+1Ι(l)+R1(19)式中P1=A1P2,R1=A1R2+B1。事實上,式(16)、式(17)形式相同,在上游端這種遞推關(guān)系是自然成立的,對I(l)、I(l)-1微段只需令PI(l)=AI(l),RI(l)=BI(l)即可,所以,式(19)成立,其中系數(shù)可沿水流方向自I(l)斷面至1斷面遞推得到。同理,若計算結(jié)果為逆流,即1斷面為上游,同樣可以得到與式(19)形式相同的關(guān)系:Sn+1Ι(l)=Ρ′Ι(l)Sn+11+R′Ι(l)(20)式(20)中系數(shù)沿水流方向自1斷面至I(l)斷面依次遞推得到。如此,依照計算出的水流方向,通過式(19)或式(20)可將河段末端的含沙量直接與上游端的含沙量建立聯(lián)系,與水流的分級解法類似,未知數(shù)被集中于汊點上。(3)分配復(fù)雜，難度汊點是河網(wǎng)中水流、泥沙的分匯點,流出汊點的河段之間應(yīng)存在著水量、沙量的分配比例問題。一般來講,在糙率確定的情況下,分流比可以通過水流模擬得到,但沙量的分配要復(fù)雜得多。假設(shè)通過匯流河段進入汊點的水沙混合均勻,其含沙量為Sm,pj,各分流口門的含沙量為Km,lSm,pj,Km,l為第l個分流口門的含沙量與汊點平均含沙量比值,則式(11)可表示為Lin(m)∑l=1Qn+1m?l?Sn+1m?l=Sm?pjLout(m)∑l=1Qn+1m?l?Κn+1m?l(21)對河網(wǎng)中的M個汊點可以列出M個方程,但未知數(shù)卻有河段數(shù)L個,所以方程組不封閉,必須補充一定的分沙模式,確定出Km,l,使得方程組得以封閉。封閉型點方程組比較簡單一點的處理辦法是,認為各分流口門含沙量就是汊點平均含沙量,即Km,l=1。這樣根據(jù)式(21)列出M個方程,其未知數(shù)降為M個,聯(lián)系式(19)或式(20)所示的關(guān)系,可形成封閉的汊點方程組為[As]{Spj}={Bs}(22)式中[As]是M×M維系數(shù)矩陣,其中元素均已知;{Bs}是已知向量。以上處理辦法是一種比較粗糙的方法,實際中由于各分流河段的比降、糙率、口門形態(tài)、口門高程、口門附近的河床形態(tài)等不同,其水力因素各異,挾沙能力也不相同。因此,有必要采用更合理的分沙模式。lout條分流河段對單一河道的分汊問題,有的研究提出了考慮地形因素的分沙模式,丁君松等認為主支汊河床上均存在一鞍點,造成主支汊引水深度存在差別,由于含沙量及泥沙級配沿垂線分布的不均勻,造成主支汊含沙量不同。將這一模型應(yīng)用于河網(wǎng)汊點,假設(shè)沿垂線含沙量分布符合張瑞謹公式:S=Spjc(1+c)(c+ξ)2(23)式中Spj為垂線平均含沙量;c為含沙量分布不均勻程度的數(shù)值,恒為正,c值越小,分布越不均勻;ξ為相對水深,河底為0,水面為1;S為ξ點處的含沙量(全沙、分組沙均可)。假設(shè)第m汊點的所有Lout條分流河段中,第1條河段進口口門高程最低,其參考點為n1,n1點至水面的相對水深為1,此進口含沙量為Sm,1,其余分流河段中任一河段的口門高程參考點為ni,ni相對n1的相對水深為ξi,則可以導(dǎo)出其進口含沙量Sm,i與Sm,1的比值為Κm?i=Sm?iSm?1=cc+ξii=2???Lout(m)(24)這樣,可以得到所有的Lout條分流河段入口含沙量之間的比例符合下式:S1∶那么,根據(jù)地形確定出各分流河段的相對水深ξi之后,式(21)轉(zhuǎn)化為[As]{S1}={Bs}(26)式中{S1}為各汊點的Sm,1組成的矢量,但[As]中的系數(shù)與式(22)不同,將與Km,i有關(guān)。實際上,此模式運用的關(guān)鍵是出流河段口門高程參考點的確定,可取進口斷面的平均高程,也可取距入口一定距離之內(nèi)的河床最高點。此種模式對資料要求較高,當(dāng)各口門高程比較接近時,式(24)中的比值將接近于1,與簡單方式是相同的。分流河段分沙比的物理機理式(24)結(jié)構(gòu)比較繁雜,而且每次沖淤計算之后都要重新對其修正。河網(wǎng)計算涉及區(qū)域大,地形等資料收集不易,而且河網(wǎng)中的汊點可能連接多個分流河道,出于計算速度和精度考慮,要求提出一個能夠反映物理實質(zhì),簡潔而又能滿足實際精度需要的分沙模式。造成各分流口門含沙量不同的主要原因是水流條件不同,而反映一定條件下水流挾沙能力的指標(biāo)為S*,所以,影響分沙比的各種因素的綜合作用很大程度就體現(xiàn)在S*的大小上。出于這種考慮,本文提出以分流河段進口斷面的挾沙力S*值確定汊點的分沙比,認為存在下式的關(guān)系,這在物理機理上是合理的。S1∶S2∶?∶SLout=S*1∶S*2∶?∶S*Lout(27)根據(jù)式(27),在每個汊點m均存在一個比例常數(shù)Nm,使得Si=Νm?S*ii=1?2???Ιout(l)(28)當(dāng)然,在不同的時間Nm是變化的。這樣,式(21)汊點方程組轉(zhuǎn)化為[As]{Κ}={Bs}(29)式中矩陣[As]中元素只與流量和挾沙力有關(guān);{Bs}為已知向量。[As]和{Bs}由式(19)或式(20)結(jié)合式(21)得到,{K}為未知向量(各個節(jié)點的Nm值)。在挾沙力公式中,含沙量與流速的高次方成正比,此種模式實際是以流速為主分配沙量,當(dāng)各分流口門流速相近時,就變成了簡單分沙模式。1.4點方程組求解不論采取哪種分沙模式,汊點方程組式(22)、式(26)或式(29)結(jié)構(gòu)類似,與水流求解過程中的汊點方程組也類似,所以可采取相同的解法。當(dāng)汊點數(shù)較多時,同樣可采取分組的解法,重新利用求解水流方程過程中的分組情況,為計算程序的編制提供方便。求解汊點泥沙方程組之后,通過逐步回代,可以得到所有斷面的含沙量Sn+1。2河網(wǎng)沉積物模型的適用性2.1非均勻沙模型懸移質(zhì)泥沙方程的求解并未區(qū)分均勻沙和非均勻沙,可將其當(dāng)作均勻沙,也可以是非均勻沙中的某一粒徑組。如果為非均勻沙模型,對每一粒徑組均可列出相應(yīng)的汊點方程組,而河床的總變形將是每一粒徑組所引起的變形之和。除此之外,還將涉及非均勻沙挾沙力及其級配的計算、活動層選取、床沙級配調(diào)整、床沙粗化或細化過程中糙率的調(diào)整等系列問題。2.2各口門處的截沙模式:丁君松模式以上給出的三種分沙模式可以看出,第一種簡單分沙模式是丁君松模式和挾沙力模式在特殊情況下的特例,總體來講處理比較粗糙。當(dāng)對泥沙模擬的精度要求不高,或者當(dāng)各口門水力條件相近,口門高程相差不大時,可以考慮采用此種模式。丁君松模式注意了含沙量沿垂線分布的不均勻和各分流口門地形(主要是高程)不同造成的含沙量不同,而挾沙力模式注重了水力條件不同造成各口門處的挾沙力區(qū)別。這兩種方式分別注重了問題的兩個方面。天然情況下,水流、泥沙運動相互影響,相互耦合,這是河流發(fā)生自調(diào)整作用的原因,河網(wǎng)也不例外。設(shè)想河網(wǎng)中連接某一汊點的分流河段進口口門高程、形態(tài)類似,而河段的河床比降、糙率不同,因而水力要素不同。按照推廣后的丁君松模式,它們的進口含沙量相同,根據(jù)河床的自調(diào)整原理,各河段將發(fā)生沖刷或淤積來調(diào)整河床形態(tài),此過程中同時也淤高或蝕低了口門高程,直至沖淤后口門進沙量與水流相適應(yīng)。從這個角度講,挾沙力模式與此是一致的,因為河床調(diào)整總是以加大或減小比降達到減小含沙量與挾沙力的差值為目標(biāo),而調(diào)整過程比較明顯地表現(xiàn)在河床上段的泥沙沖淤上。只不過丁君松模式更注重了水深對泥沙粒徑的分選作用,當(dāng)各分流口門中存在一個河底高程明顯較低的主汊時(例如長江三口),此種模式更適用,這也與它是從單一河道的分汊情況推廣過來的這一背景相符合。實際應(yīng)用中,還可以根據(jù)分汊的夾角等因素對丁君松模式進行修正。實際上,式(26)和式(29)形式相同,在同一模型中可以根據(jù)汊點的實際情況靈活選用兩種模式的任一種,只是相關(guān)系數(shù)存在差別,不影響汊點矩陣的形成和求解。2.3水體沙量計算實際的河網(wǎng)計算中,經(jīng)常有可蓄水汊點,如湖泊、蓄洪區(qū)等,只要提供水位與水面面積(容積)曲線,仍可用汊點水流連續(xù)方程(7)計算,但含沙量需要單獨處理。進水情況,含沙量可按以上的模式計算。進入汊點的泥沙將發(fā)生沉降,沉降的時間將影響流出汊點的含沙量。設(shè)泥沙沉速為ω,汊點水位為Z,出流口門高程為Zd,hd=Z-Zd。泥沙降到口門以下需要時間Tω=hd/ω,當(dāng)發(fā)生水體在汊點中貯留時間超過Tω時,可近似認為出流含沙量為0。當(dāng)T<Tω時,泥沙經(jīng)過時間T下降距離d=ω?T,則水面下深度為d的范圍內(nèi)含沙量近似為0,再以下的含沙量沿垂線分布符合張瑞謹公式(23)。從Zd向水面積分,可得出流平均含沙量。以上處理辦法既適用于均勻沙也可以是非均勻沙的某一組。3蓄滯洪區(qū)一維河網(wǎng)模型本文采用杜家臺蓄滯洪區(qū)的實測資料對所建模型進行驗證和分析。杜家臺分洪閘位于漢江下游仙桃市區(qū)以東,上距丹江口水庫502km,下距漢江河口150km。在漢江干流仙桃以下河段超過安全泄量時(在5250～9150m3/s之間,根據(jù)漢口水位高低來確定),超額洪水經(jīng)過杜家臺分洪閘分泄,由分洪道引入漢南泛區(qū)調(diào)蓄,再經(jīng)黃陵磯閘泄入長江。杜家臺總共有11個蓄滯洪區(qū),包括下東城垸、上東城垸、新農(nóng)垸、桐湖斗昌湖、洪北東垸、洪北西垸、曲口垸、興無垸、銀蓮湖、消泗外垸以及未圍部