1、洞庭湖地区水系水动力耦合数值模型第 39 卷第 1 期2008 年 1 月海洋与湖沼OCEANOLOGIAETLIMNOLOGIASINICAV_o1-39.No.1洞庭湖地区水系水动力耦合数值模型赖锡军姜加虎黄群(中国科学院南京地理与湖泊研究所南京 210008)提要针对洞庭湖地区复杂水系和复杂水情下的实际水流运动特征,建立了完全基于水力学方法的洞庭湖地区水系的一,二维耦合的全局水动力学模型,实现了不同水体数值模拟的自动有机衔接.模型包括湖泊,河网水系两部分.对长江干流荆江段,湘资沅澧四水,三口分流河道,运用能够反映河道主槽和边滩不同行洪特征的扩展一维水量模型进行模拟;对于洞庭湖湖区部分,采
2、用二维非结构的有限体积法建立水动力模型;应用重叠投影法实现模型的耦合 .模型具有动边界的自动处理,河道分区计算,分洪,溃堤过程实时模拟等功能.剖分的网格较为细致地刻画了湖区的地形,使模型能够模拟“高水湖相 ,低水河相 “的湖流特征.运用 1996 年 7 月的实际洪水过程,进行了模型的验证,较为准确地模拟了实际水流特征.关键词洞庭湖,一维,二维,水动力学模型中图分类号 P731洞庭湖地区水系洪水受长江,湘资沅澧四水来流量和长江下泄量的共同影响.该区域内的洪水研究必须综合考虑河网与湖泊之间的相互影响,相互作用问题.开发涵盖长江干流荆江段,湘资沅澧四水,三口分流河道以及洞庭湖湖区的大系统全局水动力
3、学模型,实现各种数值模拟的自动有机衔接,是研究洞庭湖地区洪水运动特征的基础.长江中下游水动力学模拟模型常见的有:一维和一二维连接的水动力学模型(龙超平等 ,2002;仲志余等,1999;吴作平等,2003;胡四一等,1998,2002;郭熙灵等,2002;谭维炎等,1996).其中以谭维炎,胡四一等人(谭维炎等,1996;胡四一等,1998)于 20 世纪 90 年代中期构建的长江中下游防洪系统数学模型较为突出.它利用一维,二维非恒定水动力学模型和水文学方法相结合的方法设计了模拟模型.之后,胡四一等(2002)又对原模型中河网算法,内外边界条件的处理,水流阻力计算等问题进行了改进和完善.该模型
4、为长江中下游复杂水情计算提供了一套非常有价值的思路.但是仍存在一些问题:(1) 该模型加入的一些水文经验模型修正计算技巧忽视了水动力学特性本身的相互作用机制.(2)一维河网计算仍显粗糙,不能体现平原区河道低水时主槽过流,高水边滩主槽共同行洪的典型的过流特征.(3) 内部边界处理和区块连接不够完善,例如河湖耦合计算采用直接显式连接等.(4)湖区二维模型网格设计简单,不能反映湖区水动力特征.为了进一步提高洞庭湖地区复杂水系和复杂水情下的实际水流运动特征的模拟能力和实用性能,减少人工干预,自动模拟各区水体之间的水动力交互作用.本文中作者建立了完全基于水力学方法的洞庭湖地区水系一,二维耦合的全局水动力
5、学模型.该模型包括湖泊,河网水系两部分,对长江干流荆江段,湘资沅澧四水,三口分流河道的水流运用能够反映河道主槽和边滩不同行洪特征的扩展一维水力模型进行模拟,采用 Preissman 四点隐式差分格式进行离散;对于洞庭湖湖区部分,采用非结构的有限体积法,以 Roe 平均的通量向量差格式(FDS)建立二维水力模型,模型采用 Roger 提出的方程组形式 ,成功克服了静水计算不和谐问题,同时实现了露滩和漫滩中国科学院知识创新工程方向性项目,KXCX3 一 SW 一 331 号和中国科学院南京地理与湖泊研究所知识创新工程特别支持重大项目“我国不同气候区湖泊演变与重大环境问题,发生机制及其控制对策研究“
6、,2003-2006 年.通讯作者:赖锡军,助理研究员,博士,Email:xjlainiglas.ac.crl收稿日期:2006.0120,收修改稿日期:200701261 期赖锡军等:洞庭湖地区水系水动力耦合数值模型 75的高效处理;对于整个水域,利用重叠一投影法在河网与湖泊水力模型的连接处动态耦合一维,二维模型,实现一,二维数值模型的自动准确衔接.本文中将重点论述模型的构建,并展示洞庭湖洪水计算的部分成果.1 基本方程1.1 水深平均的二维流动的基本方程在静水压强假定下,沿水深积分三维NavierStokes 方程可得水深平均的二维流动方程组.为解决非平底坡引起的静水计算不和谐问题,采用以
7、下守恒形式的方程组(Rogersetal,2003):+V.,():S(1)S:=r,1式中,=IhuI,:hu“+要(+2)huvhvhuvV+要(+2)Og4SoghS 詹+Vbl+Cw(02sino+fvhP一y-ghSfy+V2Ehv+cwPa(02coso,-扣 h其中,为向通量向量;Fy 为 Y 向通量向量;S 为源项向量;为自由面高度;h 为水深;h=h为基准面至床底的深度;12,v 分别为 X 和 Y 向水深平均的流速分量;g 为重力加速度;So=_dzb,为向的水底底坡;So:_Ozb,为 Y 向的水底底坡;oyr_=_S 肛:pnZl万U2一+Vz,为向的摩阻底坡为糙率);
8、.厂:,为 Y 向的摩阻底坡.c 为风的阻力系数;Pa 为空气的密度:为风速;为风速与 Y 轴的夹角.该形式以自由面高度代替水深 h 作为困变量,能解决高程急剧变化带来的计算失稳.1.2 一维非恒定流方程洪水泛滥的中下游平原地区的较大的江河水系,由于江河长期演化形成了复杂多变的河道过流断面.这些天然河道在长期的演化中形成了在低水位时主槽过流,高水位时滩地和主槽共同行洪的典型过流特征.为了能够真实地模拟预测这些河道的水流,必须扩展当前模型对复式河道非恒定流的模拟能力.将河道水流根据其过水特征在空间上划分为几个过流通道(如左边滩+主河槽+ 右边滩),各过流通道水流满足一维近似.假定全断面上的水位相
9、等,即断面上水流横向比降可以忽略不计;河段的纵向坡度很小(图 1).在这些假定基础上,对二维水深平均的水流方程沿横断面积分导出各过流通道沿横向(断面方向 1 平均的一维非恒定流动的扩展方程.连续性方程:+:0(2),.,jK动量方程:一I-妾(譬+Mk=-要一 c 北+(3)其中,下标 k 表示第 k 条过流通道;A 为过水断面面积;Q 为流量;Z 为水位;q 为流道和外界的质量交换;M为流道和外界的动量交换量;动量校正系数;S,为阻力坡降;S 能量耗散项;F 为流道的侧向剪切力;g 为重力加速度.然后,沿断面纵向积分得到各流道的河段积分方程.同时,在各流道在断面首尾端的水位相等的假设基础上,
10、忽略紊动引起的能量耗散项,经化简,导出如下的某河段全断面上水流的积分控制方程(赖锡军等,2005):a一十一 QB:Q,.(4)3tAfuD一gA(5)其中,上标“一 “表示河段 XkBXkE 上的均值;下标 k表示第 k 条流道;下标 B,E 分别表示河段的首,末断面;V 为河段上流道的蓄水容量:为流道的流量76 海洋与湖沼 39 卷分布系数,=Ok/O;K 为流量模数 ;全断面的总流量模数 K=Kk;等效河段长 =c首断向图 1 河段的过流通道Fig.1Explanatorysketchshowingaflowpathofareach2 数值求解首先分别对一,二维水流基本方程进行数值离散,
11、建立河网一维水量模型和二维水动力模型;然后利用模型耦合技术实现一,二维模型的准确衔接.2.1 河网水流数值模型2.1.1 方程离散采用四点加权 Preissmann 隐格式离散,并经线性化处理,可得到离散化的代数方程组(汪德燧 ,1989).2.1.2 节点方程节点(河道交汇区)处通常满足流量的连续性条件和能量守恒条件.为简化计,常以水位处处相等来近似能量守恒条件.则有:A7一 AJ(6)Z,=ZN,J=1,2,(7)式中,Z 为汇合处各河道的水位;Z为汇合处水位,又称节点水位:O,j 为各河道进入(取负号)或流出(取正号)节点的流量;A 为节点水面面积._/-IL_N 为节点水位变化率.当节
12、点面积很小,可忽略不计.若节点面积 A 较大,则需考虑节点调蓄能力.2.1.3 河网算法根据节点水位法求解河网的算法形成节点水位方程组,应用预条件最小化残差算法(PGMRES)(陈杨等,2003)求得各节点水位:最后回代求得各河道断面水位和流量.2.2 二维水流数值模型二维方程应用非结构网格有限体积法求解.界面通量采用基于 Roe 平均的黎曼近似解评估.先对方程在任意控制体积其边界为 A)作体积分,利用高斯定理将体积分化成面积分,有:UdV=-LFUndA+SdV(8)式中,n 为单元边界 A 的外法向单位向量,()?n为法向数值通量.对控制体单元取平均后,可得到有限体积法的空间离散化方程为d
13、U=一m,Aj+(9)dt一其中为单元体积,为单元面总数,A 为单元面的面积,为单元的源项平均值,单元面.,的法向对流通量为.=Fj(U)xn.对流通量,=【“,huu+寺 g(+),hvu+1g(+)n,(10)利用 Roe 型通量差分裂格式(FDS)(Roe,1981),则对流通量可以表示成:F=0.5(F(U)+F(U 一)-NIAIL(u 一 U 一)(11)式中,上标一,十分别表示界面的左右状态.人为Jacobi 矩阵 J:的特征矩阵.L,R 为左右特征向0U量,L=R .Al=$ltl+vnycl000j“n+1,ny010R=l“一 Cn00“n+vnv+cnv“+衄(12)式(
14、12)qb 的变量 “,v,C 按 Roe 状态平均方法计算,定义如下:+u-v+/=产=?:= 一,1,=产=1=一,+一+一1 期赖锡军等:洞庭湖地区水系水动力耦合数值模型 77c=水动力学模型2.3 模型耦合技术对于在大系统的非恒定水动力模拟,不同模型接口的准确衔接是计算成功的关键要素之一.作者利用重叠一投影法(赖锡军等,2002),在一,二维水域公共边界点的连接处建立了一维,二维耦合的整体模型.设定一,二维计算的耦合区,在其内进行一,二维精确解的投影,准确实现模型的连接(图 2).可以选用不同的方式实现,如水位耦合和流量耦合.但是在实际计算中,一,二维连接点处的水力要素和几何要素的对接
15、必须协调一致.2 一二维模型耦合Fig.2The1-Dand2-Dmodelcoupling2.4 动边界处理洞庭湖洪枯水情悬殊,水位变幅大,洲滩漫露过程交替进行.动态处理这类运动边界,就成为了洞庭湖地区水动力学模型的关键点之一.根据有限体积法特点,以水深为判别标准设定单元界面的类型,并运用相应的方法计算跨界面的法向通量,以保证水量平衡.在每一时间步判断单元界面两侧单元的水深h,h 月,当 hL 和 h 月均为 0,则通量为零;当 hL 和 h 月有一为零,则结合两侧单元的底高程,交界面类型可能为固壁边界,跌水,漫流等形态的边界,根据不同类型可选择瞬时溃坝解析解,堰流公式等方式估算通量;当 hL 和 h 均不为零时,则按正常的方法估算通量.开发模型时,设定某一水深限值占(本文中取为0.01m),以此来判断单元的干湿情况.如果相邻单元的水深均小于,则认为该单元为干单元,不进行计算.
