• 全国中文核心期刊
  • 中国科技核心期刊
  • 美国工程索引(EI)收录期刊

一维-三维耦合的明渠瞬变流模拟方法与应用

吴家阳, 李安强, 程永光

吴家阳, 李安强, 程永光. 一维-三维耦合的明渠瞬变流模拟方法与应用[J]. 水科学进展, 2022, 33(6): 989-999. DOI: 10.14042/j.cnki.32.1309.2022.06.013
引用本文: 吴家阳, 李安强, 程永光. 一维-三维耦合的明渠瞬变流模拟方法与应用[J]. 水科学进展, 2022, 33(6): 989-999. DOI: 10.14042/j.cnki.32.1309.2022.06.013
WU Jiayang, LI Anqiang, CHENG Yongguang. 1D-3D coupling method for transient flow simulations in open channel and its applications[J]. Advances in Water Science, 2022, 33(6): 989-999. DOI: 10.14042/j.cnki.32.1309.2022.06.013
Citation: WU Jiayang, LI Anqiang, CHENG Yongguang. 1D-3D coupling method for transient flow simulations in open channel and its applications[J]. Advances in Water Science, 2022, 33(6): 989-999. DOI: 10.14042/j.cnki.32.1309.2022.06.013

一维-三维耦合的明渠瞬变流模拟方法与应用

基金项目: 

湖北省重点研发计划资助项目 2020BCA087

国家重点研发计划资助项目 2021YFC3000102087

详细信息
    作者简介:

    吴家阳(1989-), 男, 湖南长沙人, 高级工程师, 博士, 主要从事计算流体动力学、高性能计算方面研究。E-mail: wujiayang@cjwsjy.com.cn

  • 中图分类号: TV131

1D-3D coupling method for transient flow simulations in open channel and its applications

Funds: 

the Key R&D Program of Hubei Province, China 2020BCA087

the National Key R&D Program of China 2021YFC3000102087

  • 摘要: 为精细化模拟明渠中水电站瞬变过程引发的涌浪, 克服以往低维模型无法细致刻画涌浪的横向和垂向运动、以及高维模型计算效率偏低等问题, 基于水位预测矫正法, 提出了一维-三维耦合的明渠瞬变流模拟方法。通过潮汐波浪运动、局部溃坝波和顺直明渠瞬变流算例验证了方法的有效性和准确性, 模拟结果与渐进理论解、商用软件模拟结果等吻合较好。将该方法应用在实际工程中, 计算了葛洲坝水利枢纽瞬变过程中库区的涌浪过程, 分析了瞬变过程的最大、最小涌浪波高和涌浪的衰减规律, 结果与频谱分析结果吻合较好, 说明了本方法能准确模拟实际的明渠瞬变流过程。
    Abstract: Based on water stage prediction-correction method, an 1D-3D coupled method is proposed to accurately simulate the surge wave caused by the transient process of hydropower stations in open channel, overcoming the inability of low-dimensional models to fully resolve the horizontal and vertical fluctuation of the surge wave, and low computational efficiency of high-dimensional models. Typical verification cases, including the tidal wave motion, the partial dam break flow and the open channel transient flow are simulated to analyze the accuracy of the proposed method. The simulated results are in good agreement with asymptotic theoretical solution and those yielded by commercial software. To futher demonstrate its capability in practical cases, the surge wave triggered by the transient process of the Gezhouba hydropower station is simulated. The maximum and minimum surge wave heights including its attenuation characteristics are analyzed. The results are in good agreement with those yielded by spectrum analysis, showing that the proposed method is capable of simulating practical transient flows in open channel.
  • 布置于明渠上的水电站在发生增、甩负荷等瞬变过程时往往会在上游水库引发较大涌浪,产生机组出力和引航道水位剧烈波动等问题,对水电站的运行安全产生较大危害[1]。采用数值或物理模型手段研究涌浪在明渠内的传播规律和控制措施,对减缓机组瞬变过程中的安全危害意义重大。目前,明渠瞬变流问题大都采用断波法[2]或数值求解一维圣维南方程组和二维浅水方程组[3-4]解决。当瞬变过程发生在较为宽阔的河道型水库(如三峡水库)或湖泊(如鄱阳湖水利枢纽)内时,为了准确刻画机组附近复杂的横向和垂向流动过程,往往需要应用高维模型模拟,反之为了获得较高的计算效率,则需要尽可能采用低维模型计算。不加区分地采用低维或高维模型模拟枢纽附近和库尾河道等不同位置的瞬变过程,既不现实也无必要,亟需根据计算域几何特征和流场特点,有的放矢地选取不同维度的模型模拟不同计算区域的瞬变过程,在不牺牲计算精度的同时保证计算效率。

    多维耦合模型在平原河网、入海河湾和大型河库等区域的水动力学求解中应用较为广泛。Chen等[5]、王智勇等[6]采用水位预测矫正法,构建了一维-二维耦合的河湖系统整体水动力模型,并采用特征理论证明了该方法的收敛性;诸裕良等[7]建立了一种一维、二维全隐河网海湾水动力联网数学模型,采用接口断面法实现一维、二维模型间的耦合;Morales- Hernández等[8]提出了一种适用于明渠浅水流动的一维-二维耦合方案,其中一维、二维模型分别采用隐式有限差分法和有限体积法数值求解;陈文龙等[9]建立了基于侧向联解的一维-二维耦合水动力学模型,有效克服了基于堰流公式的传统方法难以处理模型间动量交换的缺点;顾巍巍等[10]建立了一维-二维耦合的洪水演进数学模型,并采用堰流公式实现一维、二维模型连接处的水流交互;Liu等[11]提出了一种适用于复杂地形和非规则边界的一维-二维水动力学耦合模型,并应用于河道-蓄滞洪区系统的水动力求解;王秀杰等[12]考虑洪水和风暴潮的共同作用,构建了天然河道漫溃堤洪水在防洪保护区的一维-二维水动力耦合模型;程涛等[13]利用MIKE软件构建了适用于不规则计算域和地形的溃坝洪水一维-二维耦合模型。可见,目前多维耦合模型主要应用于一维和二维模型间的耦合求解,而三维模型因其模型实现难度大、计算效率较低和自由液面处理复杂等因素在模型耦合中尚不多见。

    本文提出一种一维-三维耦合的浸没边界-格子玻尔兹曼明渠瞬变流模拟方法。其中,瞬变流场中的狭长河道采用一维圣维南方程组描述,水库、宽阔湖面、机组进水口、泄水闸等高维流动特征明显的水域采用三维Navier-Stokes(N-S)方程组描述;一维圣维南方程组和三维N-S方程组分别采用隐式Pressimann差分方法和格子Boltzmann方法求解,而一维-三维模型间的耦合采用水位预测矫正法[5, 7]实现。

    在假定明渠底坡较缓、渠道断面压力沿垂线按静水压力分布、渠道断面流速均匀分布等条件下,一维明渠水流可由圣维南方程组描述:

    $$ B \frac{\partial Z}{\partial t}+\frac{\partial Q}{\partial x}=q $$ (1)
    $$ \frac{\partial Q}{\partial t}+\frac{\partial}{\partial x}\left(\frac{Q^2}{A}\right)+g A \frac{\partial Z}{\partial x}+g A \frac{n^2 Q|Q|}{A^2 R^{4 / 3}}=0 $$ (2)

    式中:ZQ分别为明渠的水位和流量;AB分别为明渠断面的过水面积和宽度;q为旁侧入流或出流的单宽流量;n为明渠河床糙率系数;R为水力半径;g为重力加速度;x为沿明渠水流方向的河道里程;t为时间。式(1)和式(2)适用于一般非棱柱型明渠的水动力学模拟。

    采用广泛应用的四点Pressimann隐格式求解一维圣维南方程组,时间偏导数取相邻节点时间偏导数的平均值,空间偏导数取相邻节点一阶向前差分的加权平均值。控制方程离散后得到的非线性方程组采用Newton-Raphson迭代法求解。Pressimann隐格式的详细推导过程详见文献[12]。

    三维黏性不可压缩流体的控制方程为N-S方程组:

    $$ \nabla \cdot \boldsymbol{u}=0 $$ (3)
    $$ \frac{\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u} \cdot \nabla \boldsymbol{u}=-\frac{1}{\rho} \nabla p+\nu \Delta \boldsymbol{u}+\boldsymbol{f} $$ (4)

    式中:upf为分别流速、压强、外力;ν为流体的运动黏度;ρ为流体密度。

    本文中,N-S方程组由格子Boltzmann方法数值求解。在格子Boltzmann方法中,最基本的物理量为分布函数fα(x, t),其统计解释为在t时刻x坐标处沿着α方向运动的粒子质量密度。流体密度和动量可由分布函数的各阶矩计算:

    $$ \sum\limits_\alpha f_\alpha(x, t)=\rho $$ (5)
    $$ \sum\limits_\alpha \boldsymbol{e}_\alpha f_\alpha(x, t)=\rho \boldsymbol{u} $$ (6)

    式中:eα为离散速度方向矢量。

    fα(x, t)满足如下离散格子Boltzmann方程:

    $$ f_\alpha\left(x+\boldsymbol{e}_\alpha \delta_t, t+\delta_t\right)-f_\alpha(x, t)=-\frac{1}{\tau}\left[f_\alpha(x, t)-f_\alpha^{\mathrm{eq}}(x, t)\right]+\delta_t \boldsymbol{F}_\alpha $$ (7)

    式中:fαeq(x, t)为平衡态分布函数;τ为松弛时间,与流体的运动黏度有关;δt为时间步长;Fα为外力项。在三维情形下,常用的离散速度模型为D3Q19模型,其离散速度定义详见文献[14]。

    平衡态分布函数一般取为麦克斯韦统计分布,D3Q19离散速度模型的平衡态分布函数可表示为

    $$ f_\alpha^{\mathrm{eq}}=\omega_\alpha\left[\rho+3 \boldsymbol{e}_\alpha \cdot \boldsymbol{u}+\frac{9\left(\boldsymbol{e}_\alpha \cdot \boldsymbol{u}\right)^2}{2}-\frac{3 \boldsymbol{u} \cdot \boldsymbol{u}}{2}\right] $$ (8)

    式中:ωα为权系数,取值可参考文献[14]。

    为了追踪自由液面,参考文献[15]中的方法,将三维网格点划分为流体节点、气体节点和液面节点3类,其中,流体节点完全被水体填充,气体节点完全被气体填充,而液面节点可视为被水体部分填充。液面节点在计算过程中的任一时刻均能形成封闭的曲面,将流体节点和气体节点完全分隔开,如图 1所示。追踪自由液面的运动过程分为追踪液面节点的运动、处理液面节点边界条件和更新流场节点类别3个步骤。

    图  1  三维格子Boltzmann方法中3类网格节点示意
    Fig.  1  Schematic diagram of three types of grid node in three-dimensional lattice Boltzmann method

    本文通过计算每一网格点所含水体质量实现自由液面的追踪。对于每一网格点,定义其质量(m)和体积分数(ε),体积分数定义为网格点所含水体质量与水体密度之比(ε=m/ρ)。根据该定义,气体节点ε=0,流体节点ε=1,自由液面节点0 < ε < 1。与Volume-of-Fluid(VOF)方法类似,液面追踪通过计算任意网格点与周围网格点的质量交换实现。考虑到格子Boltzmann方法中分布函数fα(x, t)的物理意义,可以利用分布函数在相邻网格点间的迁移实现质量交换。基于上述思想,对于tt时刻位于x位置的网格点,相邻网格点既可能是流体节点,也可能是自由液面节点,当其与相邻流体节点进行质量交换时,沿x方向的质量变化(Δmα(xtt))可表示为

    $$ \Delta m_\alpha(x, t+\Delta t)=f_{\bar{\alpha}}\left(x+\boldsymbol{e}_\alpha, t\right)-f_\alpha(x, t) $$ (9)

    式中:αα的反方向。式(9)中等号右端第一项为沿α方向进入网格点的流体质量,第二项为沿α方向离开网格点的流体质量。当其与相邻自由液面节点进行质量交换时,则需要考虑自由液面节点被水体部分填充的事实,将体积分数考虑在内,质量变化可表示为

    $$ \Delta m_\alpha(x, t+\Delta t)=\left[f_\alpha\left(x+\boldsymbol{e}_\alpha, t\right)-f_\alpha(x, t)\right] \frac{\varepsilon\left(x+\boldsymbol{e}_\alpha, t\right)-\varepsilon(x, t)}{2} $$ (10)

    式(9)和(10)均为对称结构,若坐标x的网格点沿α方向获取了质量Δmα,则坐标x+eα的网格点将沿α方向损失相同质量Δmα。因为气体网格点与其他类型网格点间不存在质量交换,所以式(9)和(10)能确保整个系统质量守恒。对于tt时刻位于x位置的自由液面节点,其质量变化(m(xtt))为Δmα(xtt)沿D3Q19模型所有离散速度方向的和:

    $$ m(x, t+\Delta t)=m(x, t)+\sum \Delta m_\alpha(x, t+\Delta t) $$ (11)

    式中:m(x, tt)为自由液面节点下一时步的质量,可根据其计算结果实现节点类别的转变。若m(x, tt)>(1+κ)ρ(x, tt),则转变为流体节点;若m(x, tt) < -κρ(x, tt),则转变为气体节点。其中,阈值κ=10-3,以防止新生成的自由液面节点在下一时步再次改变节点类型。

    需要说明的是,对于新生成的流体节点,其质量可能大于水体密度,导致ε>1;对于新生成的气体节点,其质量可能为负值,导致ε < 0。因此,需进行质量的二次分配,在保证系统质量守恒前提下,确保所有网格点的体积分数在[0, 1]区间内。对于上述2类特殊节点,需进行二次分配的质量(mex)分别为mex=mρmex=m,为了在二次质量分配后不再出现体积分数超出[0, 1]现象,mex将均匀分配给周围新生成的自由液面节点。

    本文将水位预测矫正法[5, 7]推广至一维和三维模型间的耦合。如图 2所示,若某一汊点有M个一维和三维河段与之相连,则在该汊点需满足如下连接条件:

    图  2  一维、三维模型耦合示意
    Fig.  2  Schematic diagram of one-dimensional and three-dimensional model coupling
    $$ \sum\limits_{i=1}^M Q_i=0 $$ (12)
    $$ Z_1=Z_2=\cdots=Z_{M-1}=Z_M $$ (13)

    式中:Qi为流入和流出汊点的净流量(流入为正,流出为负);Zi为与汊点相连的第i个河段的水位,i=1, 2,…,M

    可见,汊点连接条件实质上为各连接河段提供边界条件。在水位预测矫正方法中,所有模型在汊点处均采用水位边界条件。若整个明渠系统在n时刻的水位、流量已知,则可假定汊点在n+1时刻的水位为Z′,从而计算整个明渠系统在n+1时刻的水位、流量,进而计算出流入、流出汊点的净流量。对于非恒定流,汊点净流量为Z′的函数,若能找到某一迭代方法逐步修正Z′以使得汊点净流量趋于0,则可实现一维、三维模型间的耦合。根据文献[10-11],Z′可按式(14)修正:

    $$ \Delta Z^{\prime}=\frac{\Delta t \sum\limits_{i=1}^M Q^{\prime}{ }_i}{A_{\mathrm{c}}} $$ (14)

    式中:ΔZ′为水位矫正值; Qi为流入和流出汊点的流量; Ac为控制迭代收敛速度的系数。水位预测矫正法的收敛性证明及Ac建议取值范围详见文献[7]。

    若给定n+1时刻耦合节点的水位初始估计值(Zkk为迭代步),一维、三维模型耦合步骤总结如下(δ为某一小量):

    (1) 将明渠流场演化至k迭代步,并计算所有耦合节点的净流量;

    (2) 针对某一耦合汊点,若汊点净流量小于δ,则停止该耦合节点的迭代,耦合汊点水位边界条件取为Zk,否则执行步骤(3);

    (3) 根据式(14)计算所有耦合汊点的ΔZ′;

    (4) 矫正所有耦合汊点的水位,并更新与耦合汊点相连的所有河段的水位边界条件;

    (5) 将迭代步更新为k=k+1,重复执行步骤(1)。

    采用潮汐波浪运动、局部溃坝波和顺直明渠瞬变流3个算例分别验证一维圣维南方程组求解方案、三维自由液面N-S方程组求解方法和一维-三维耦合的河库瞬变流模拟方法的有效性和准确性。

    潮汐所引起的波浪运动在海岸工程中十分常见,本小节采用顺直河槽中的潮汐波浪运动算例验证一维圣维南方程组数值解法的有效性和准确性。本例中,顺直明渠总长L=14 000 m,渠底高程(H(x))由式(15)给定:

    $$ H(x)=10+\frac{40 x}{L}+10 \sin \left[\pi\left(\frac{4 x}{L}-\frac{1}{2}\right)\right] $$ (15)

    式中:x为从河槽左侧端点起算的明渠里程。计算初始条件和边界条件分别由式(16)和式(17)给定:

    $$ h(x, 0)=H(x) \quad u(x, 0)=0 $$ (16)
    $$ h(0, t)=H(0)+4-4 \sin \left[\pi\left(\frac{4 t}{86400}+\frac{1}{2}\right)\right] \quad u(L, t)=0 $$ (17)

    式中:h(x, t)和u(x, t)分别为明渠水深和断面平均流速。在上述条件下,明渠内潮汐波浪的波长足够短,本例存在如下渐进理论解。

    $$ h(x, t)=H(x)+4-4 \sin \left[\pi\left(\frac{4 t}{86400}+\frac{1}{2}\right)\right] $$ (18)
    $$ u(x, t)=\frac{(x-14000) \pi}{5400 h(x, t)} \cos \left[\pi\left(\frac{4 t}{86400}+\frac{1}{2}\right)\right] $$ (19)

    本例中,计算空间步长Δx=70 m,时间步长Δt=10 s。t=9 120 s时刻明渠沿程水深和流速的计算结果以及渐进理论解如图 3所示。可见,模拟结果与渐进理论解吻合程度较好,沿程水深的最大统计误差均在1%以内,沿程断面平均流速的最大统计误差也在5%以内,从而验证了一维圣维南方程组四点Pressimann隐式差分解法的有效性和准确性。

    图  3  潮汐波浪运动t=9 120 s计算结果与渐进理论解的比较
    Fig.  3  Comparison between the simulated results and the asymptotic theoretical solutions of tidal wave motion at t=9 120 s

    采用三维局部溃坝算例验证三维N-S方程组格子Botlzmann数值解法的有效性和准确性,该算例广泛应用于MacCormack格式、Pressimann隐格式等明渠非恒定流数值解法的验证。计算条件设置如图 4所示,计算域为200 m×200 m的正方形区域,河床床底高程均为0 m,初始时刻被概化为一条直线的竖直坝体布置在计算域的中心线位置,计算域上游和下游水深分别为10 m和5 m。计算开始时,坝体在图 4所示的位置发生溃坝,溃坝段长度为75 m,溃口一直延伸至河床床底。计算参数如下:重力加速度g=9.8 m/s2,水的运动黏度ν=1.0×10-6 m2/s,水体密度ρ=1.0×103 kg/m3

    图  4  局部溃坝波算例计算条件设置
    Fig.  4  Simulation setups of partial dam break wave case

    图 5为不同时刻水深沿顺河向溃坝段中心线的分布,图 6为不同时刻溃坝自由液面的形态。为了与其他数值格式计算结果进行对比,图 5也绘制了液面梯度法[16]和二维浅水格子Boltzmann方法(SWE-LBM)[17]的计算结果。可见,本文方法与其他数值格式的计算结果总体吻合较好,考虑到液面梯度法忽略了湍流效应,并忽略了控制方程中的扩散项,导致波前存在一定程度的不吻合现象。

    图  5  溃坝段中心线的沿程水深分布
    Fig.  5  Water depth distribution along the center line of dam break breach
    图  6  不同时刻局部溃坝波自由液面形态
    Fig.  6  Free surface configuration of partial dam break wave case at different time

    为验证一维-三维耦合的河库瞬变流模拟方法的有效性和准确性,本节计算简单顺直明渠中的瞬变流过程。计算域为一假想水电站的尾水矩形明渠,明渠长和宽分别为100 m和10 m,渠底底坡为0。尾水渠道下游水深恒定为15 m,上游取为流量边界条件,且流量过程由水电站机组导叶控制。本例模拟水电站机组增、甩负荷过程中尾水明渠内的涌波过程,并与Delft 3D软件的模拟结果进行对比。在水电站机组增、甩负荷过程中,假设机组导叶按照线性规律开启或关闭,开启或关闭历时均为8 s,机组额定工况条件下应用流量为600 m3/s。增负荷工况下,初始时刻尾水明渠水体维持静止,之后导叶开启,瞬变过程开始。直至尾水明渠内水流达到恒定流状态,该恒定流状态也将作为甩负荷工况的计算初始条件。

    为了验证一维-三维耦合的河库瞬变流模拟方法,本例将机组尾水明渠一分为二,上游50 m渠段采用一维圣维南方程组的Pressimann隐式差分法计算,下游50 m渠段采用三维N-S方程组的格子Boltzmann方法计算。一维、三维模型计算空间步长分别为Δx=0.25 m和Δx=0.20 m,计算时间步长均为t=0.01 s。

    图 7为水电站机组增、甩负荷过程中上游水位和下游流量的历时曲线。可见,本文方法计算结果与Delft 3D软件模拟结果总体吻合较好,且本文方法捕捉涌波陡峭拐点的能力略强于商用软件,说明其具有较小的数值扩散性。

    图  7  机组尾水明渠上游水位和下游流量的历时曲线
    Fig.  7  Histories of upstream water level and downstream discharge of tailrace open channel

    葛洲坝水利枢纽兴建于20世纪七八十年代,是一座低水头、大流量、径流式水电站。作为上游三峡水利枢纽工程的反调节航运枢纽,与三峡枢纽工程联合运行,可消除三峡电站日调节时的下游不稳定流,以改善河道航运条件,并利用该河段水位差发电。

    由于泥沙淤积,在葛洲坝枢纽坝址形成了2座岛屿,将坝址江面分为大江、二江和三江3条支汊。枢纽工程共布置21台机组,其中,大江厂房14台,额定流量均为825 m3/s; 余下7台机组布置在二江厂房,包括2台额定流量为1 130 m3/s和5台额定流量为825 m3/s的机组; 三江仅具备冲沙和航运功能。机组以额定工况恒定运行时,库前水位维持正常蓄水位66 m。

    本例模拟大江14台机组同时甩负荷而二江7台机组以额定流量运行时,上游江面的涌浪过程。大江机组甩负荷时间为33 s,机组应用流量假设按线性规律从满发流量11 550 m3/s(14×825 m3/s)减小至0 m3/s。

    计算范围为葛洲坝枢纽坝址至三峡枢纽工程共40 km的库区和河道,其中葛洲坝枢纽坝址上溯5 km范围采用三维N-S方程组的格子Boltzmann方法模拟,而余下的35 km河道则采用一维圣维南方程组的Pressimann隐式差分法计算。上游边界设为定流量边界,下游边界设为流量边界条件,且出流过程由21台机组的应用流量控制。为简化起见,本例不考虑枢纽冲沙闸和船闸过流能力,且认为机组进水口之间没有间隙,同时不考虑支流及区间入汇。河道糙率一律取为0.012,计算时间步长取为0.172 8 s,一维、三维模型的计算空间步长分别为11.67 m和8.64 m,网格总数达到千万量级,恒定流计算耗时14.4 h。

    为了对葛洲坝水库各处的水位进行监测,本文将在图 8所示位置布置流场测点,其中,A测点位于大江冲沙闸前;B点位于大江机组进水口前;C点位于二江泄水闸前;D点位于二江机组进水口前;E点位于三江冲沙闸前;F点位于大江防洪淤堤上游端点;G点位于三江防洪淤堤上游端点;H点位于支流河口。

    图  8  流场监测点布置位置示意
    Fig.  8  Layout of flow field monitoring points

    各时刻葛洲坝水库水位和水平流速计算结果如图 9(分图左为库水位分布,分图右为水平流速分布)所示。由于江面起伏相较于计算域水平尺度太小,为了更好地显示江面起伏情况,图中水位均沿纵向放大了100倍。可见,大江机组甩负荷后,涌浪随即以机组为原点以弧形方式向四周传播。当t=1.0 min时,涌浪首次来到三江防洪淤堤的右岸并发生反射,致使二江机组进水口水位开始上涨。随后,涌浪向大江防洪淤堤传播,并在随后一段时间内在在大江和三江防洪淤堤之间往复来回行进,该现象在模拟初期尤为明显,在二江左右岸间存在明显的水位高低交替变化。

    图  9  大江机组甩负荷后库区流场变化示意
    Fig.  9  Flow field in the load rejection process of Dajiang hydro-plant units

    涌浪在大江、三江防洪淤堤间往复传播的同时也在向上游行进。t=2.0 min时,大江冲沙闸首次受到涌波的影响,由于大江河道较为狭窄,其水位上涨非常明显。t=8.0 min时,大江冲沙闸水位首次到达波谷;t=3.5 min和t=4.5 min时,大江和二江机组分别达到波谷。

    图 10为各流场监测点的水位历时曲线,其中时间按小时计,且自大江机组甩负荷开始计时。不同频率涌浪波动之间相互叠加的现象在图中清晰可见。从图 10中可以得到枢纽各部位受涌浪影响的起始时间、最大和最小涌浪波高、涌浪完全衰减耗时等对枢纽安全运行至关重要的参数。例如,对于三江冲沙闸,最大和最小涌浪分别出现在机组甩负荷后1.2 h和2.55 h左右,最大和最小涌浪水位分别为68.0 m和65.7 m。

    图  10  各流场监测点水位过程
    Fig.  10  Water level histories of flow field monitoring points

    图 11将大江和二江机组进水口水位进行叠加,高频振动成分在图中清晰可见,大江水位的波峰和波谷正好分别对应二江水位的波谷和波峰。如前所述,该现象源自于在大江、三江防洪淤堤间往复传播的涌浪,因为二江左右岸距离较短,所以该振动以高频形式附加在低频基波上,本例中的低频基波源自于大江机组进水口和上游流量边界间往复传播的涌浪。

    图  11  大江、二江机组进水口水位—历时曲线
    Fig.  11  Water level histories of Dajiang and Erjiang hydro-plant unit′s inlets

    利用快速傅里叶变换对大江机组、二江机组、大江冲沙闸和三江冲沙闸的水位—历时曲线进行频谱分析,区分水位波动的各种频率成分并对其成因进行解释。频谱分析结果表明,大江机组、二江机组、大江冲沙闸均包含4种频率成分,分别为1.92×10-4 Hz(频率1)、8.32×10-4 Hz(频率2)、2.50×10-3 Hz(频率3)、7.78×10-3 Hz(频率4);三江冲沙闸包含2种频率成分,分别为1.92×10-4 Hz(频率1)、8.32×10-4 Hz(频率2)。本例中河道水深大致在28~31 m范围内波动,根据波速公式可得涌浪传播速度约为16 m/s。大江进水口至上游边界和支流顶端边界的距离分别为40 km和9 km左右。根据波速公式,往返于两者间的涌浪周期分别为5 000 s和1 125 s,对应频率分别为2.00×10-4 Hz和8.89×10-4 Hz,而这与频率1和频率2基本吻合; 其次,大江进水口至三江冲沙闸的距离为3.55 km,往返于其间的涌浪周期为443.7 s,对应频率为2.25×10-3 Hz,这与频率3基本一致;最后,对于往返于大江、三江防洪淤堤间的涌浪(淤堤间距1.05 km),其传播周期为131.3 s,对应频率为7.62×10-3 Hz,这与频率4基本一致。可见,计算结果基本符合物理规律,说明一维-三维耦合的河库瞬变流模拟方法基本具备模拟实际明渠瞬变流的能力。

    采用水位预测矫正法,提出了一种一维-三维耦合的明渠瞬变流计算模型,其中一维模型采用传统有限差分法计算,三维模型采用基于自由液面的格子Boltzmann方法模拟,主要结论如下:

    (1) 本模型可模拟具有任意地形、边界条件的实际明渠瞬变流过程,可计算瞬变流过程中的水位、流量和流速等水力变量。

    (2) 根据模拟需求,通过在计算域的不同位置运用不同维度的模型进行计算,本模型在不损失整体计算精度的前提下,大幅提高实际明渠瞬变流模拟效率。

    (3) 采用潮汐波浪运动、局部溃坝波和顺直明渠瞬变流等算例验证了方法的有效性和正确性,模拟结果与渐进理论解、商用软件模拟结果均吻合较好。

    (4) 将本文方法应用在实际工程中,计算了葛洲坝水利枢纽瞬变过程中库区的涌浪过程,分析了瞬变过程中的最大、最小涌浪波高和涌浪的衰减规律。通过对涌浪水位过程进行频谱分析,得到枢纽工程各部位水位波动的频率组成,并解释了各频率成分的形成原因,说明了本方法能准确模拟实际工况下的明渠瞬变流过程。

    今后应重点研究本方法的并行加速策略,以大幅提高实际问题的模拟效率。

  • 图  1   三维格子Boltzmann方法中3类网格节点示意

    Fig.  1   Schematic diagram of three types of grid node in three-dimensional lattice Boltzmann method

    图  2   一维、三维模型耦合示意

    Fig.  2   Schematic diagram of one-dimensional and three-dimensional model coupling

    图  3   潮汐波浪运动t=9 120 s计算结果与渐进理论解的比较

    Fig.  3   Comparison between the simulated results and the asymptotic theoretical solutions of tidal wave motion at t=9 120 s

    图  4   局部溃坝波算例计算条件设置

    Fig.  4   Simulation setups of partial dam break wave case

    图  5   溃坝段中心线的沿程水深分布

    Fig.  5   Water depth distribution along the center line of dam break breach

    图  6   不同时刻局部溃坝波自由液面形态

    Fig.  6   Free surface configuration of partial dam break wave case at different time

    图  7   机组尾水明渠上游水位和下游流量的历时曲线

    Fig.  7   Histories of upstream water level and downstream discharge of tailrace open channel

    图  8   流场监测点布置位置示意

    Fig.  8   Layout of flow field monitoring points

    图  9   大江机组甩负荷后库区流场变化示意

    Fig.  9   Flow field in the load rejection process of Dajiang hydro-plant units

    图  10   各流场监测点水位过程

    Fig.  10   Water level histories of flow field monitoring points

    图  11   大江、二江机组进水口水位—历时曲线

    Fig.  11   Water level histories of Dajiang and Erjiang hydro-plant unit′s inlets

  • [1] 陈慧敏, 陈龙, 郑国晨, 等. 多江汊上河床式电站过渡过程涌波特性[J]. 水电能源科学, 2016, 34(9): 173-177. https://www.cnki.com.cn/Article/CJFDTOTAL-SDNY201609043.htm

    CHEN H M, CHEN L, ZHENG G C, et al. Surge waves characteristics of water retaining type hydropower station on multiple channel river[J]. Water Resources and Power, 2016, 34(9): 173-177. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SDNY201609043.htm

    [2] 王刚, 刘沛清. 无压瞬变流的改进断波法及其应用[J]. 长江科学院院报, 2001, 18(3): 3-6. https://www.cnki.com.cn/Article/CJFDTOTAL-CJKB200103000.htm

    WANG G, LIU P Q. Improved segment wave algorithm of unsteady flow in open channel and its application[J]. Journal of Yangtze River Scientific Research Institute, 2001, 18(3): 3-6. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-CJKB200103000.htm

    [3]

    TSENG M H. Improved treatment of source terms in TVD scheme for shallow water equations[J]. Advances in Water Resources, 2004, 27(6): 617-629. doi: 10.1016/j.advwatres.2004.02.023

    [4]

    DELIS A I, SKEELS C P. TVD schemes for open channel flow[J]. International Journal for Numerical Methods in Fluids, 2015, 26(7): 791-809.

    [5]

    CHEN Y C, WANG Z Y, LIU Z W, et al. 1D-2D coupled numerical model for shallow-water flows[J]. Journal of Hydraulic Engineering, 2012, 138(2): 122-132. doi: 10.1061/(ASCE)HY.1943-7900.0000481

    [6] 王智勇, 陈永灿, 朱德军, 等. 一维-二维耦合的河湖系统整体水动力模型[J]. 水科学进展, 2011, 22(4): 516-522. doi: 10.14042/j.cnki.32.1309.2011.04.010

    WANG Z Y, CHEN Y C, ZHU D J, et al. 1D-2D coupled hydrodynamic simulation model of river-lake system[J]. Advances in Water Science, 2011, 22(4): 516-522. (in Chinese) doi: 10.14042/j.cnki.32.1309.2011.04.010

    [7] 诸裕良, 严以新, 李瑞杰, 等. 河网海湾水动力联网数学模型[J]. 水科学进展, 2003, 14(2): 131-135. doi: 10.14042/j.cnki.32.1309.2003.02.001

    ZHU Y L, YAN Y X, LI R J, et al. A coupled hydrodynamic numerical model for bay-river network[J]. Advances in Water Science, 2003, 14(2): 131-135. (in Chinese) doi: 10.14042/j.cnki.32.1309.2003.02.001

    [8]

    MORALES-HERNÁNDEZ M, GARCÍA-NAVARRO P, BURGUETE J, et al. A conservative strategy to couple 1D and 2D models for shallow water flow simulation[J]. Computers & Fluids, 2013, 81: 26-44. doi: 10.3969/j.issn.1001-3695.2013.01.005

    [9] 陈文龙, 宋利祥, 邢领航, 等. 一维-二维耦合的防洪保护区洪水演进数学模型[J]. 水科学进展, 2014, 25(6): 848-855. http://skxjz.nhri.cn/article/id/2474

    CHEN W L, SONG L X, XING L H, et al. A 1D-2D coupled mathematical model for numerical simulating of flood routine in flood protected zone[J]. Advances in Water Science, 2014, 25(6): 848-855. (in Chinese) http://skxjz.nhri.cn/article/id/2474

    [10] 顾巍巍, 江雨田, 张卫国, 等. 一维-二维耦合的水库下游地区洪水演进数学模型及应用[J]. 水电能源科学, 2018, 36(5): 43-45, 132. https://www.cnki.com.cn/Article/CJFDTOTAL-SDNY201805011.htm

    GU W W, JIANG Y T, ZHANG W G, et al. 1D-2D coupled mathematical model for numerical simulation of flood routing in downstream of reservoir and its application[J]. Water Resources and Power, 2018, 36(5): 43-45, 132. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SDNY201805011.htm

    [11]

    LIU Q, QIN Y, ZHANG Y, et al. A coupled 1D-2D hydrodynamic model for flood simulation in flood detention basin[J]. Natural Hazards, 2015, 75(2): 1303-1325. doi: 10.1007/s11069-014-1373-3

    [12] 王秀杰, 胡冰, 苑希民, 等. 洪水与风暴潮共同作用下的溃堤洪水一维、二维耦合模型及应用[J]. 南水北调与水利科技, 2017, 15(5): 43-49. https://www.cnki.com.cn/Article/CJFDTOTAL-NSBD201705007.htm

    WANG X J, HU B, YUAN X M, et al. A 1D-2D coupled model of the dam-break flood under the joint action of flood and storm surge and its application[J]. South-to-North Water Transfers and Water Science & Technology, 2017, 15(5): 43-49. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-NSBD201705007.htm

    [13] 程涛, 何万成, 赵海镜. 基于一二维耦合技术的溃坝洪水模拟研究[J]. 水力发电, 2022, 48(3): 21-23, 81. https://www.cnki.com.cn/Article/CJFDTOTAL-SLFD202203005.htm

    CHENG T, HE W C, ZHAO H J. Research on dam-break flood risk analysis based on 1-D and 2-D hydraulic coupling calculation model[J]. Water Power, 2022, 48(3): 21-23, 81. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SLFD202203005.htm

    [14]

    WU J Y, CHENG Y G, ZHANG C Z, et al. Three-dimensional simulation of balloon dynamics by the immersed boundary method coupled to the multiple-relaxation-time Lattice Boltzmann method[J]. Communications in Computational Physics, 2015, 17: 1271-1300.

    [15]

    KÖRNER C, THIES M, HOFMANN T, et al. Lattice Boltzmann model for free surface flow for modeling foaming[J]. Journal of Statistical Physics, 2005, 121(1/2): 179-196.

    [16]

    ZHOU J G, CAUSON D M, MINGHAM C G, et al. The surface gradient method for the treatment of source terms in the shallow-water equations[J]. Journal of Computational Physics, 2001, 168(1): 1-25.

    [17]

    LIU H F, ZHOU J G, BURROWS R. Lattice Boltzmann simulations of the transient shallow water flows[J]. Advances in Water Resources, 2010, 33(4): 387-396.

图(11)
计量
  • 文章访问数:  278
  • HTML全文浏览量:  52
  • PDF下载量:  64
  • 被引次数: 0
出版历程
  • 收稿日期:  2022-06-23
  • 网络出版日期:  2022-09-27
  • 刊出日期:  2022-11-29

目录

/

返回文章
返回