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

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

由于河流网络的河流垂直和水平方向，水流方向不能确定，因此有可能有门、坝等建筑物。当研究领域较大时，利用物理模型难以控制和调节由于水和土壤变化而发生的变化。数学模型已成为研究和模拟河流网络水和沉积物问题的唯一手段。尽管已经简化为一维问题,但由于在分汊点处要考虑水流、泥沙的衔接,而且复杂河网中大量的河段和汊点数目将引起计算速度和精度下降,这些都增加了河网水沙模拟的复杂性。目前,隐式解法中分级、分组解法以其稳定性好、收敛快等优点,已成为实际中求解河网非恒定流方程组的主要方法。河网水沙运动模拟相比于单一河道水沙模拟的主要困难表现于,要处理好汊点处的沙量连续,即分沙问题,还要解决汊点数量多而引起的计算速度和精度问题。多位学者提出了单一河道分汊口的悬移质分沙模式,用于分汊口的分沙特性分析,但能够用于大型河网泥沙数值模拟的还不多见。目前的河网水沙模型多是以汊点处的分流比直接确定分沙比,显得过于粗糙。此外,许多模型采取逐断面求解离散后的不平衡输沙方程,类似于水流方程的直接解法,必然是以速度和精度的降低为代价。针对以上问题,本文分析提出了多种分沙模式用于河网水沙模拟,并且实现了泥沙方程组的分级解法以提高计算效率和精度,对河网中经常存在的可储水汊点处的水沙分汇问题也进行了初步的分析研究。1建立非固定流水数学模型1.1基本控制方程式与点连接条件(1)泥沙运动方程及含沙量河网区一般位于中下游的平原或河口三角洲区,水体所挟带泥沙以细颗粒悬移质为主,河床冲淤也以悬移质起主要作用,因此,本文在模型中也只考虑悬移质泥沙。通常单一河道悬沙数学模型基本控制方程可表述如下:水流连续方程∂Q∂x+B∂Ζ∂t=ql(1)水流运动方程∂Q∂t+∂∂x(βQ2A)+gA(∂Ζ∂x)+gn2Q|Q|AR4/3=0(2)悬移质泥沙运动方程∂AS∂t+∂QS∂x=-αBω(S-S*)(3)河床变形方程∂AS∂t+∂QS∂x+ρ′B∂Ζ0∂t=0(4)水流挟沙力采用张瑞瑾公式S*=Κ(U3gRω)m(5)式中x为流程;Q为流量;ql为旁侧单宽入流量;Z为水位;B为河宽;t为时间;β为动量修正系数;A为过水断面面积;R为水力半径;Z0为平均河底高程;U、H分别为断面平均流速及水深;ω为悬移质沉速;ρ′为泥沙干密度;α为悬移质泥沙恢复饱和系数;S、S*分别为断面平均含沙量及挟沙力;K、m等系数需要根据实测资料率定。(2)河网点与单一河道的衔接河网中,汊点是两条或两条以上的河段交汇点,计算中,边界点也可以作为汊点处理。河网与单一河道的不同之处在于其中分布着大量汊点,各单一河段通过汊点交换水沙和能量,水流、泥沙必须满足质量和能量的守恒,所以对于河网中的汊点(包括内部汊点和边界汊点),水流、泥沙必须满足一定的衔接条件。qn+dm的增殖比例关系进出每一汊点的流量必须与汊点实际水量增减率相平衡: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为河网中的汊点数;L(m)为与汊点m相连接的河段数;Qn+1m‚l为与汊点m相接的第l条河段流进(或流出)该汊点的流量,Qn+1m为除与汊点m相接的河段外的其它流进(或流出)该汊点的流量;Vm为汊点m的蓄水量;Sm为汊点m的水面面积;ΔZm为汊点m的水位增量。各河段水位增加压计算汊点动量守恒条件与是否考虑汊点处各河段端点处流速水头、阻力损失等有关,一般情况下,可近似地认为汊点处各河段端点水位相同。即Ζ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相对应的含沙量;Sn+1m为与Qn+1m相应的含沙量;∂Ζb‚m∂t为汊点m的河床高程变形情况;ρ′为泥沙干密度;Am为汊点m的平面面积。一般情况下,汊点的冲淤可以忽略,式(10)右边一项可略去。根据实际水流方向区别进入和流出的沙量,式(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)分别为汇入河段数和流出河段数;流量按照实际的水流方向取值,均为正值。1.2点方程的转化由式(1)和式(2)组成的方程采用任何一种隐格式离散,结合水流衔接条件式(7)、式(9)和适当的边界条件即构成封闭的河网方程组,一般采用三级解法求解,其基本过程是通过迭代消元将未知数集中于汊点上,得到上下端点的流量增量与上下端点水位增量之间的关系:ΔQ1=E1ΔΖ1+F1+Η1ΔΖΙ(l)(12)ΔQΙ(l)=E′1ΔΖ1+F′1+Η′1ΔΖΙ(l)(13)式中1和I(l)分别为第l条河段上下端点(本文规定1断面为下游出口,I(l)断面为上游进口,下同);式中系数只和各断面上一时段水力因素有关。代入汊点连接方程之后,河网方程组可以转化成汊点方程组:[A]{ΔΖ}={B}(14)式中A为系数矩阵,各元素仅与系数E1、H1、E′1、H′1及汊点水面面积S/Δt有关;{ΔZ}为汊点水位增量矢量;{B}为汊点处各河段端点流量Qnm‚l,其他入流量Qn+1m及F1、F1′组成的矢量。式(14)的阶数与河网中的汊点数相同,但若汊点很多时,汊点方程组的系数矩阵也很大,给求解带来了困难。针对这个问题,文献提出了汊点分组解法,基本思路是利用矩阵分块,根据河网的实际情况和计算的需要,将河网中的汊点分为多组,能够有效地降低系数矩阵的阶数、节省存储量,从而提高计算的速度和精度。利用分级解法或汊点分组解法确定出各汊点的水位增量,逐步回代可得河段端点流量增量以及各河段内部断面的水位、流量增量,这样就可进行下次迭代或下个时段的计算。1.3河流网络沙模型的构建(1)有限分析法式(3)对数值解法的精度要求较高,采用一般的差分格式,精度不能保证,含沙量的计算结果有时可能会出现不合理的数值,甚至为负值;采用高精度的数值格式,不仅计算复杂,而且计算量太大,不适于河网计算。为了避免含沙量计算结果的不合理现象,同时不致于使计算量猛增,对悬移质泥沙运动方程的离散可借鉴有限分析法的思想,相邻时层之间用差分法求解,在同一时层上求式(3)的分析解,求解关于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河段内的平均流速、单宽流量、挟沙力及含沙量。由式(15)可见,在水力要素已知的情况下,沿着水流方向,Δxi河段出口断面的含沙量决定于进口断面的含沙量,以及两断面上一时层的含沙量。(2)递推关系的分析一般的水沙数学模型采用非耦合解法,对单一河道而言是先通过非恒定流计算得出各个断面的水力要素和挟沙力,然后就可利用上游含沙量从上至下依次推求至下游出口断面。但对于河网,逐河段逐断面进行计算就显得过于繁琐。因为对某一汊点而言,与其相连的既有汇入河段,也有流出河段,而根据沙量的守恒,必须计算出所有汇入河段出口断面的沙量之后,才能进行各输出河段的计算。逐河段计算的实现过程就是不断地枚举、判断,直至最终实现所有河段的遍历,效率不高且过于繁琐。注意到式(15)中,下游断面仅对上游相邻断面有依赖性,如果能够仿照水流分级解法的思想,利用递推关系消去河段中间各断面的未知量,使河段进口和出口断面直接发生联系,无疑将使问题简化。假定河段实际的水流方向与上文的规定相同,式(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)中系数将为已知值。假如对I(l)以下任意断面i,存在以下递推关系: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)形式相同,在上游端这种递推关系是自然成立的,对I(l)、I(l)-1微段只需令PI(l)=AI(l),RI(l)=BI(l)即可,所以,式(19)成立,其中系数可沿水流方向自I(l)断面至1断面递推得到。同理,若计算结果为逆流,即1断面为上游,同样可以得到与式(19)形式相同的关系:Sn+1Ι(l)=Ρ′Ι(l)Sn+11+R′Ι(l)(20)式(20)中系数沿水流方向自1断面至I(l)断面依次递推得到。如此,依照计算出的水流方向,通过式(19)或式(20)可将河段末端的含沙量直接与上游端的含沙量建立联系,与水流的分级解法类似,未知数被集中于汊点上。(3)分配复杂，难度汊点是河网中水流、泥沙的分汇点,流出汊点的河段之间应存在着水量、沙量的分配比例问题。一般来讲,在糙率确定的情况下,分流比可以通过水流模拟得到,但沙量的分配要复杂得多。假设通过汇流河段进入汊点的水沙混合均匀,其含沙量为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)对河网中的M个汊点可以列出M个方程,但未知数却有河段数L个,所以方程组不封闭,必须补充一定的分沙模式,确定出Km,l,使得方程组得以封闭。封闭型点方程组比较简单一点的处理办法是,认为各分流口门含沙量就是汊点平均含沙量,即Km,l=1。这样根据式(21)列出M个方程,其未知数降为M个,联系式(19)或式(20)所示的关系,可形成封闭的汊点方程组为[As]{Spj}={Bs}(22)式中[As]是M×M维系数矩阵,其中元素均已知;{Bs}是已知向量。以上处理办法是一种比较粗糙的方法,实际中由于各分流河段的比降、糙率、口门形态、口门高程、口门附近的河床形态等不同,其水力因素各异,挟沙能力也不相同。因此,有必要采用更合理的分沙模式。lout条分流河段对单一河道的分汊问题,有的研究提出了考虑地形因素的分沙模式,丁君松等认为主支汊河床上均存在一鞍点,造成主支汊引水深度存在差别,由于含沙量及泥沙级配沿垂线分布的不均匀,造成主支汊含沙量不同。将这一模型应用于河网汊点,假设沿垂线含沙量分布符合张瑞谨公式:S=Spjc(1+c)(c+ξ)2(23)式中Spj为垂线平均含沙量;c为含沙量分布不均匀程度的数值,恒为正,c值越小,分布越不均匀;ξ为相对水深,河底为0,水面为1;S为ξ点处的含沙量(全沙、分组沙均可)。假设第m汊点的所有Lout条分流河段中,第1条河段进口口门高程最低,其参考点为n1,n1点至水面的相对水深为1,此进口含沙量为Sm,1,其余分流河段中任一河段的口门高程参考点为ni,ni相对n1的相对水深为ξi,则可以导出其进口含沙量Sm,i与Sm,1的比值为Κm‚i=Sm‚iSm‚1=cc+ξii=2‚⋯‚Lout(m)(24)这样,可以得到所有的Lout条分流河段入口含沙量之间的比例符合下式:S1∶那么,根据地形确定出各分流河段的相对水深ξi之后,式(21)转化为[As]{S1}={Bs}(26)式中{S1}为各汊点的Sm,1组成的矢量,但[As]中的系数与式(22)不同,将与Km,i有关。实际上,此模式运用的关键是出流河段口门高程参考点的确定,可取进口断面的平均高程,也可取距入口一定距离之内的河床最高点。此种模式对资料要求较高,当各口门高程比较接近时,式(24)中的比值将接近于1,与简单方式是相同的。分流河段分沙比的物理机理式(24)结构比较繁杂,而且每次冲淤计算之后都要重新对其修正。河网计算涉及区域大,地形等资料收集不易,而且河网中的汊点可能连接多个分流河道,出于计算速度和精度考虑,要求提出一个能够反映物理实质,简洁而又能满足实际精度需要的分沙模式。造成各分流口门含沙量不同的主要原因是水流条件不同,而反映一定条件下水流挟沙能力的指标为S*,所以,影响分沙比的各种因素的综合作用很大程度就体现在S*的大小上。出于这种考虑,本文提出以分流河段进口断面的挟沙力S*值确定汊点的分沙比,认为存在下式的关系,这在物理机理上是合理的。S1∶S2∶⋯∶SLout=S*1∶S*2∶⋯∶S*Lout(27)根据式(27),在每个汊点m均存在一个比例常数Nm,使得Si=Νm?S*ii=1‚2‚⋯‚Ιout(l)(28)当然,在不同的时间Nm是变化的。这样,式(21)汊点方程组转化为[As]{Κ}={Bs}(29)式中矩阵[As]中元素只与流量和挟沙力有关;{Bs}为已知向量。[As]和{Bs}由式(19)或式(20)结合式(21)得到,{K}为未知向量(各个节点的Nm值)。在挟沙力公式中,含沙量与流速的高次方成正比,此种模式实际是以流速为主分配沙量,当各分流口门流速相近时,就变成了简单分沙模式。1.4点方程组求解不论采取哪种分沙模式,汊点方程组式(22)、式(26)或式(29)结构类似,与水流求解过程中的汊点方程组也类似,所以可采取相同的解法。当汊点数较多时,同样可采取分组的解法,重新利用求解水流方程过程中的分组情况,为计算程序的编制提供方便。求解汊点泥沙方程组之后,通过逐步回代,可以得到所有断面的含沙量Sn+1。2河网沉积物模型的适用性2.1非均匀沙模型悬移质泥沙方程的求解并未区分均匀沙和非均匀沙,可将其当作均匀沙,也可以是非均匀沙中的某一粒径组。如果为非均匀沙模型,对每一粒径组均可列出相应的汊点方程组,而河床的总变形将是每一粒径组所引起的变形之和。除此之外,还将涉及非均匀沙挟沙力及其级配的计算、活动层选取、床沙级配调整、床沙粗化或细化过程中糙率的调整等系列问题。2.2各口门处的截沙模式:丁君松模式以上给出的三种分沙模式可以看出,第一种简单分沙模式是丁君松模式和挟沙力模式在特殊情况下的特例,总体来讲处理比较粗糙。当对泥沙模拟的精度要求不高,或者当各口门水力条件相近,口门高程相差不大时,可以考虑采用此种模式。丁君松模式注意了含沙量沿垂线分布的不均匀和各分流口门地形(主要是高程)不同造成的含沙量不同,而挟沙力模式注重了水力条件不同造成各口门处的挟沙力区别。这两种方式分别注重了问题的两个方面。天然情况下,水流、泥沙运动相互影响,相互耦合,这是河流发生自调整作用的原因,河网也不例外。设想河网中连接某一汊点的分流河段进口口门高程、形态类似,而河段的河床比降、糙率不同,因而水力要素不同。按照推广后的丁君松模式,它们的进口含沙量相同,根据河床的自调整原理,各河段将发生冲刷或淤积来调整河床形态,此过程中同时也淤高或蚀低了口门高程,直至冲淤后口门进沙量与水流相适应。从这个角度讲,挟沙力模式与此是一致的,因为河床调整总是以加大或减小比降达到减小含沙量与挟沙力的差值为目标,而调整过程比较明显地表现在河床上段的泥沙冲淤上。只不过丁君松模式更注重了水深对泥沙粒径的分选作用,当各分流口门中存在一个河底高程明显较低的主汊时(例如长江三口),此种模式更适用,这也与它是从单一河道的分汊情况推广过来的这一背景相符合。实际应用中,还可以根据分汊的夹角等因素对丁君松模式进行修正。实际上,式(26)和式(29)形式相同,在同一模型中可以根据汊点的实际情况灵活选用两种模式的任一种,只是相关系数存在差别,不影响汊点矩阵的形成和求解。2.3水体沙量计算实际的河网计算中,经常有可蓄水汊点,如湖泊、蓄洪区等,只要提供水位与水面面积(容积)曲线,仍可用汊点水流连续方程(7)计算,但含沙量需要单独处理。进水情况,含沙量可按以上的模式计算。进入汊点的泥沙将发生沉降,沉降的时间将影响流出汊点的含沙量。设泥沙沉速为ω,汊点水位为Z,出流口门高程为Zd,hd=Z-Zd。泥沙降到口门以下需要时间Tω=hd/ω,当发生水体在汊点中贮留时间超过Tω时,可近似认为出流含沙量为0。当T<Tω时,泥沙经过时间T下降距离d=ω?T,则水面下深度为d的范围内含沙量近似为0,再以下的含沙量沿垂线分布符合张瑞谨公式(23)。从Zd向水面积分,可得出流平均含沙量。以上处理办法既适用于均匀沙也可以是非均匀沙的某一组。3蓄滞洪区一维河网模型本文采用杜家台蓄滞洪区的实测资料对所建模型进行验证和分析。杜家台分洪闸位于汉江下游仙桃市区以东,上距丹江口水库502km,下距汉江河口150km。在汉江干流仙桃以下河段超过安全泄量时(在5250～9150m3/s之间,根据汉口水位高低来确定),超额洪水经过杜家台分洪闸分泄,由分洪道引入汉南泛区调蓄,再经黄陵矶闸泄入长江。杜家台总共有11个蓄滞洪区,包括下东城垸、上东城垸、新农垸、桐湖斗昌湖、洪北东垸、洪北西垸、曲口垸、兴无垸、银莲湖、消泗外垸以及未围部