布置于明渠上的水电站在发生增、甩负荷等瞬变过程时往往会在上游水库引发较大涌浪，产生机组出力和引航道水位剧烈波动等问题，对水电站的运行安全产生较大危害^[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]}实现。

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

1.1   一维明渠圣维南方程组及求解方法

1.1.1   控制方程

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

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

1.1.2   数值解法

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

1.2   三维自由液面N-S方程组及求解方法

1.2.1   控制方程

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

式中：u、p、f为分别流速、压强、外力；ν为流体的运动黏度；ρ为流体密度。

1.2.2   数值解法

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

式中：e_α为离散速度方向矢量。

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

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

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

式中：ω_α为权系数，取值可参考文献[14]。

1.2.3   三维自由液面处理方法

为了追踪自由液面，参考文献[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)的物理意义，可以利用分布函数在相邻网格点间的迁移实现质量交换。基于上述思想，对于t+Δt时刻位于x位置的网格点，相邻网格点既可能是流体节点，也可能是自由液面节点，当其与相邻流体节点进行质量交换时，沿x方向的质量变化(Δm_α(x，t+Δt))可表示为

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

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

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

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

1.3   一维-三维计算方法的耦合

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

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

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

下载: 全尺寸图片 幻灯片

式中：Q_i为流入和流出汊点的净流量(流入为正，流出为负)；Z_i为与汊点相连的第i个河段的水位，i=1, 2，…，M。

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

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

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

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

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

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

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

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

2   方法验证

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

2.1   潮汐波浪运动算例

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

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

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

本例中，计算空间步长Δ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

下载: 全尺寸图片 幻灯片

2.2   局部溃坝波算例

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

图  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

下载: 全尺寸图片 幻灯片

2.3   顺直明渠瞬变流算例

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

下载: 全尺寸图片 幻灯片

3   葛洲坝水库涌浪特性模拟

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

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

3.1   瞬变流计算条件

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

下载: 全尺寸图片 幻灯片

3.2   大江机组甩负荷涌浪模拟结果

各时刻葛洲坝水库水位和水平流速计算结果如图 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基本一致。可见，计算结果基本符合物理规律，说明一维-三维耦合的河库瞬变流模拟方法基本具备模拟实际明渠瞬变流的能力。

4   结论

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

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

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

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

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

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