淹没堆积体附近水流特性研究

2023-02-21李馨楠赵红兮

人民黄河 2023年2期

李馨楠，向波，张婧，王东，赵红兮

（1.西华大学流体及动力机械教育部重点实验室，四川成都 610039；2.西华大学流体机械及工程四川省重点实验室，四川成都 610039；3.四川省公路规划勘察设计研究院有限公司，四川成都 610041）

有关研究团队已经发表一系列非淹没条件下堆积体的研究成果［3－5］，但是对淹没堆积体的相关水沙研究尚属空白。绕流结构体附近水流紊动比较强烈，尤其是回流区旋涡、气泡较多，给试验的精密测量带来困难，因此许多学者采用数值模拟技术对淹没丁坝的水流特性展开研究。堆积体和丁坝均属于绕流结构体中的一种，只是形态有所不同。与丁坝相似，堆积体存在于河道中，会影响其附近的水流特性，堆积体的存在会使上游的壅水现象更加明显，沿程水深分布将受到流量、堆积体形态、底坡坡度等因素的影响，呈现较为明显的规律。蔡喆伟等［6］为研究新型结构淹没丁坝的水流特性，基于有限体积法与自由液面捕捉法，通过流体计算软件Fluent建立三维数学模型对不同流速、水深与坝长条件下新型结构齿形丁坝的水流特性进行研究。宁健等［7］通过建立数学模型，模拟在水流作用下的泥沙运动，预测了泥沙侵蚀、对流扩散以及沉积等过程，同时考虑河床坡度对泥沙输移以及临界希尔兹数的影响，得到丁坝附近流场、床面变形以及冲刷坑形态等，与试验结果吻合较好。陈宇豪等［8］为研究勾头丁坝的水流特性，选用Flow－3D中的RNGk－ε紊流模型建立了勾头丁坝水流计算三维数学模型，并采用淹没式丁坝的物理模型试验验证了所建模型能较好地模拟过流丁坝水流。王小明等［9］基于OpenFOAM开源程序包模拟了不同透空率潜坝的三维流场，研究分析了透水潜坝周围三维水流特性受透空率和流量的影响规律。

堆积体的存在会对下游水电站的安全运行造成威胁，而其工程处置的首要问题应是明晰堆积体附近的水流特性，淹没丁坝的研究思路可为本研究提供参考。笔者采用数值模拟技术对淹没堆积体附近水流的水力特性进行研究，以期为后续关于堆积体泥沙起动及泥沙输移的试验研究打下基础，并为堆积体影响下峡谷地区河道治理方案的优化提供理论支持。

1 数值模型及其验证

1.1 控制方程

本研究采用标准k－ε湍流模型，其水流连续性方程为

动量守恒方程为

湍流模型控制方程为

式中：ρ为水流密度，kg／m3；t为时间，s；xi和xj分别代表笛卡儿坐标系的x、y轴；ui和uj分别为i、j方向上的流速，m／s；μ为液体动力黏滞系数；Si为动量守恒方程源项；P为压力，N；k为湍动能；ε为湍动能耗散率；μt为湍动黏度；Gk为由平均速度梯度得到的湍动能的产生项；Gb为浮力引起的湍动能的产生项；YM为可压湍流中脉动扩张的作用；G1ε、G2ε、G3ε为经验常数；σk、σε为与k和ε对应的Prandtl数；Sk、Sε为用户定义的源项。

1.2 网格划分及边界条件

根据试验水槽尺寸，在ICEM中建模和划分网格，图1为模型的三维视图。堆积体由上、下游对称面和侧面构成，上、下游面均采用半圆锥体结构，其平面投影为扇形，堆积体模型侧向坡度S＝45°，沿水流方向上的最大长度l＝50 cm，横向宽度b＝25 cm，高度Hcon＝10 cm，进占比为b／B＝0.5，采用非结构性四面体网格进行网格划分，将堆积体模型周围局部加密，网格总数为85万，网格划分结果如图2、图3所示。本文采用速度和压力进口边界，出口边界选择压力出口边界，采用VOF（Volume of Fluid）法计算自由水面。对于水气两相流，网格内水流所占比例采用体积分数α表示，在每个单元体中，水和空气的体积分数之和为1。速度进口由进口面的水相进口表示，压力进口由进口面的气相进口表示（设为恒压）。本文采用SIMPLE算法进行求解，即水的体积分数α可由耦合格式计算。

9月25日下午，中欧水资源交流平台地下水管理专题会议召开。水利部副部长胡四一、丹麦环境大臣奥肯出席开幕式并致辞。

图1 三维模型

图2 计算区域网格划分结果

图3 堆积体附近局部网格加密情况

1.3 模型验证

验证试验在西华大学流体及动力机械教育部重点实验室进行。试验的供水系统为循环式供水系统，矩形水槽长L＝18 m，宽B＝0.5 m，高H＝0.42 m，水槽底部坡降为0.15%，水深用LH－1自动水位仪测量，流速采用螺旋测速仪测量。

试验流量Q＝35.8 L／s，测量布置如图4所示，x、y、z分别为横向、纵向、垂向坐标。研究采用以下无量纲化的处理方式：x∗＝x／l和y∗＝y／l分别表示断面和测线的相对位置，以z∗＝z／h0代表z方向的相对位置，u∗＝u／u0代表相对流速。

图4 水槽试验测量布置

选取y∗＝0.3、y∗＝0.5、y∗＝0.7三条测线水深的横向变化进行对比，实测值与模拟值如图5所示，可以看出两者较为吻合，平均误差为2.88%，堆积体下游回流区附近的误差相对较大。选取了断面x∗＝－0.5、x∗＝0及x∗＝0.5的中心测线（y∗＝0. 5）的纵向流速的垂向变化进行对比分析，结果如图6所示，平均误差为1.13%，同样，堆积体下游回流区附近的误差相对较大。其原因是在回流区伴随大量气泡、涡旋的产生，导致水位和流速的测量结果误差相对偏大。总体来说，模拟值与实测值整体较为吻合，水深的横向变化和纵向流速的沿程变化规律与水槽试验结果一致。本文淹没堆积体数值模拟研究的变量包括流量、进占比和坡度，计算工况见表1。

表1 模拟计算工况

图5 纵向水面线验证结果

图6 纵向流速验证结果

参考前人关于淹没丁坝的研究，笔者引入δ（淹没系数）作为堆积体处于淹没状态时衡量其淹没程度的指标，定义如下：

式中：h为对应工况下堆积体顶部中轴线x∗＝0断面的水深，m；h0为对应工况下的均匀流水深，m。

2 水面线分析

从3种淹没系数（小、中、大）中选择具有代表性的工况进行分析，其中流量因素以D1、D2、D3为例，进占比因素以M1、M2、D3为例，侧向坡度因素以N1、N2、N3为例，以所有测线的平均水深代表各横断面的水深，计算不同工况下沿程相对水深变化情况，如图7所示。

图7 断面平均水深沿程分布

整体上看，3种不同因素影响下的水面线沿程变化趋势较为一致。在堆积体上游段，水面线因受到堆积体的阻水作用而壅高，水面线较为平稳。进入堆积体渐缩段，水面线下降较为迅速，水流的惯性作用使水深在x∗＝0.3处最小。水流经过堆积体中轴线后，过流断面逐渐恢复，水面线开始缓慢回升。堆积体对水流的影响越往下游越小，水深逐渐恢复到均匀流水深。由此可以看出，处于淹没状态时，堆积体附近水面线沿程变化的一般规律与非淹没状态下无明显差异［10］。比较分析不同因素影响下的水面线可以发现，上游水面线的壅水高度与流量成反比，与进占比和坡度成正比。从水面线变化幅度来看，堆积体进占比、侧向坡度、流量3种因素的影响程度依次加强。与非淹没情况相比，堆积体侧向坡度因素的阻流影响较强，而进占比的影响较弱。

相对壅水水深h′以x∗＝－4.5～－0.5各断面水深h与均匀流水深h0的比值表示，即h′＝h／h0，h′越大，则上游壅水越明显。相对最小水深h′min以断面x∗＝0.3的平均水深与均匀流水深的比值表示，h′min越小，则堆积体段的水深下降幅度越大。将各工况条件下的淹没系数与相对壅水水深及相对最小水深进行数据拟合，结果表明淹没系数与相对壅水水深及相对最小水深具有良好的指数关系：

由图8、图9可以看出，淹没系数越大，相对壅水水深越小，相对最小水深越大，且淹没系数对相对壅水水深的影响随着淹没系数的增大而减弱。由此可知，堆积体处于淹没状态时，对堆积体顶部以下的水流仍然起到调节作用，其调节作用随着淹没系数的增大而减小，这与丁坝处于淹没状态时的调节机理相似［11］。当淹没系数足够大时，可以认为水面线几乎不受堆积体影响。

图8 淹没系数与相对壅水水深关系曲线

图9 淹没系数与相对最小水深关系曲线

3 流场特性分析

3.1 紊动特性分析

水流的紊动强度对河道水流挟沙能力的影响较大，本文用紊动能来表征水流整体的紊动情况，其计算公式为

式中：TKE为紊动能；u、v、w分别为横向、纵向、竖向的脉动流速。

以工况D2为例，分析不同流层的紊动能分布情况，选取底层（z∗＝0.02）、中层（z∗＝0.51）和上层（z∗＝0.75）为代表，分析不同流层的紊动能分布规律，如图10所示。整体上看，堆积体上游壅水区的紊动能数值均较小，此区域容易造成泥沙淤积。底部过流断面进占程度较高，水流在近底层受到明显压缩，水流扩散导致的高紊动区出现在底部，且随流层升高，水流的扩散作用减弱，紊动强度也随之减弱。下游出现的高紊动区主要是主流扩散导致的，剪切流所引起的高紊动虽然存在，但和水流扩散的高紊动连成一片，之后随着流层升高，回流现象开始出现［12］，剪切流所引起的高紊动才更加清晰地被观察到。随着回流现象的出现，高紊动区域的分布和强度开始受到回流现象的主导，同时，堆积体束窄作用减弱，水流扩散导致的高紊动已不再明显，高紊动区已逐渐转移到回流区附近，随着回流影响范围的增大，其紊动能也逐渐增大，越往下游紊动强度越弱。由此可知，高紊动区受到回流的影响会发生改变，流层越往上，受回流的影响越大。泥沙运动与水流紊动强度有密切的关系，泥沙发生运动，河床剪切力并不是唯一的影响因素，而紊动引起的上举力会促使水流进一步挟带和运移河床表面的泥沙颗粒［13］，因此在回流区更易发生泥沙输移，且水流挟带泥沙的能力更强。

图10 不同水深紊动能分布

3.2 涡量分析

涡系的存在影响了水流中物质的输移特性且影响水流的能量耗散过程，研究堆积体附近的涡量场使得对绕流结构体作用下的冲淤特性具有更加清晰的认识。根据文献［14］推导的涡量表达式，在竖直方向上涡量采用如下公式计算：

式中：ux、uy分别为横向、纵向的时均流速。

图11为工况D1下z∗＝0.15（底层）、z∗＝0.45（中层）和z∗＝0.75（上层）3个流层水平剖面的涡量分布。整体上可以看出，涡量较大的区域均产生于堆积体表面和回流区与主流区的交界段，这与淹没丁坝坝后水流结构PIV试验研究［15］得出的结论一致。涡量大是形成冲坑的主要原因［16］，因此在堆积体表面和回流区与主流区的交界段容易造成较为严重的冲刷。堆积体的束窄作用使绕流处存在较大的流速梯度，因此堆积体表面的正涡值较大。堆积体“下大上小”的形态特征导致水流在近底层（见图11（a））受到明显压缩，流速逐渐增大，主流与下游低流速区的流速梯度引起剪切流，主流区水流带动回流区水流向下游运动，产生了向下游扩散的负涡。由于底层尚未出现回流，因此下游沿水流方向涡量相对较小。如图11（b）、（c）所示，随着流层的升高，受到主流区和回流区水流之间强烈的剪切作用，回流现象逐渐明显［12］，因此该区域形成了正涡值较大的涡量带。同时，堆积体束窄作用减弱，堆积体表面涡量影响范围也随之减小。两岸边壁的正负涡量绝对值虽不及上述两区域大，但也会对岸坡造成冲刷影响。文献［14］指出不可压缩流体中能量耗散率与流体中涡量绝对值的大小直接相关，它与涡量的平方成正比，由此可知水流在经过堆积体与绕流的交界处和回流区与主流区的交界段后，其能量会迅速耗散。