大坝泄水时大量空气被卷吸进入泄洪水流并进入下游消能池内，气泡由于消力池中的高承压作用迅速溶解，从而达到相对于大气压的总溶解气体(total dissolved gas, TDG)过饱和状态。在某些水深较大、紊动较弱的水体中，特别是在大型深水库区，过饱和TDG的释放是一个很缓慢的过程。在此过程中可能对水生生物尤其是鱼类造成损伤，导致鱼类患“气泡病”(gas bubble disease)甚至死亡^[1]。

目前已有大量室内实验和原型观测的研究结果表明，高坝泄洪生成的过饱和TDG在下游河道中的释放过程受诸多因素的影响，主要与水深、紊动强度、水温、风速、含沙量以及河流形态等因素相关^[2-4]。而目前针对过饱和TDG的数值模拟研究主要基于Eulerian法。美国陆兵军工团根据对美国Columbia河上TDG饱和度的原型观测结果，首次提出纵向一维恒定流预测模型^[5-6]，模型中将过饱和TDG的沿程释放的影响因素归结为释放系数。之后文[7-8]在原型观测和机理实验的基础上，分别提出了综合考虑分子扩散和紊动扩散的影响、考虑河流底坡和紊动动能的影响的过饱和TDG释放系数估算模型。随着梯级电站的建设，大量学者也开始对梯级电站库区内的过饱和TDG释放开展研究。冯镜洁^[8]分别考虑了水体内部过饱和TDG释放和水气界面的传质作用，建立了适用于大型深水库的宽度平均立面二维过饱和TDG输移释放模型，该模型在实际工程中得到了很好的验证；在此基础上，马倩^[9]进一步改善了模型，在模型中加入了温度方程，考虑了水温对过饱和TDG输移释放的影响，计算结果与原观数据拟合较好。Witt等^[10]发展了降阶过饱和TDG模型用于预测美国哥伦比亚河中部流域7个水电站尾水渠的过饱和TDG分布。目前国内外关于过饱和TDG释放过程的数值模拟研究主要集中在一维和二维计算，但河道内的地势往往弯曲复杂，单靠一维或二维模型的计算很难完全模拟出整个库区或长范围河道内的TDG分布变化。

目前已有研究表明，TDG过饱和水体中鱼类患气泡病的严重程度与TDG过饱和程度、鱼苗种类、鱼的生长发育阶段、鱼苗体长、暴露时间、水体水深、水体温度等因素相关。目前在工程应用上对这方面的研究主要集中在TDG过饱和度和鱼类的暴露时间。当前还没有对鱼类可适应的TDG饱和度的统一标准，美国环境保护部门(US Environmental Protection Agency)提出河流中的TDG饱和度应限制在110%以下，黄翔等^[11]通过室内实验测试了岩原鲤鱼苗在不同TDG饱和度水箱中的半数死亡时间。结果显示，随着TDG饱和度的增加，岩原鲤的半数死亡时间缩短，并据此研究成果提出115%为适用于岩原鲤的TDG过饱和耐受性阈值。对于暴露时间，Weitkamp等^[12]将实验鱼置于TDG饱和度为120%~126%的天然水体表面8、12、16 h后，实验鱼的死亡率分别为1%、4%、12%；四川大学针对我国西南地区部分鱼类的TDG过饱和耐受性进行了系列实验研究^[13-14]，通过研究发现，导致实验鱼半数死亡的TDG饱和度随持续暴露时间的延长而降低。目前的研究手段还未能完全模拟鱼类在不同TDG时空分布下的运动轨迹与暴露时间，为研究在过饱和TDG水体中的暴露时间对鱼类的影响，本文采用TDG的持续滞留时间为参考因子以反映鱼类在TDG水体中的暴露时间。采用Eulerian法对过饱和TDG进行模拟时，模拟结果往往很难得到TDG的持续滞留时间，而Lagrangian法则可以很好地通过对粒子运动的模拟，计算得到粒子在某个区域中的运动轨迹和滞留时间等参数，以此反映计算区域内的TDG持续滞留时间。为了更好地分析过饱和TDG对鱼类的影响，有必要将Eulerian法与Lagrangian法结合起来。目前，Eulerian-Lagrangian模型已在气泡运动、悬浮颗粒流等领域得到了很好的应用^[15]，但尚未见Eulerian-Lagrangian方法应用于对过饱和TDG对鱼类影响的研究。

综合以上分析，本文建立了Eulerian-Lagrangian模型分析大渡河上铜街子电站库区的过饱和TDG对鱼类的影响。大渡河是长江流域岷江的最大支流，水利资源丰富，目前已形成以下尔呷、双江口、猴子岩、长河坝、大岗山、瀑布沟等水电站为主要梯级格局的3库22级开发方案。大渡河流域内存在多种特有鱼类，大渡河梯级电站的运行对鱼类产生了不利的影响，为更好地保护鱼类资源，有必要对大渡河梯级电站中的过饱和TDG及其对鱼类的影响开展研究。本文采用Eulerian模型对大渡河铜街子电站库区过饱和TDG的释放展开三维数值模拟，并通过Lagrangian模型释放粒子以对粒子所受的TDG饱和度和持续暴露时间进行追踪，从而基于鱼类耐受性试验成果，分析过饱和TDG的影响。

1 计算模型

Eulerian-Lagrangian法用基于Eulerian守恒方程来描述流体流动，用Lagrangian方程来描述粒子在流场中的运动轨迹。由于Lagrangian法考虑了流体的瞬时流动特性对粒子的影响，更适于模拟粒子随时间的变化情况^[15]。为此，本研究所用的Eulerian-Lagrangian模型主要分为两部分，Eulerian法用于对流场和TDG饱和度的时空分布进行模拟，Lagrangian法用于对计算区域内的粒子进行模拟计算，得到粒子在前述Eulerian模型场内的运动轨迹、运动速度和粒子的TDG饱和度等。

1.1 TDG输移释放模型

TDG输移释放采用Eulerian法，包括水动力学方程和TDG输运方程。

1.1.1 水动力学方程

水动力学方程采用了钢盖假定的单相流模型进行求解。模型包括连续方程和动量方程。

连续方程：

$ \nabla \cdot \boldsymbol{U}=\bf{0}. $

(1)

动量方程：

$ \frac{\partial \boldsymbol{U}}{\partial t}+\nabla \cdot(\boldsymbol{U} \otimes \boldsymbol{U})=-\nabla \frac{p}{\rho}+\nabla \cdot(\mu \nabla \boldsymbol{U}) . $

(2)

其中：U为流速(m/s)；p为压强(Pa)；ρ为密度(kg/m³)；μ为运动黏度(m²/s)；t为时间(s)。

采用RANS紊流模型进行求解，RANS紊流模型的核心在于求解时均化的Reynolds方程，既可以节约计算量，同时在工程应用上也有很好的效果。本研究的模型采用Realizable k-ε紊流模型，相关方程如下：

$ \frac{\partial(\rho k)}{\partial t}+\frac{\partial\left(\rho k u_{i}\right)}{\partial x_{i}}=\frac{\partial}{\partial x^{j}}\left[\left(\mu_{\mathrm{m}}+\frac{\mu_{\mathrm{t}}}{\sigma_{\mathrm{k}}}\right) \frac{\partial k}{\partial x_{j}}\right]+G_{k}-\rho \boldsymbol{\varepsilon}, $

(3)

$ \begin{array}{c} \frac{\partial(\rho \varepsilon)}{\partial t}+\frac{\partial\left(\rho \varepsilon u_{i}\right)}{\partial x_{i}}=\frac{\partial}{\partial x^{j}}\left[\left(\mu_{\mathrm{m}}+\frac{\mu_{\mathrm{t}}}{\sigma_{\varepsilon}}\right) \frac{\partial \varepsilon}{\partial x_{j}}\right]+ \\ \rho \mathrm{C}_{1} E \varepsilon-\rho \mathrm{C}_{2} \frac{\varepsilon^{2}}{k+\sqrt{v \varepsilon}}. \end{array} $

(4)

其中：k为流体紊动动能(m²·s^-2)；ε为流体紊动动能耗散率(m²·s^-3)；σ_k为紊动能所对应的Prandtl数；σ_ε为紊动能耗散率所对应的Prandtl数；v为运动黏度(m²/s)；μ_m为分子黏性系数(N·s·m^-2)；μ_t为紊流黏性系数(N·s·m^-2)；G_k为由平均速度梯度引起的紊动动能产生项。

1.1.2 TDG输运方程

对过饱和TDG的释放过程概化为两个部分：在水体内部，以源项模拟过饱和TDG在水体内部的释放过程；在自由水面处，将过饱和TDG向大气的扩散过程概化为水气界面的传质过程。其中，TDG的内部释放过程由TDG浓度对流扩散方程表示：

$ \frac{\partial G}{\partial t}+\nabla \cdot(\boldsymbol{U} G)=\nabla \cdot\left(\left(\nu+\frac{\nu_{t}}{S_{c}}\right) \nabla G\right)+S_{\mathrm{inner}}, $

(5)

$ {S_{{\rm{inner}}}} = {k_{{\rm{inner}}}}({G_{{\rm{eq}}}} - G). $

(6)

其中：G为总溶解气体饱和度(%)：Sc为Schmidt数；ν为分子扩散系数(m²/s)；ν_t为紊动扩散系数(m²/s)；S_inner为水体内部TDG释放源项；G_eq为当地平衡饱和度(%)；k_inner为水体内部TDG释放系数(s^-1)，由实验或原型观测率得到。

过饱和TDG的释放还需要考虑水气界面传质，在自由水面上，过饱和TDG以分子扩散的形式通过水气界面传质进入大气。因此，本文模型在计算区域的水体表面边界添加通量方程以考虑过饱和TDG的自由界面传质部分。

$ F_{\mathrm{w}} =k_{\text {surface }}\left(G_{\text {patm }}-G\right), $

(7)

$ k_{\text {surface }} =0.085 k^{1 / 2}+0.001\;4 . $

(8)

其中：F_w为过饱和TDG在水面边界处的通量；G_patm为当地大气压下的TDG饱和度；传质系数k_surface的计算采用李然等^[16]提出的水气传质系数与表面紊动动能的定量关系。

1.2 粒子运动模型

粒子运动模型基于离散单元法(discrete element method)，模型中不考虑粒子与粒子、粒子与壁面之间的碰撞作用。离散单元法根据Newton第二定律追踪单个粒子的运动，粒子i运动方程表示如下：

$ m_{i} \frac{\mathrm{d} \boldsymbol{U}_{\mathrm{pi}}}{\mathrm{d} t}=m_{i} \boldsymbol{g}+\boldsymbol{F}_{i} , $

(9)

$ \boldsymbol{I}_{i} \frac{\mathrm{d} w_{i}}{\mathrm{~d} t}=\sum\limits_{j=1}^{n_{i}}\left(\boldsymbol{T}_{t i j}+\boldsymbol{T}_{n i j}\right). $

(10)

其中：m_i、U_pi、I_i和w_i分别为粒子i的质量、速度、惯性力和旋转速度；F_i为粒子受到的包括拖曳力、接触力等；T_tij和T_nij分别为正切力和滚动摩擦力。

粒子在各TDG饱和度下的滞留时间t_TDG由粒子在各TDG饱和度下的计算时间步长dt累积相加而得。

1.3 数值求解

本文基于OpenFOAM开源平台开展模拟求解。其中，TDG释放模型以pisoFoam求解器为基础，在计算流场、压力场的基础上添加TDG输运方程，以PISO算法反映压力变化对速度场的影响。对粒子运动的模拟计算以Lagrangian法的icoUncoupleKinematicPacelFoam求解器为基础，在计算工况的入流边界上释放100个粒子。在方程离散格式上，瞬时项采用Gauss线性格式；对流项采用GausslimitedLinearV格式；Laplace项使用Gauss线性校正插值格式。

2 模型验证

模型关于水动力学模拟的验证采用孙兴^[17]参考Laufer^[18]所做的平面槽道流动实验对OpenFOAM中的pisoFoam求解器所模拟出的流场进行验证，计算采用了RANS紊流模型。通过计算值与实验值的对比分析可知，采用pisoFoam求解器计算得到的数值模拟结果与实验值具有良好的一致性，相关参数与实验结果的误差可以控制在6%之内，说明该求解器对平面槽道流动的RANS模拟可以得到比较准确的数值结果。

对于模型中的TDG饱和度分布的验证，采用黄膺翰^[19]所做的室内顺直水槽实验作为TDG释放模型的验证工况。

2.1 室内实验介绍

实验在四川大学进行，水槽实验装置主体为有机玻璃水槽，槽身尺寸为15 m×0.5 m×0.3 m，坡降为0.45‰。在水槽上游端上方安装保证来流水体的TDG饱和度保持不变的装置。实验测量TDG饱和度、温度、流量、流速、流向、水深等参数。其中，TDG饱和度和温度由分别放置在水槽上下游的2台Point Four测定仪同时测量。流速测量采用由中国南京圣荣仪器设备有限公司生产的LGY-2型智能流速仪。实验设置5个不同流量的工况，测量工况及测量结果如表 1所示。

2.2 模拟结果验证

采用节2.1中的模型对实验水槽中的过饱和TDG水流的运动进行模拟研究。

2.2.1 网格划分及边界条件

计算区域包括整条水槽，如图 1所示，计算区域采用3D结构化网格，网格在3个方向的尺寸分别为0.05、0.03、0.007 m。计算区域内设有入流边界、出流边界、顶部边界和边壁4个边界条件，各边界条件如表 2所示。通过在计算模型的顶部边界上采用自编的codedFixedValue边界条件以在自由液面添加表示TDG自由表面传质通量的方程(7)。

本研究中使用的主要边界条件如下：

1) fixedValue：计算参数边界值为指定的固定值。

2) zeroGradient：计算参数边界的法向梯度为0。

3) slip：对于标量场使用零梯度条件，对于矢量场的切向分量为零梯度条件，法向分量为fixedValue。

4) codedFixedValue：自定义边界。

2.2.2 模型验证成果

以验证工况5为例，图 2为验证工况5在靠近水面处的TDG分布示意图，可以看出，TDG在横向上呈均匀分布，随着水流的流动，TDG在水槽内沿程释放，饱和度逐渐降低。各验证工况的沿程TDG平均值分布如图 3所示，可以看出，TDG饱和度的沿程分布呈现明显的线性关系，相关系数高达0.99。通过拟合曲线计算得到各个验证工况下的TDG沿程释放系数，计算结果如表 3所示。从表中数据看出，验证工况1的入口流速最小，过饱和TDG释放最慢，TDG沿程释放率为0.29 m^-1；随着入口流速的增大，各工况的TDG沿程释放率变大，验证工况5的TDG沿程释放率为0.52 m^-1。以上结果表明，随着水流流速的增加，水体的紊动强度相应增大，这将促进过饱和TDG的释放。

各验证工况的下游出口处TDG饱和度时均值如表 3所示，对比试验测量数据，计算值与实验测量值的最大差值仅为1.1%，在可接受差值范围内。计算结果说明本研究中所用的TDG释放模型可以用于对TDG释放的数值模拟计算中。

3 典型库区TDG分布及对鱼类影响模拟 3.1 计算工况选择

选取大渡河铜街子水电站库区为研究对象，库区上游为龚嘴水电站。铜街子电站正常蓄水位为474.0 m，总库容2.60亿m³，有效库容5 620万m³，水库最低运行水位为469.0 m，具有日调节性能。铜街子库区的常年洪水流量和多年平均流量分别为6 280和1 490 m³/s。上游龚嘴电站为混凝土重力坝，最大坝高85.6 m，电站设计水头48 m，最大水头53.08 m，最小水头34.7 m，水电站装机容量70万kW，保证出力17.9万kW，多年平均发电量34.2亿kW·h。工程以发电为主，兼顾漂木。龚嘴水电站主要泄水建筑物包括4个表孔溢流道和1个冲沙底孔，其中3个表孔位于主河槽7#~9#坝段，简称河床溢流坝，为经常启用的泄水孔口；另一表孔位于原导流明渠4#坝段，简称明渠溢流坝。溢流坝孔口宽度为12 m，闸墩厚6 m，单孔溢流坝坝段宽18 m。龚嘴电站、铜街子电站及铜街子库区如图 4所示。

四川大学于2009年8月3日—7日期间对龚嘴水电站泄水进行了过饱和TDG原型观测。其中，在龚嘴坝下发电尾水处布置1个测点；在坝下830 m处的龚电大桥断面布置3个测点。原型观测的观测因子包括过饱和TDG、溶解氧DO和水温，主要观测仪器为Point Four TGP测定仪。

分别选取原型观测期间所观测到的龚嘴电站坝下的最大和最小TDG饱和度工况进行模拟计算，所选工况的入库流量在铜街子库区常年洪水流量与多年平均流量范围中。入流TDG饱和度采用泄洪水流与发电尾水充分掺混后的TDG饱和度，并认为TDG饱和度在深度上也达到充分掺混。计算工况如表 4所示。根据原型观测结果计算得到铜街子库区的内部TDG释放系数的取值为0.16 h^-1。

3.2 边界条件与网格划分

选用龚嘴坝下至铜街子坝前，约28 km长的库区范围为计算对象，计算区域如图 4所示。计算区域采用3D结构化网格，沿水流流动方向上，网格长度约为20 m，纵向上的网格约为8 m，深度方向上的网格平均大小约为2 m，总网格数量约为1.1×10⁵。为便于计算，对地形做了概化处理，计算区域表层边界采用钢盖假定。计算模型边界条件如图 5所示。

3.3 TDG空间分布模拟结果

在采用稳定入流的情况下，由于TDG的自由表面释放和水体内部释放，在计算库区内，TDG饱和度不断降低，工况1在大约20 h后，出口处的TDG饱和度达到稳定，此时出口断面的平均TDG饱和度为113.5%；工况2在大约10 h后，出口处的TDG饱和度达到稳定，此时出口断面的平均TDG饱和度为123.5%。图 6和图 7分别为两组工况的库区计算工况在靠近自由表面处的沿程TDG饱和度分布示意图，从图中可以看出，TDG饱和度在库区内沿程逐渐降低，对比出流TDG饱和度，工况1降低了11.1%，工况2降低了16.3%。由于采用混合均匀的TDG饱和度为入流条件，因此在计算结果上，并没有明显反映出TDG饱和度在横向和垂向上的显著差异。对比两个工况的计算结果可以看出，由于工况2的流速较大，TDG在库区中的释放速度较快，工况2的TDG饱和度降低值高于工况1。

3.4 粒子滞留时间模拟结果

在计算工况的入流边界分别释放100个粒子。其中工况1有84个有效粒子，工况2有90个有效粒子。工况1与工况2的各个粒子的运动轨迹及其TDG饱和度分别如图 8和图 9所示。粒子在计算区域内跟随水流运动，所受的TDG饱和度在运动方向上逐渐降低。对不同时刻下的各粒子所受的TDG饱和度取平均值，如图 10所示，可以看出，随着粒子运动时间的增加，粒子所受的TDG饱和度不断降低，受地势影响，TDG降低速度先增大再变小，运动时间3~6 h内的下降速度最大。工况1中的粒子TDG饱和度从入口处的124.6%降低至113.7%，工况2中的粒子TDG饱和度从入口处的139.9%降低至123.3%。各粒子在库区的滞留时间和在各TDG饱和度水体中的持续暴露时间如图 11所示，可以看出，各粒子在库区内的运动时间的分布较均匀。对各工况下的粒子信息进行统计分析，工况1中，粒子在库区中的平均滞留时间为14.9 h，粒子在大于115%TDG饱和度水体中的平均持续暴露时间为7.9 h，在大于120% TDG饱和度水体中的平均持续暴露时间为3.5 h。工况2中，粒子在库区中的平均滞留时间为11.6 h，粒子在大于125%TDG饱和度水体中的平均持续暴露时间为6.1 h，在大于130% TDG饱和度水体中的平均持续暴露时间为4.6 h，在大于135% TDG饱和度水体中的平均持续暴露时间为2.8 h。

4 基于鱼类耐受性的TDG过饱和风险分析

根据大渡河鱼类资源现状调查结果，大渡河水系内以鲤科鱼类为主，有31种，占大渡河水系鱼类的60%，其他鱼类占40%。自然水体中经济鱼类主要有：齐口裂腹鱼、重口裂腹鱼、虎嘉鱼、鲫鱼、黑鱼、黄鳝等。西南大学于2016年对大渡河进行生态调查，认为调查区域内有8种长江上游特有鱼类，即长薄鳅、齐口裂腹鱼、重口裂腹鱼、青石爬鮡、黄石爬鮡、壮体鮡和山鳅^[20]。为此，本文采用已有关于齐口裂腹鱼的TDG过饱和急性致死试验结果作为鱼类影响判别依据，其中以不同TDG饱和度下半致死时间(LT₅₀)代表河道内鱼类受过饱和TDG影响的程度和范围。半致死时间是指在一定的TDG饱和度下，引起实验鱼苗死亡率达到50%时所需的时间^[21]。

参考齐口裂腹鱼的TDG过饱和急性致死时间的试验结果^[22]，根据前述Lagrangian模型中计算出的粒子滞留时间结果统计分析TDG的持续滞留时间，如表 5所示。为考虑最不利情况，统计分析中采用了各工况有效粒子暴露时间的最大值与平均值。考虑到鱼类本身对TDG过饱和水流具有规避能力，本文所探讨的为最不利情况下的鱼类所受到的过饱和TDG的影响研究。

从结果可以看出，当上游梯级电站以工况1的泄洪方案持续泄洪时，各个TDG饱和度段的持续滞留时间均低于3种鱼类的LT₅₀实验值，工况1泄洪方案对库区内的鱼类的影响较小，鱼类在整个库区段内相对比较安全。当上游梯级电站以工况2的泄洪方案持续泄洪时，对于齐口裂腹鱼，其在TDG饱和度小于135%内的区域游动时，粒子的最大持续暴露时间较长；在TDG饱和度大于120%内的区域游动时，粒子的平均持续暴露时间大于半致死时间，说明鱼类在此区域内受伤的概率较大。计算结果表明，采用工况2泄洪方案持续泄洪时，鱼类受过饱和TDG影响的风险较大。针对本文中的不同泄洪流量与不同来流TDG饱和度的两组工况的计算结果，认为可以通过采取分段泄洪的调度方式，减缓大坝下游的TDG饱和度及对鱼类的影响。

5 结论

为模拟梯级电站库区中过饱和TDG的释放过程及饱和度分布，同时分析鱼类在TDG过饱和水体中所受到的影响，本文建立了Eulerian法和Lagrangian法相结合的Eulerian-Lagrangian过饱和溶解气体模型，其中，采用Eulerian模型对大渡河铜街子电站库区过饱和TDG的输移释放展开三维数值模拟，并通过Lagrangian模型释放粒子以追踪粒子所受的TDG饱和度和持续暴露时间，进而基于鱼类耐受性实验成果，分析鱼类受过饱和TDG胁迫的致死风险。针对铜街子库区的模拟表明，当上游龚嘴电站泄洪流量1 463 m³/s时，铜街子库区全部为安全区域，库区内鱼类受过饱和TDG影响较小；当龚嘴电站泄洪流量2 575.8 m³/s时，鱼类在TDG饱和度大于120%内的区域游动时，TDG持续滞留时间大于半致死时间，此时铜街子库区内鱼类受过饱和TDG影响的风险较大。本文为定量研究过饱和TDG对鱼类影响提供了一种新方法，同时为泄洪方案的优化研究提供了新的参考依据。未来将进一步研究当上游电站非恒定泄洪情况下，TDG饱和度随时间的变化规律和鱼类受暴露时间变化，进而定量评价非恒定泄洪的过饱和TDG对鱼类的影响风险。