张帮稳,吴保生,章若茵

(清华大学水沙科学与水利水电工程国家重点实验室,北京 100084)

三峡水库泥沙淤积问题不仅关系到水库的使用寿命和功能发挥,而且对水库及下游的生态环境具有重要的影响[1]。近年来长江上游梯级水库的建设导致三峡入库径流量和泥沙量发生变化,严重影响了三峡水库综合效益的发挥,为此,三峡水库在“蓄清排浑”运用方式的基础上,开展了中小洪水调度和汛末提前蓄水调度,目的是提高发电和航运效益,减轻坝下游的防洪压力[2]。但汛期中小洪水调度和汛末提前蓄水会抬高库区平均水位,降低库区流速,加重库区泥沙淤积。为缓解中小洪水调度及汛末提前蓄水调度带来的泥沙淤积影响,三峡水库于2012年和2013年汛期利用洪峰和沙峰的异步运动特性开展了沙峰调度试验,取得了较好的排沙效果[3- 4]。

研究者基于实测资料对三峡库区洪峰和沙峰的运动特性进行了分析,发现三峡水库蓄水后,库水位抬高,库区水流流速减慢,沙峰输移时间较蓄水前大幅增加,沙峰多滞后于洪峰,且沙峰滞后于洪峰的时间沿程逐渐增加,有利于进行水库的沙峰调度[5- 7]。但目前关于洪峰和沙峰异步运动特性多以实测资料进行分析,基于试验研究及数值模拟开展相对较少,对洪峰和沙峰异步运动特性缺乏较为深入的认识。利用数值模拟方法对不同情况下洪峰和沙峰运动特性进行系统分析,可为应对未来愈加复杂的水沙变化和地形变化情势下水库的沙峰排沙调度提供参考。水沙数学模型在研究水沙输移、泥沙淤积及河床演变规律中有着重要的作用。基于断面平均或垂向平均的一维和二维数值模型[8- 9]能够对工程水沙问题进行一定程度的研究,但是泥沙淤积、冲刷和输移是一个非常复杂的水动力现象,特别是在复杂的自然地形,例如河漫滩、弯曲河段和沙波地形等产生二次流和流体分离等三维水动力特性,严重影响了数值模拟结果的准确性。随着计算机科学技术和高性能计算机的发展,三维模型[10- 11]逐渐应用到工程实际当中,目前大多数采用三维数值模拟方法针对局部河段或库区的泥沙输运和淤积进行研究,但对洪水传播过程中洪峰和沙峰异步运动特性研究较少。主要原因是洪峰和沙峰的运动特性与研究河流或水库的尺度(时间尺度和空间尺度)有关,较短的距离和较短的时间难以显示出洪峰和沙峰运动规律。三峡库区地形较为复杂,总长超660 km,对此进行长距离洪水演进的三维数值模拟具有较大的难度。

掌握洪峰与沙峰异步运动特性、开展入库泥沙实时监测与预报,是进行沙峰排沙调度的基础。为进一步提高沙峰预报精度,掌握更精准的沙峰调度时机、制订更科学合理的沙峰排沙调度方案,本文采用三维水沙数值模型SCHISM(Semi- implicit Cross- scale Hydroscience Integrated System Model)对三峡库区汛期洪水传播过程中洪峰和沙峰运动进行模拟,初步分析三峡水库不同蓄水位下洪峰和沙峰在万县—三峡大坝库区的异步运动特性,可为三峡水库的沙峰调度提供参考依据。

1 数值模拟方法

1.1 水流的基本方程

SCHISM水动力学模型[12- 13]的控制方程采用基于Reynolds时均N- S方程,满足静压假定和Boussinesq涡黏性假定。在笛卡尔坐标系下,对于不可压缩流体N- S方程的连续性方程为

(1)

动量守恒方程:

(2)

(3)

式中:x和y分别表示笛卡尔水平坐标,z为垂向坐标,向上为正；u,v,w分别表示3个方向的流速；t为时间；f为柯氏力系数；η为自由水面；zb为河床底高程；ρ0和ρ分别表示参考密度和混合流体的密度；g为重力加速度；Kmh和Kmv分别为水平与垂直涡黏性系数,其中垂向涡黏性系数根据紊流模型进行封闭,水平涡黏性系数采用常数化处理；pa为自由水面大气压强。

自由水面采用水位函数法处理,对连续方程(1)沿水深方向积分,可得自由水面方程

(4)

模型在河床底部的动力学边界条件由底床摩擦剪应力和底层水体的雷诺应力平衡给出

(5)

式中:tbx和tby分别为床面的x、y方向摩擦剪应力。

对于垂向紊动涡黏性系数Kmv,利用紊流闭合模型GLS(Generic Length Scale)[14]进行求解,包括紊动能k方程和通用紊动长度ψ方程。GLS紊流模型的控制方程为:

(6)

(7)

(8)

(9)

式中:k为紊动能；ψ为通用紊动长度参数；μ为盐度、温度等物质的垂向扩散系数；υk和υy分别表示紊动能和通用紊动长度的垂向扩散系数；ε为紊动耗散项；Fw为壁函数,在k- ε模式中为1；σk、σψ、cψ1、cψ2和cψ3为模型系数；M2和N2分别表示由于剪切变形和密度分层而引起的紊动能产生项,通用紊动长度ψ和紊动耗散项ε作为紊流模型中的关键参量,其表达式如下:

ψ=(cμ)mknlp,ε=(cμ)mknlp

(10)

式中:l为紊动掺混长度；cμ为常数0.3。在GLS模型k-ε模式双方程紊流模式的参数取值见表1。

表1 GLS模型k- ε模式双方程紊流参数值Table 1Parameter values of the turbulence in k- ε equation based on GLS model

1.2 悬移质泥沙输移基本方程

在悬移质泥沙动力学中,对流扩散理论将泥沙视为单一的连续介质,并假设悬移质泥沙运动与水流运动在垂直方向上存在速度差,且等于泥沙颗粒的沉降速度。基于对流- 扩散理论的三维悬移质泥沙输运方程表达式为

(11)

式中:Ci为i组含沙量；ws,i为i组泥沙颗粒的沉降速度；Ksv为泥沙垂向扩散系数,通常假定与水流紊动黏性系数呈倍数关系,可通过紊流模型求解或采用经验关系估计，Ksv=Kmv/σs,σs为Schmidt数,通常取值在0.6～1.2之间；Ksh为泥沙水平扩散系数,考虑到水平扩散的量级远小于垂向扩散,通常忽略不计；Dq为泥沙的沉积通量;Eq为泥沙的冲刷通量。三峡水库中细颗粒泥沙存在絮凝现象,在泥沙运动力学中,采用絮凝因数F反映絮凝对泥沙颗粒的沉速的影响[15],其表达式为

F=ωs,i/ω0

(12)

式中:ω0为泥沙颗粒静水中的沉速。广泛采用的絮凝因数F=αdβ,其中,α取值为0.013[16]或0.005 5[17],β取值为-1.9[15]或-1.02[16]。张地继等[16]采用α=0.005 5,β=-1.02研究了万县—庙河库区沙峰的衰减规律,本研究中采用相同的絮凝因数系数。

1.3 河床冲刷演变基本方程

根据泥沙质量守恒方程,悬移质泥沙输移过程中可将床面边界条件视为床面附近泥沙通量的处理,包括床面泥沙的沉积通量Db和冲刷通量Eb,其表达式分别为:

Db=ωs,ic1,i

(13)

Eb=E0,i(1-qi)(τf/τcr,i-1)τf>τcr,i

(14)

式中:c1,i为数值模拟中最底部一层网格的i组泥沙浓度；E0,i为经验冲刷率系数,取决于局部床面泥沙颗粒条件,取值大小范围为10-4～10-2kg/(m2·s-1)；qi为床面层i组泥沙体积分数；τcr,i为i组泥沙颗粒临界起动剪切应力；τf为水流底部的剪切应力。河床泥沙的冲淤和悬移质泥沙的净输移导致河床表面的地形变化,其公式为

(15)

式中:Δh为床面高程变化量;ρs为泥沙颗粒的密度。

2 数值模型的建立及验证

为了研究三峡库区汛期洪水传播过程中洪峰和沙峰异步运动特性,本文选择万县站—三峡大坝的库区段作为研究对象(图1)。若把寸滩或清溪场作为入口边界条件,则河段过长,计算耗时过大,计算效率极低。万县至大坝长280 km,不仅长度适中,并且也是目前少有的长距离、大范围实际水库的三维数值模拟,难度较高；另一方面原因,万县—三峡大坝库区之间仅有靠近大坝庙河站可以测量流量和含沙量,长度不足以模拟洪峰和沙峰的传播特性,因此只有万县站满足本文的研究要求,可作为入口条件。模拟中采用的资料分别来源于万县站、奉节站、巫山站、庙河站和茅坪站,如图1(a)所示,其中水文站有万县站和庙河站,具有日平均水位、流量和含沙量信息,水位站有奉节站、巫山站和茅坪站。

图1 研究河段的地形高程、水文站和网格设置Fig.1 Topographic elevation,hydrological station and grid setting of the study reach

2.1 模型网格设置

万县—三峡大坝的地形较为复杂,宽阔和狭窄的河段交替变化,水平方向上进行网格设置时主要混合使用2种网格类型,在宽阔河段混合采用四边形和三角形网格,主河道采用四边形网格,岸滩采用三角形网格,在狭窄河段采用三角形网格,网格尺度约为20～23 m,共得到495 662个网格节点和856 056个网格单元。计算初始地形条件由2011年实测地形插值得到,插值后的局部地形如图1(b)所示。此外,为了更好地模拟河道至坝前水深变化幅度大的特点并且提高计算效率,垂向上采用分层LSC2坐标[17],最大水深处可达36层,最浅处为16层,平均约为19层,每层厚度约5～6 m。图1(c)为局部河段横断面的垂向网格示意图。为了数值模拟的稳定性,时间步长设置30 s。整个计算模型在清华大学高性能计算集群上采用280个核进行计算,模拟时间为8 d。

2.2 模型边界条件

本文以三峡库区2013年汛期7月1日—8月1日作为模拟时间段,万县站的流量和含沙量过程作为模型的入口边界条件,茅坪站的水位过程作为出口边界条件(图2)。三峡库区泥沙主要以悬移质输移为主,万县站7月实测的泥沙粒径分别为0.002 mm、0.004 mm、0.008 mm、0.016 mm、0.032 mm和0.064 mm,对应的泥沙级配分别为10.3%、12.3%、20.9%、26.0%、18.5%和12.0%,考虑到泥沙0.002 mm和0.004 mm 的粒径很小,统一归为0.003 mm的泥沙进行计算,床沙依据大断面实测的泥沙级配进行插值设置,并且选择与悬移质同样的5组粒径泥沙。数值模拟泥沙输移过程中暂时没有考虑温度和盐度的影响。

图2 数值模型边界条件Fig.2 Boundary conditions in numerical model

2.3 数值模型验证

为验证数值模型模拟洪水传播和泥沙输移过程的正确性,分别对比分析了实测与模拟的奉节站和庙河站的水位,以及庙河站的流量、平均流速和含沙量。

图3对比了奉节站和庙河站水位的模拟值和实测值,可以看出数值模拟的结果和实测值吻合较好,能够较好地反映洪水传播过程中水位的变化。图4对比了庙河站断面平均的流量、流速和含沙量的模拟值和实测值,可以看出洪水流量模拟值的变化和实测值的变化过程基本一致,局部的差别可能是由于万县—三峡大坝之间支流入汇和复杂的局部地形糙率的影响,经过分析2013年7月支流平均流量的总和占万县站流量的1%,本文中暂时没有考虑支流径流对主河道洪水传播的影响。平均含沙量的模拟值与庙河站实测时刻的含沙量点较为接近,和2013年水文年鉴中日均含沙量有局部的差别。这是因为水文年鉴中日均含沙量是基于实测某一时刻的含沙量回归插值得到的[18],本身存在一定的误差,也可能是限于万县—三峡大坝之间实测资料的限制,模拟区域局部的冲淤不能很好地反映出来。

图3 水位模拟值和实测值的对比结果Fig.3 Comparison of water level between simulated and measured results

图4 庙河站平均流量、流速和平均含沙量的模拟值和实测值对比Fig.4 Comparison of average flow discharge,velocity and sediment concentration between simulated and measured results at Miaohe station

图5给出了庙河站断面流速和含沙量分布的瞬时模拟结果。从图中可以看出流速和含沙量分布与断面的几何形态有关,断面含沙量的分布不均匀,呈现出分层的现象。图6给出了庙河站垂向流速和含沙量的测量值和模拟值的对比结果,测点的位置见图5,测点之间间隔120 m。从图中可以看出垂向流速的数值模拟结果和测量值吻合较好,垂向含沙量的模拟结果和测量值存在一定的误差,数值模拟断面两侧的测点底部含沙量偏小,可能与数值模型没有考虑推移质的泥沙运动和断面几何形态有关。三维水沙数值模型能够进行断面垂向流速和含沙量的分析,相对于一维和二维数值模拟的精度更高。

图5 庙河站断面流速和含沙量模拟结果的瞬时分布Fig.5 Instantaneous distribution of simulated flow velocity and sediment concentration at Miaohe station

图6 庙河站断面垂向流速和含沙量模拟值和实测值的对比结果Fig.6 Comparison of vertical flow velocity and sediment concentration between simulated and measured results at Miaohe station

3 三峡库区坝前不同蓄水位下洪峰和沙峰异步运动特性

三峡工程于2003年6月进入围堰蓄水期,坝前水位汛期按135 m、枯季139 m运行；2006年汛后初期蓄水后,坝前水位按汛期144 m、枯季156 m运行；自2008年汛末三峡水库进行175 m试验性蓄水以来,工程进入175 m试验性蓄水期。为了更好地分析三峡水库蓄水以来汛期不同蓄水位下洪峰和沙峰异步运动特性的规律,本文设置3个研究方案。方案2为2013年7月实测洪水期间坝前水位；方案1和方案3分别为2013年7月实测洪水期间坝前水位减去和加上10 m。图7给出了3种研究方案的坝前水位变化。基于数值模拟的结果分别提取万县—三峡大坝之间重要水文站的流量和含沙量随时间的变化,用以分析洪峰和沙峰沿程的异步运动特性。

图7 不同方案的坝前水位变化Fig.7 Temporal variation in water level in front of the dam in different schemes

图8给出了整个研究区域三维流场瞬时的模拟结果,从图中可以看出上游库区的水流流速较大,坝前段水流流速较小,主要由于坝前段水深较大。

图8 研究区域三维流场瞬时的模拟结果Fig.8 Instantaneous results of the three- dimensional simulated flow field in the study zone

图9分别给出了不同方案下重要水文站的流量和含沙量随时间变化的模拟结果对比。从沿程各水文站流量的模拟结果可以看出,蓄水位的变化对洪峰传播时间的影响较小；从含沙量模拟结果可以看出,蓄水位的变化对沙峰传播时间和大小的影响较大,随着蓄水位增加,沙峰滞后洪峰的时间逐渐增加,其实质是洪水在河道型水库向下游传播过程中,随着水深的增加,流速降低导致沙峰传播越来越慢,沙峰滞后洪峰的时间也越来越大；同时由于水流流速减小,水流的挟沙能力降低,泥沙沿程不断地落淤导致造成沙峰坦化,沙峰的峰值逐渐减小。从奉节站、巫山站和巴东站的含沙量随时间变化的曲线形态可以看出在巫山站—巴东站库区间存在局部的泥沙冲刷。

图9 3种不同方案下沿程各水文站流量和含沙量的模拟结果对比Fig.9 Temporal variation of flow discharge and sediment concentration at each hydrographic station under different schemes

为了定量地反映坝前蓄水位的变化对三峡库区洪峰和沙峰异步运动特性的影响,表2分别给出了不同方案下沿程各水文站的洪峰和沙峰到达时间及沙峰滞后洪峰的时间,表明坝前蓄水位的变化对洪峰传播时间的影响较小,对沙峰传播时间的影响较大；并且沿程各水文站沙峰滞后于洪峰的时间随着坝前蓄水位增加越来越大,即沙峰滞后洪峰的时间随着水流流速的降低越来越大。方案2和方案3相对于方案1坝前沙峰滞后于洪峰的时间分别增加了6.4%和16%。

表2 不同方案下沿程各水文站洪峰和沙峰到达时间及滞后时间Table 2Arrival time and lag time of flood peak and sediment peak reaching at each hydrographic station under the different schemes

4 结 论

本文采用三维水沙数值模型SCHISM对万县—三峡大坝280 km的库区进行了洪水传播和泥沙输移的大尺度数值模拟,主要对比分析了不同坝前蓄水位下沿程各水文站洪峰和沙峰异步运动特性,结果表明:

(1) 三维数值模型SCHISM能够较好地模拟三峡库区长距离的洪水传播和泥沙输移过程,与实测值验证结果较好。

(2) 三峡水库坝前蓄水位的变化对洪峰传播时间的影响不明显,对沙峰传播时间的影响较为显著；坝前水位的增加导致水流流速减小,沙峰传播减慢,引起库区主要水文站沙峰滞后于洪峰的时间越来越大。

(3) 在万县—三峡大坝库区洪水传播过程中,随着水深的增加,水流流速减小,水流挟沙能力降低,泥沙不断的落淤导致沙峰峰值减小。

致谢:本研究得到了长江水利委员会水文局许全喜教高、武汉大学张为副教授、郑珊副教授的帮助和支持,数值模型在清华大学探索100集群上计算,特此致谢！