河道及河口一维及二维嵌套泥沙数学模型

2001年10月
水 利 学 报SH UI LI X UE BAO 第10期
收稿日期:2000208230
基金项目:国家自然科学基金及水利部联合资助重大项目(59890200).
作者简介:张修忠(1972-),男,山东临沂人,博士生.
文章编号:055929350(2001)1020082206河道及河口一维及二维嵌套泥沙数学模型
张修忠1,王光谦1
(11清华大学水沙科学教育部重点实验室,北京　100084)
摘要:建立了一种河道及河口一、二维嵌套的泥沙数学模型,对基本的控制方程、方程的离散和求解方法、嵌套连接条件以及非均匀沙的处理等问题进行了研究.以非恒定非均匀不平衡输沙理论作为本文建模的基础,为方便处理二维计算域的不规则边界,采用有限元数值离散格式.验证算例对河道做一维简化,对口外海域做二维处理,通过交界面的水位、流量和含沙量等的传递,在每一迭代步内进行耦合计算.数值模拟结果与实测资料吻合较好,且计算省时,表明本文建立的嵌套模型是一种解决某些实际工程问题的可靠的和高效的工具.关键词:河口;泥沙输运;嵌套连接;有限元离散
中图号:T V149 文献标识码:A
泥沙数学模型作为研究和解决河流、水库和近海等水域的水流运动和泥沙冲淤问题的有效工具,已得到了较为普遍的应用.一维模型计算省时,可快速方便地进行长河段、长时期的洪水和河床演变预报,但无法给出各物理量在平面范围的分布,因而在模拟河床细部变形、河口和港湾等水域的流动和冲淤问题时,显得无能为力.水深积分的二维模型克服了一维模型的缺陷,但因计算量剧增,模拟长河段、长系列、平面大范围的水流运动和河床演变问题时很不经济,即使是短时期问题也不易做到实时预报.因此,将一维和平面二维模型嵌套连接,发挥其各自的优势,对于解决许多生产问题是必要的和有意义的.文献[1]在这方面做了比较细致的研究工作,文献[2]应用一、二维嵌套技术成功的模拟了黄河口的演变.
1　水流泥沙数学模型及其求解方法
111　河道一维非恒定流水沙方程　河道水流运动的圣维南方程:
5A 5t +5Q 5x
=0(1)5Q 5t +55x (Q 2A )=-gA 5ζ5x -gA Q 2K 2(2)
悬移质不平衡输运方程及河床变形方程:
5(AS k )5t +5(QS k )5x
=-αωk B (S k -S 3k )(3)γ′5A sk 5t
=αωk B (S k -S 3k )(4)式中:A 、B 、Q 、
ζ分别为河道的过水面积、河宽、流量和水位;K 为流量模数,由谢才公式计算;S k 、S 3k 、
ωk 、A sk 分别为第k 粒径组泥沙的含沙量、挟沙力、沉速及冲淤面积;α为恢复饱和系数;
γ′为淤积物干容量;x 、t 为空间和时间变量.
112　口外平面二维水沙运动基本方程　对于平面大范围的自由表面流动,由于水深尺度一般远小于水面尺度,可以引入浅水假定以简化基本守恒方程.假定压力沿水深服从静压分布,对基本方程(N 2S 方程)沿水深积分,可得到如下守恒型的浅水方程:
55t (h )+55x j (q j
)=0(5)55t (q i )+55x j (u j q i )=f δij q j -gh 5ζ5x i -λq i +55x j (νt 5q i 5x j
)(6)悬移质不平衡输沙方程和海床变形方程:
55t (Ψ)+55x j (u j Ψ)=αωk s 3k -βΨ+55x j (εs 5Ψ5x j
)(7)γ′5Z b 5t =∑N s k =1
αωk (s k -s 3k )(8)式中:u j 、q j 为x j 方向的平均流速和单宽流量;f 为柯氏力系数;δ为系数矩阵,除δ12=1和δ21=-1外其余元素均为0;λ=g u j u j Π(C 2
h );C 为谢才阻力系数,可由曼宁公式计算;β=αωk Πh ;Ψ=hs k ;涡粘性系数νt 由νt =κu 3h Π6.计算,κ为卡门常数,u 3为摩阻流速;泥沙紊动扩散系数εs 假定与水流涡粘性系数相等;h 表示水深;水位函数ζ由水深和床底高程确定,即ζ(x j ,t )=h (x j ,t )+Z b (x j ,t );i ,j =1,2.
113　水流挟沙力　潮汐河口挟沙力可由下式表示
[3]:s 3=K V 2gh (9)
在风、浪和潮流联合作用下,流速应该是风、浪和潮的合成流速,即:
V =| V T + V b |+|V w |
(10)式中:V T 为潮流速度;V b 为风吹流的平均速度;V w 为波流的平均速度;V b =0102W ,W 为平均风
速;V w =012ch c Πh ,c 为波速;h c 为波高;K 为率定系数.
114　基本方程的有限元离散11411　河道单元的离散　河道单元的流动守恒方程和泥沙输运方程可写成如下统一形式的对流方程:5φ5t +5(U φ)5x
=F (11)式中:φ=[A ,Q ,AS k ]T ,F =[0,-gA (5ξ5x +Q 2K
2),-αωk B (S k -S 3k )]T ,U =Q ΠA .对流方程的有限元离散可写成:
M 5φi 5t =C ij φj +F i (12)
式中:M 表示集中质量矩阵,M =∫ΩN i N j d Ω;C ij 表示对流矩阵,C ij =-∫Ωw i
5UN j 5x
d Ω;F i 代表源项,F i =∫Ωw i F d Ω.11412　口外平面二维单元的离散　有限元在本质上属于非结构化网格离散方法,便于处理复杂边界问题.因此,本文对控制方程采用有限元法离散,方程(5)～(7)的弱解形式经分部积分后可得如下的空间半离散方程:
M 5h i 5t
=C ij h j (13)M 5<i 5t
=C ij <j +D ij <j +F 1+M ・F 2<i (14)
<=[q x ,q y ,hs ]T ,F 1=[F x ,F y ,F z ]T ,F 2=[-λ,-λ,-αωΠh ]T
式中:C
ij为对流矩阵;D ij为扩散矩阵;由下列各式表示:
C ij=-∫Ωw i5uN j5x+5νN j5y dΩ
D ij=-∫Ωνt5N i5x5N j5x+5N i5y5N j5y dΩ
F x=-∫Ωgh5ζ5x w i dΩ　F y=-∫Ωgh5ζ5y w i dΩ　F z=∫Ω(αws3)w i dΩ
式中:N、w
i分别为插值函数与权函数.
若上述离散中的权函数与插值函数相等,则构成经典的G alerkin有限元法.对于对流占优问题, G alerkin法等价于中心差分格式,因缺乏足够的耗散,往往导致数值振荡.为此,本文采用高分辨率格式对流项进行重构,即通过引入几乎相等的扩散与反扩散以保证格式的高精度,同时利用限制因子保证影响系数的非负性及解的保单调性[4].
115　离散方程的求解　为使计算收敛或加快收敛,离散中对源项M・F2<进行负坡线性比.离散后的方程(13)、(14)为常微分方程,可采用多种显式或隐式方法求解.本文对时间导数项采用C2N格式离散,对离散后的代数方程采用QMR[5]方法迭代求解,该方法具有节省内存,收敛快的优点.
本文顺序求解离散后的方程.对于二维海域,先由二维对流扩散输运方程(6)计算流速,由对流输运方程(5)计算水深,由悬沙对流扩散方程(7)计算含沙量,最后由河床变形方程(8)计算节点的冲淤深度.对于河道一维计算,运动方程的单元离散转化为求解流量Q的方程,连续方程的单元离散方程则转化为求解水位ζ的方程.
2　嵌套连接条件
一、二维嵌套模型是通过交界面连接的,由于一维模型只给出物理量的断面平均值,二维模型给出节点的水深平均值,因此在交界面上存在各物理量断面平均值和节点垂线平均值的相互转化和衔接问题.水沙运动在交界面上的连续性是模型嵌套连接的基本原则,因此,一、二维嵌套的连接条件是:
水位相等,即:ZB=∫B0z d y.式中:Z为断面平均水位;z为节点的水位;B表示交界断面的水面宽度.
流量相等,即:Q=∫B0uh d y.式中:Q为通过交界面的流量:u、h分别是节点的垂线平均流速和水深.
悬移质输沙量相等,即:SQ=∫B0uhs d y.式中:S、s表示交界面的断面平均含沙量和节点垂线平均含沙量.
此外,还有阻力、挟沙力、河床变形等连续条件.本文口门水位由二维控制,流量由河道一维计算给出,含沙量进行相互传递.
3　非均匀沙水流挟沙力级配及床沙级配计算
311　挟沙力级配计算　鉴于水流中的泥沙源于上游水流挟带和床沙紊动扩散进入,因此由水流条件和床沙组成推求非均匀沙的分组挟沙力的做法是较为合理的.本文采用李义天通过建立输沙平衡状态下的床沙质级配和床沙级配间的关系以及垂线平均悬沙浓度和河底悬沙浓度之间的关系得到的挟沙力级配公式[6].
312　床沙级配计算　床沙级配随河床冲淤而变化,对阻力、输沙率及河床冲淤影响显著.若已知各
粒径组泥沙的冲淤厚度ΔH s
k及总的冲淤厚度ΔH s,则床沙级配调整计算可分为以下两种情况.
(1)ΔH s>0,即发生淤积的情况,床沙活动层级配由下式计算:
ΔP t
+Δt bk =[ΔHs k +ΔP t bk (H t m -ΔHs )]ΠH t +Δt
m 式中:ΔP t bk 、ΔP t
+Δt bk 分别为t 时刻和t +Δt 时刻的床沙活动层级配;H t
m 、H t
+Δt m 为相应时刻的床沙活
动层的厚度.
(2)ΔH s <0,即发生冲刷的情况,床沙活动层级配由下式计算:
ΔP t
+Δt bk =[ΔHs k +ΔP t bk H t m +|ΔHs |ΔP remk ]ΠH t +Δt m
式中:ΔP t bk 、ΔP t
+Δt bk 、H t
m 、H t
+Δt m 同前,ΔP remk 为若干个记忆层内的床沙平均级配.
床沙活动层是指河床发生冲淤变化过程中,河床表层参与河床冲淤变形的那一层床沙,水流挟带
的泥沙与床沙的交换在这里发生,河床冲淤变形也在这一层里发生[7].床沙活动层厚度是指一个冲淤
计算时段内感受到水流作用并且泥沙组成发生变化的床沙厚度[7,8].受河床变形和来水来沙条件的影
响,活动层的厚度和组成不断变化.由于问题的复杂,要从数学上严格定义和表达活动层厚度目前还比较困难.尽管有许多学者对这一问题进行了研究,给出了一些计算方法.但这一问题不能考虑过细,一方面缺乏床沙级配沿垂向变化的实测资料;另一方面考虑过细不一定能提高精度.故本模型采用较常用的处理方法,即取一固定值2m.
313　床沙级配的分层记忆模式　为了模拟床沙组成的变化过程,将床沙划分为床沙活动层及其下面
的记忆层[8]两部分.记忆层可根据实际情况分n 层.计算中,当河床发生淤积时,记忆层层数增加,
增加层的级配为t 时刻的床沙活动层级配ΔP t bk .当河床发生冲刷时,根据冲刷量的大小,记忆层数
相应减少,且级配作相应的调整.
4　模型的验证
411　计算条件　漳卫新河是海河流域南系的一条尾闾河道,担负着漳河、卫河的泄洪排涝任务.自1973年扩大治理以来,由于入海径流少,辛集闸闸下河道被潮汐动力所控制,源源不断的海相来沙使河道严重淤积.据94年实测地形资料分析[9],淤积河道长达26km ,淤积总量达到1262万m 3,河道行洪能力下降47%.本文验证计算的河道一维计算域取自漳卫新河的辛集闸至河口,长3716km ;口外海域的下边界至-20m 等深线,纵向长30km ,横向宽20km.河道地形资料采用94年大断面资料,口外地形采用1∶50000海图.
412　边界条件　河道进口给定流量、含沙量过程,口外开边界条件采用潮位控制,岸边界采用水流无滑移条件.口外各角点水位由实测潮位根据潮波传播相位差推延得到,再根据域内测点流速过程验证情况稍作调整,以90年实测大潮潮型概化计算潮型.一、二维连接断面采用流速边界,并按曼宁公式进行分配.
413　有关参数的选取[10]　在现有的认识条件下,河口水沙预测的关键是选取可靠的基本参数,如糙
率、挟沙力系数和泥沙恢复饱和系数等.为此,需对河口现状水流泥沙条件进行验证,它一方面是对数学模型本身的检验,另一方面也是率定水流泥沙基本参数,为各方案科学预报提供依据.
河道糙率采用01025,河口二维海域糙率采用0102;水流挟沙系数采用海河口数据K=100;淤积
物干容重取0165t Πm 3;根据验证计算确定泥沙恢复饱和系数α冲=011,α淤=0125;波高取大口河测波
站平均波高.
图1　计算与实测潮位过程对比图2　计算与实测流速过程对比
414　验证计算结果　图1～3给出了94年8月26～27日河口处的潮位、流速和含沙量计算与实测的对比,图中零时刻对应于26日14时.可以看出计算与实测潮位、流速吻合良好,表明水流计算参数的选取是合理的,计算方法也是可靠的.含沙量过程计算与实测有一定差别,主要是由于在潮流和波浪共同作用下泥沙参数的选取还有待进一步改进.
图4给出了83年～93年河道累计淤积量计算与实测的对比,全河段累计淤积量计算值约590万m3.图5～6给出了河口局部涨急和落急流场,可以看出,涨潮流速明显大于落潮流速,与实测资料一致.这也是涨、落潮输沙不平衡,河道淤积的一个重要原因.受资料限制,口外海床变形未作验证.
图3　计算与实测含沙量过程对比图4　计算与实测河道累计淤积量对比
图5　涨急局部流速矢量场图6　落急局部流速矢量场
5　结语
对口外海域进行平面二维计算,对河道采用一维模拟;或者对河道流动复杂段应用二维模型,流动简单或顺直河段应用一维模型的一、二维耦合算法,既具有一维模型的快速方便,又能获得局部河段或平面大范围的细部信息.这样可以用较少的机时复演和预测长河段的河床变形及其重点段的细部变形,是一种解决某些实际工程问题的有效方法.
致谢:论文得到大连理工大学土木系金生教授的指导和帮助,在此表示衷心的感谢.
参　考　文　献:
[1]　Wu W M,Li Y T.One2and T w o2Dimensional nesting mathematical m odel for river flow and sedimentation[C],5th
International sym posium on river sedimentation,1992,K arlsruhe547-554.
[2]　Zhang S Q.One2D and T w o2D combined m odel for estuary sedimentation[J],Int.J.Sediment Research,1999,14
(1):37-45.
[3]　刘家驹,张镜潮.淤泥质海岸航道、港地、淤积计算方法及其应用推广[J].水利水运科学研究,1993
(4).
[4]　张修忠,王光谦,金生.浅水流动有限元分析及其高分辨率格式[J].长江科学院院报,2001(1).
[5]　Freund R W,Nachtigal N M.An im plementation of the QMR method based on coupled tw o2term recurrence[J].SI2
AM.J.Sci C om put.,1994,15(2):313-337.
[6]　李义天.冲淤平衡状态下的床沙质级配初探[J].泥沙研究,1987,(3).
[7]　李义天,胡海明.床沙混合活动层计算方法探讨[J].泥沙研究,1994,(1):64-71.
[8]　吴卫民,等.河床床沙组成数值模拟方法[J].武汉水利电力大学学报,1994,27(3):320-327.
[9]　王文治,梁永立.漳卫新河冲淤变化及发展趋势的分析[R].水利部天津勘察设计研究院,1996.
[10]　金生.漳卫新河河口泥沙冲淤计算[R].大连理工大学,2000.
12D and22D nesting sediment transport model for rivers and estuaries
ZH ANG X iu2zhong1,W ANG G uang2qian1
(11T singhua Univer sity,Beijing　100084,China)
Abstract:A12D and22D combined sediment transport m odel for rivers and estuaries is presented.The basic equation,the numerical method,the coupling conditions and the treatments for non2uniform sediment are stud2 ied.The m odel is based on unsteady non2uniform and non2equilibrium sediment transport theory.The finite el2 ement method is abopted to solve the g overning equations for its capability of accepting com plex geometry.The m odel is verified by the simulations of the flow and sediment transport in the estuaries of the Zhangweixin River,
in which the river area is treated as12D and the sea area is treated as22D.By trans ferring the water level,dis2 charge and sediment concentration at the interface,the coupling calculations are conducted in each iterative step.The results are in g ood agreement w ith the measured data and a lot of CPU time is saved,which shows that the proposed m odel is reliable and high efficiency in solving practical engineering problems.
K ey w ords:estuary;sediment transport;nesting linking;finite element discretization。