倒虹吸工程河段洪水与河床变形的数值模拟（刘培斌,孙东坡,杨明）
关键词：渠河交叉；二维水沙数值模拟；壅水；床面冲刷
1 研究问题
南水北调中线总干渠沿线与许多河流交叉，其中在河北穿越七里河的交叉建筑物型式为渠穿河倒虹吸工程。倒虹吸设计长度初选为700m，设计洪水标准为百年一遇，设计洪峰流量2410m3/s。渠穿河倒虹吸工程的修建对该处河道水流及河床演变产生的影响，倒虹吸工程的位置、尺寸、埋置深度是干渠设计所关心的重大问题。利用二维水沙数学模型，可以较好地模拟反映渠、河交叉工程附近水流与河床变形状况，主要研究：在该河修建渠穿河倒虹吸后，交叉工程附近河段流场流速、壅水及河床冲淤变形程度。通过河流模拟，分析、评价交叉工程设计方案并提出工程修改建议。
2 交叉段河道特性
在交叉工程附近七里河分为两汊，属宽浅型沙质河床：河道北槽较大，主流靠左岸。中泓处的河床质表层为粗沙，d50在1.0mm左右，滩地及两岸为中沙及壤土。
该河为冲积性河流，河床形态、演变规律与洪水造床作用和常年水沙特性密切相关。据河道查勘与河床地形图分析：1963年洪水前期引起河床较强的冲刷，洪水后河道又普遍回淤；河段滩槽高差一般为2～3m,河槽宽850～1200m。交叉工程附近河床组成沿流程分选明显，由粗至细变化较大。1980年以来河道受人类活动影响严重，由于河道内多处挖取粗沙，原主槽回淤的泥沙被大量挖走，有的挖深达3～４m。多年小水作用及人类挖沙已经使现河道形成人为不连续窄深槽，窄槽宽度一般为100～300m。这造成在现状河道小洪水易于归槽且侧侵蚀较为明显，见图1。
| 3 平面二维水沙数值模拟 3.1 数学模型的基本理论 对于宽浅型河流，水深平均的二维水沙控制方程可较好地反映河流中挟沙水流运动特征。本模型的水流基本方程由三维时均雷诺方程沿水深积分得到，并以混长紊流模型求解紊动切应力：模型采用悬移质泥沙扩散方程与河床变形方程求解河床冲淤变化，由床面冲淤临界切应力判断床面泥沙冲淤状态及床面稳定条件。该数学模型已在一些复杂工程中得到成功应用[３]，能较准确地模拟、预测一般冲积性河流上，河流工程附近的水沙运动与河床变形。 |
|
3.1.1 控制方程
| 水流连续方程 |
| (1) |
| 水流运动方程 |
| (2) |
| (3) |
| 悬移质输运扩散方程 | |
| (4) |
| 河床变形方程 | |
| (5) |
| 水流挟沙力方程 | |
S*=k(U3/gRω)m | (6) |
对于散粒沙河床，床面稳定控制的辅助方程可表示成局部区域稳定控制条件
τ*c≥τ*0或τ*c/τ*0≥1 | (7) |
| 式中τ*c=τc/(γS-γ)d=f(U*cd/v) | (8) |
τ*0=τ0/(γS-γ)d=γRJ/(γS-γ)d | (9) |
式中ξ—水位；H—水深；u,v-x,y向水深平均流速；U—合速度；U*—摩阻流速，R—水力半径；e,ν—水流涡粘系数、运动粘滞系数；β—对流项修正系数，S—含沙量，z—河床高程，J—水力坡度；D，α—泥沙扩散系数、恢复饱和系数；ω—泥沙沉速；C—Chezy系数；S*—挟沙力，k、m—挟沙力系数、指数，γS、γ′—泥沙容重与干容重；τ*0—床面无量纲水流切应力，τc—无量纲谢尔兹临界切应力；τ0—床面水流切应力，τc—谢尔兹临界切应力；d—床面分层粒径。
3.1.2 数值计算格式
(1)离散网格及变量分布 模型采用非均匀网格，可以在研究量变化梯度较大的重要局部区域设置细密网格，在研究量变化较平缓的非主要区域设置较稀疏的网格。在离散网格上，标量(ξ,H，S)被安排在单元中央，矢量(速度u,v)安排在单元的四周，交错网格上物理变量的位置相互错开。分别用四个一维数组(xu(),yv(),xh(),yh())来确定水位、各流速分量的坐标[３]。
(2)离散格式 本数学模型采用较成熟的“交替方向隐式差分逐行求解”方法，其特点是：将时间步长(TIME STEP)分成前后两个半步，在前半个时间步长取某一个方向为隐式；为保持对称，在后半个时间步长改变隐式方向。每一个时间步长，都这样交替的改变隐式方向计算。在前后两个时间半步对控制方程进行离散，为了物理概念上的清晰和格式的稳定有效，引进控制体同时在进行对流项离散时引起迎风格式。在前半个时间步长，将连续方程与Y向动量方程联立，对u,ξ进行隐式求解；在得到水流条件后隐式求解关于含沙浓度的传移输运方程。在后半个时间步长，将连续方程与Y向动量方程联立，对ν,ξ进行隐式求解；也在获得水流条件以后隐式求解泥沙传移输运方程。
3.1.3 初始条件及边界条件
(1)初始条件 对于给定的计算区域，在时间t=0时，令：
ξ｜t=0=ξ0(x,y);u｜t=0=u0(x,y);ν｜t=0=ν0(x,y) |
(2)开边界条件ξ(x,y,t)=ξopb(x,y,t)或Q=Qopb(t) 以及Sξ(x,y,t)=Sopb(x,y,t)。
其中ξopb，Qopb以及Sopb分别为开边界上已知的水位、流量以及含沙量，一般由实测水文资料确定。
对于具有宽滩深槽复杂地形的河道，模型进口开边界上宜给定洪水水位过程或流量过程，便于处理全断面流量分配。上游开边界所需相应的悬移质含沙量及过程如无法实测资料，则借用附近河流相应洪水的含沙量及过程。下游开边界处的含沙量可按第二类边界条件确定。陆地边界即河流岸边，满足固壁非穿越及无滑移条件，边界上的法向、切向流速为零。
(3)动边界处理 在宽浅河流中的边滩和江心洲随着水位的变化，其水边线也不断改变，形成所谓动边界。模型相应采用水边界全区自动跟踪的处理方法，设置一个跟踪指标数组Iwet()。先将最大可能的淹没区域包纳在计算域区，在计算过程中根据计算单元内的水深来判断该单元是淹没或是露出，即Iwet()应该赋0还是赋1(0表示露出，为陆地；1表示淹没，为应计算的水域)，凡是陆地单元均不纳入计算范围。这种动边界智能跟踪的处理方法特别适合边滩和江心洲淹没与出露频繁的水域，使流场边界条件模拟更为真实。
3.2 七里河数学模型的建立
3.2.1 数学模型范围及网格剖分
根据数值模拟的研究任务和交叉工程附近河段的河道形态，河床组成及洪痕调查情况，充分考虑上、下游边界的水流条件，确定模型上边界取在距交叉工程上游4.0km处；模型下边界取在距交叉工程下游4.5km处，沿河道模型总长约8.5km。侧边界在考虑洪水的最大可能淹没范围及地形条件，基本沿陡坎近岸高地选取。模型最大宽度为2.5km，计算区域总面积为21.25km２。
为了反映交叉工程河段河道地形，适应不同区域流场计算精度要求，交叉工程附近的重点区域，采用加密网格为50×50m;对于一般非重点计算区域，网格尺寸为100×50m。模型全部节点为4500个。
3.2.2 数学模型调试与验证
(1)调试依据与要求 模型调试主要依据是1963年8月和1996年8月该河发生的两场重要洪水的调查洪痕(以下简称“63.8”洪水、“96.8”洪水)；水文分析所得计算河段进出口附近断面的水位流量关系以及河床冲淤变形有关的查勘资料与河床探测资料。河道地形采用1966年1/10000地形图并参考1994、1996两年汛后河道的纵横断面测量资料。考虑“96.8”洪水河道侧侵蚀严重的特点，对河道地形进行了适当修正。“63.8”洪水属特大洪水，“96.8”洪水属一般大洪水，分别选用这两场洪水验证模型流场，可使模型能正确模拟不同洪水条件，不同阻力特征的流场。
(2)调试结果 二维水沙数学模型反演“63.8”、.96.8”洪水的计算成果经后处理，绘制成流场流速矢量图。流场矢量图所反映的流场的主流走向、流态基本与调查情况一致。这里仅给出“96.8”洪水流场中主要洪痕点水位调查值与计算值的对比，见图2。河床糙率上段取0.030～0.035，下段取0.025～0.032。
通过对洪水流场主流走向、流态、流势、流速分布等方面模拟与洪痕点水位校验对比，数学模型反演的“63.8”和“96.8”洪水与实际洪水调查情况基本一致，正确反映了两类洪水的洪水河势及行洪特点。数学模型在河道地形处理，糙率选择，冲淤强度控制参数方面的选择是合理的。
3.3 数学模型计算成果与分析
本次计算根据交叉工程设计方案，对设计洪水进行了数值模拟，然后根据流场状况，进一步提出修改工程方案并进行数值模拟和分析对比。
3.3.1 自然状态计算成果
在未修工程的现状条件下，对设计洪水进行了自然状态河段的数值模拟。由洪峰期计算河段平面流速分布(图3)可以看出，在设计洪水时，交叉工程断面附近水流收缩集中，河槽冲刷严重。受江心洲影响，河道分成两汊又渐次合一(北汊较大)，南岸高地下游有部分回流区。
|
|
| 3.3.2 设计方案计算成果 交叉工程设计倒虹吸段长度为700m，按设计洪水进行了设计工程方案的计算。为了反映洪水演进过程中有关特征值的变化，特别是为确定最高洪水位，最大冲刷深度等有关成果，根据设计洪水过程的特点，设定洪水计算过程中相应的计算成果输出时刻，一般选取23～29个时刻，包含洪水演进的各特征时期。这里仅给出洪峰期流场流速矢量分布(图4)以及交叉工程断面处流速及冲淤分布(图5)。 从计算结果看，总干渠交叉工程修建后，在倒虹吸工程河段的主流分为两股，河道北侧深槽略大一些。两股深槽处单宽流量和流速均较大，最大流速达3.8m/s。冲刷也主要集中在两股深槽，最大冲刷在北侧深槽，最大冲刷深度约为1.94m。从冲刷发展过程看，随洪水发展两股深槽的冲刷也不断发展变化，南槽持续冲刷而北槽则在洪峰过后又略有回淤。 |
|
3.3.3 修改方案计算成果
根据自然状态和设计工程状态下流场流态、壅水、冲刷深度及部位的计算成果分析，得知南岸滩地水流流速不高，近期调查南槽仅略有增大趋势，在此适当约束水流不会产生很大影响。因此可以考虑进行修改、优化工程布置，这里提出了分别在南、北岸倒虹吸工程进口、出处将倒虹吸段缩窄50m、和将倒虹吸工程全部安排在北槽的三种修改方案，相应倒虹吸工程段长度均为650m。对不同修改方案在设计洪水条件下进行了数值模拟计算，不同方案时在交叉工程附近引起的壅水程度、冲刷程度及发展规律各不相同。交叉工程断面处的水位、流速、河床变形可见表1、图6。
表1 不同修改方案交叉工程断面特征值
Characteristic values of cross-section in different modified alternatives
| 修改方案 | 倒虹吸长度(m) | 水位(m) | 平均冲深(m) | 最大壅水高度(m) | 最大流速(m/s) |
| 1 | 650 | 85.24 | 1.37 | 0.40 | 3.40 |
| 2 | 650 | 85.42 | 1.74 | 0.55 | 3.56 |
| 3 | 650 | 85.93 | 0.96 | 0.66 | 2.56 |
| 注：修改方案1是在原设计基础上在南侧进口缩进50m；修改方案2是在原设计基础上在北侧进口缩进50m；修改方案3是将倒虹吸段650m全部放在北槽。 |
|
|
图5 交叉工程断面流速与冲淤分布(设计洪水) | 图6 不同方案下交叉工程断面河床对比(设计洪水) |
Velocity and erosion profile in cross section (designed flood) | Comparison of riverbed of cross section in different alternatives |
3.3.4 修改方案成果分析
由于受河道地形的影响，在设计方案和修改方案中，洪水主流略偏靠北侧深槽。受交叉建筑物挤压，南北两汊流速均较大，北汊最大单宽流量可达15m2/s以上，应注意边岸工程的防护。修改方案中压缩段基本属于近滩范围，故这对倒虹吸工程处河道实际主流过水影响不大。交叉工程断面附近壅水最大高度(修改方案1)为0.40m，壅水最大范围约1000m。由于流速较大，所以交叉工程断面在不同流量级洪水时，均有主槽全断面均有冲刷，而两侧滩地则有冲有淤。北汊河槽主要是河床下切，落峰时略有回淤；南槽除河床下切外还有侧向侵蚀，河床被不断冲深、拓宽。洪水过程中河床一直冲刷，但洪峰后期冲刷速度转缓。
由于河床变形剧烈，洪水演进中一般Zmax出现较早，而此时流量Q未达到最大值，冲刷强度还继续增加，而后次第出现Qmax、νmax，此时冲刷强度达到最大，此后冲刷虽继续进行但发展速率逐渐减小。在洪水后期，全断面便开始由冲转变为微淤。
比较不同方案的计算结果，可以看出修改方案1既可改善北侧水流条件又适当考虑保留南槽过流能力，能更好适应上游河势改变可能引起的主流调整。
4 结语
1.采用平面二维水沙数学模型研究七里河交叉工程河段洪水及河床变形问题是可以满足工程要求的。
2.通过历史洪水的反演与洪痕校核，进行了模型参数的调试和糙率调整，计算表明这对于正确模拟和预测不同条件洪水与河床变形是必须的。
3.根据设计洪水条件，对不同工程方案进行了模拟计算，所提出的有关流场流态、主流走向、壅水、冲刷状况，交叉工程断面处壅水冲刷的有关特征值可做为交叉建筑物初步设计的依据。
4.根据不同工程方案的流场流态、流速及工程断面处的壅水、冲刷等状况的分析，建议采用修改方案1(北侧倒虹吸缩短50m)并适当布设导流堤等边岸防冲措施。最佳方案的选定还须进行全面的技术经济比较。
参 考 文 献
[1] 曹祖德，王运洪。水动力泥沙数值模拟。天津大学出版社，1994.[2] T.J.Chung. Finite Element Analysis G-H International Book Company,1978.
[3] 孙东坡等。河流模拟理论与实践。武汉水利电力大学出版社，1998，40～49.
