1. 引言
在山溪型河道中，部分现存桥梁行洪能力不足，为确保行洪安全，需采取相应的防洪措施。通过构建水动力模型，利用洪痕水位对模型进行校验，可有效模拟河道及卡口桥梁的综合阻水效果。一维水动力模型较常用于山区河道的水面线计算，同时对于规则的平板桥，可以通过涵洞加堰流方式进行概化，但对于造型复杂的桥梁，一维模型难以准确地概化其阻水效果。三维模型可精确模拟水流经过复杂构筑物的水流形态，然而对于较长河道来说计算资源消耗巨大。因此考虑将一维模型和三维模型相结合，利用一维模型可高效地计算长河道水面线，再通过三维模型对局部复杂构筑物进行精细概化，能兼顾计算效率和精度。
2. 研究方法
2.1. 河流水动力模型
为得到满足城镇防护标准要求的河道水面线，可采用基于圣维南方程的一、二维水动力模型，常用软件有MIKE11、MIKE21、HEC-RAS等。Szydłowski [1] 采用有限体积法的显式格式，求解了描述非定常自由表面流的沿水深平均的N-S方程即二维浅水方程，能较准确地得到河道中的桥墩附近水流结构的详细信息。Ghostine [2] 等分别采用一维和二维模型对河道干支汇流进行研究，一维模型在干支流汇流处使用的方程往往依赖于实验结果给出的经验参数，故适用范围有限；二维模型计算结果相比一维模型，在22.5˚、45˚和90˚三种汇流角工况下均与实验吻合度更好。王圆圆 [3] 等采用MIKE21模型模拟了施工导流围堰工程，并通过物理模型试验进行验证，给出施工导流围堰的优化方案。Bhallamudi [4] 等采用MacCormack格式对变换坐标系中求解二维水动力方程，可较为准确的计算河道超临界和次临界混合流。Ladopoulos [5] 利用现代奇异积分方程方法，再通过有限差分法确定了特殊明渠过渡段的自由表面轮廓及水流的水力特性。二维计算常用软件有MIKE21 [6] 、SMS等。一般而言对于展向和纵向远大于垂向尺度的洪水漫流计算，采用二维浅水波近似方法已经足够满足水流特性分析需要。然而，对于垂向变化较大、弯道存在螺旋二次流等情况下，例如因溃坝、山洪爆发发生的洪水流经弯道或者复杂的桥梁构筑物，则需要用到变量参数更为完备的三维模型。何勇 [7] 等构建了弯曲型交汇河口三维水动力模型，模拟交汇前后的水位变化以及不同位置流速垂向变化。Paul [8] 等采用三维模型模拟弧顶溢洪道的流量，计算了不同P/Hd(溢洪道高度与设计水头的比值)工况下的流量，结果表明数值模型与物理模型结果趋势上比较一致，数值模型的精度随P/Hd值增大而提高。
2.2. 桥梁壅水计算方法
在一维水动力模型中，如果跨河桥梁与河道断面相比，存在水流形态剧烈收缩和扩张，那么利用一维模型精确预测桥梁上游收缩段和下游扩张段的能量损失变得较为困难。John [9] 利用二维模型对一维模型中桥梁参数进行比对分析，提高了桥梁壅水计算精度，不过局限于普通垂向一致的桥梁。在MIKE11中，对于桥梁壅水计算推荐采用涵洞加堰的方式，其中涵洞可以模拟任何形状的桥孔，并可适用于多种水流条件，包括完全淹没流、部分淹没流、临界流、孔流、自由出流等。一般而言小流量不会造成明显壅水，只有当流量到达一定量级以上，才会因为桥梁阻水导致上游水位显著壅高。
针对山溪河道特性，采用率定过的一维模型能获取较为准确的洪水位、流速等水力要素，与二维模型相比一维模型建立更为便捷，计算速度更快。对于山区城镇防洪规划设计而言，河宽往往比河长小的多，主要关注河道的洪水位、流速等水力要素，因此多采用一维模型进行分析计算。用于描述一维河道模型的圣维南方程，有连续性方程和动量方程两个基本方程组成，可以认为是质量守恒定律和牛顿第二定律的一种微分表达形式：
(1)
式中：Q为流量(m3/s)；A为断面面积(m²)；q为旁侧入流单宽流量(m2/s)；α为动量校正系数；g为重力加速度(m/s2)；h为水位(m)；C为谢才系数；R为水力半径(m)。
利用Abbott六点隐式差分格式求解圣维南方程，求解时，将河道离散成水位、流量相间的计算点。将桥梁壅水公式中计算出的Δh转化为圣维南公式中动能方程组的追赶系数的改变，最后求解线性差分方程组。
当河道中存在较为复杂的空腹拱形桥梁，一维模型则无法模拟桥梁附近河道断面的流速分布情况，即使用二维模型也无法模拟桥梁附近水域的垂向流速变化，因此为得到桥梁上下游更加完备的流场空间分布信息，采用FLOW-3D三维水动力模型软件 [10] ，可以有效模拟各拱洞出流及相互距离对河道壅水的影响。FLOW-3D通过VOF(自由流体表面跟踪算法)方法追踪液–液或液–固交界面 [11] ，该方法跟踪流体体积的变化而不是追踪自由表面的粒子运动，因此可以提高边界条件和界面跟踪的精确度。模型基于笛卡尔坐标系，对完整NS方程的数值求解，可用于水利工程、水工建筑、水环境领域的流体三维仿真。其主要控制方程包括连续性方程、动量方程，湍流模型采用RNG k-ε模型 [7] ，自由水面则通过流体体积函数来确定。
连续性方程为
(2)
动量方程为
(3)
式中：ui为速度分量，Ai为可流动面积分数分量，RDIF为湍流耗散项，RSOR为质量源项，VF为流体体积分数，ρ为流体密度，p为压力，Gi为重力加速度分量，fi为黏性加速度分量。
2.3. Froude数
Froude数是水的惯性力与重力之比，计算公式为Fr = U2/gL，L表示特征长度，Fr > 1为急流，Fr < 1为缓流，计算Fr数选取L有水深型和面积型 [12] ，针对山区河流，L可用平均水深来表示。通过Fr数计算，可以分析河道急缓流变化，进而剖析造成桥梁壅水的成因。
3. 工程实例
3.1. 研究区域概况
仕阳溪属交溪水系，位于浙南山区浙闽交界处泰顺县境内，主源发源于九峰乡海拔1237.3 m九峰尖，在泗溪镇处与南溪、西溪、北溪汇合后称东溪，向西南至雪溪乡纳三魁溪，于南坑纳雪临溪进入仕阳镇境内后合称仕阳溪。仕阳溪干流石龙桥以上流域面积377.3 km2，主河道长42.04 km，河道坡降6.77‰。石龙桥为多孔拱桥，建于河道较窄处，对洪水位抬升影响较大。该桥有一个主拱，在左右拱肩上共有7个小拱桥洞，主拱半径242.6 m，跨度59.5 m，小拱桥洞洞宽为5 m，见图1。
Figure 1. Overview map of watershed aboved Shilong bridge location
图1. 石龙桥位置示意及以上流域概化图
3.2. 模型构建和验证
3.2.1. 一维模型构建和验证
一维模型采用MIKE11模型软件，其计算准确往往依赖于经验确定糙率参数，通过历史洪水的率定，可大大提高模型参数的可靠度。2005年泰利台风期间，仕阳镇发生一场较大洪水，见图2，最大淹没水深达1.7 m，通过洪水调查获得沿线主要控制断面洪痕。同时根据同期实测降雨换算为短历时暴雨参数，并推算洪水流量为3460 m3/s，作为模型边界条件，对仕阳镇主要河段进行参数率定。
Figure 2. Schematic map of Shiyang Town’s 20-year design flood covered range
图2. 仕阳镇20年一遇洪水淹没范围示意图
为模拟桥梁孔洞的过水能力，采用上文2.2提及的涵洞加堰的形式进行概化，并通过调整断面综合糙率和涵洞概化桥梁的水头损失系数，对验证洪水位成果进行逐步逼近，最终获得模拟水位与洪水调查水位(简称洪调水位)差值在0.02 m以内，见表1。可见采用涵洞概化桥梁的方式计算河道水面线阻水效果是合适的。
Table 1. Contrast table of simulated and investigated flood water level
表1. 模拟水位与洪调水位对比表
由于一维模型对拱桥局部壅水模拟存在困难，三维模型可以对局部构筑物进行精确模拟，综合考虑一维模型和三维模型的特点，采用一维模型模拟全程河道水面线，并将桥下水位摘取，作为三维桥梁壅水模型的边界，见图3。
Figure 3. Schematic diagram of hydraulic calculation sections of Shiyang Town Reach
图3. 仕阳镇河段水利计算断面示意图
3.2.2. 三维模型构建
三维模型包括河道、岸上地形及桥梁模型，河道、岸上地形通过地形图抽取高程点，导入FLOW-3D可生成stl格式文件。桥梁模型利用cad三维建模软件在横截面上进行拉伸，同样可以导出stl格式文件。将桥梁旋转平移至合适的位置与地形文件叠加，构成三维模型计算区域，见图4。在FLOW-3D里，对于河道明渠流模型的边界设置，上游采用流量边界(volume flow rate)并采用自然入口算法来评估水位，即每个边界单元的水位与直接下游内部单元的水位相同，边界流类别选亚临界流。
Figure 4. Three-dimensional model of river near bridge
图4. 桥梁附近河道三维模型示意图
下边界采用自由出流边界(outflow)：随着洪水演进，末断面水位陡降，流速10 m/s以上，不符合实际水位。
下边界采用水位边界：计算至200多步时，步长小于限定值，自动终止计算，计算不稳定，主要原因是固定水位边界在大流量下可能存在壁面反射波，不利于计算稳定。
下边界采用水位流量关系：由一维模型提取控制断面的水位流量关系，作为三维模型的下边界。避免反射波，可以达到计算稳定。
3.2.3. 模拟范围、网格设置
一维模型模拟范围，在仕阳镇区域范围，沿河道布设断面14个，共计2.7 km，平均间距209 m。在石龙桥处加设涵洞与堰组合水工构筑物。
三维模型计算区域限定为长500 m、宽230 m、高19 m的长方体，布设200 × 120 × 20的网格，共计48万个网格，在纵向桥梁上下游、横向有桥孔位置及垂向水气交界处进行加密，布设150 × 80 × 16网格，见图5。
为分析网格尺寸对计算精度的影响，对以上网格进行加密，并比较计算成果。将长、宽、高三个方向网格数分别增加一倍，布设400 × 240 × 40的网格，共计384万个网格。加密前后网格生产的桥梁三维表面形态如下图所示。加密前桥梁孔间隔板存在镂空，加密后的桥孔隔板密实完整，整个桥梁的三维表面更加光滑，见图6。
经计算，沿程水位结果在网格加密前后较为接近，见图7左图；沿程各测点位置的Fr数在网格加密前后也较为接近，见图7右图。以20年一遇工况为例，计算最终沿程水位标准偏差为0.066，若假定加密后计算成果为
，加密前成果为
，采用
，计算纳什效率系数计算得0.97。
Figure 5. Mesh generation of local bridge 3D model
图5. 局部桥梁三维模型网格划分
Figure 6. The 3D surface form of the bridge with mesh refinement (left and right of figure are mesh before and after refinement)
图6. 网格加密前后桥梁三维表面形态(图左右分别为网格加密前后)
Figure 7. Comparison of water level and Fr number with different grid number
图7. 加密网格前后沿程水位和Fr数比较分析图
3.2.4. 计算条件与参数设置
上游给定流量边界，下游给定水位流量关系曲线。仕阳镇石龙桥控制断面流域面积大于50 km2，可采用瞬时单位线法计算设计洪水。通过计算得到不同重现期下设计洪峰流量，其中20年一遇设计洪水洪峰流量为3370 m3/s。
3.3. 计算方案
3.3.1. 不同流量级，壅水高差
选取洪水重现期5年、10年、20年和50年对应的设计洪水流量为1960 m3/s、2710 m3/s、3370 m3/s、4160 m3/s，对应下游洪水位290.78 m、291.98 m、292.88 m、293.79 m。分别计算各重现期洪水下桥梁的壅水，随着重现期加大桥梁处流量加大，桥梁壅水高度也越大，见图8左图。对于不同重现期洪水，桥梁下游的Fr数随着重现期增大而增大，整体上桥梁下游的Fr数大于桥梁上游的Fr数，见图8右图。而桥梁上游的Fr数随着流量变化先增大后减小，主要由于水位上涨引起桥孔阻水，进一步抬升水位，水深加大而流速无明显加大，造成Fr数减小。
Figure 8. Comparison of water level and Fr number with different return periods
图8. 不同重现期沿程水位和Fr数比较分析图
3.3.2 桥梁下游湍流强度分布
在桥梁下游25 m处，河道横断面在两岸各形成二次流漩涡，尾流处湍流强度局部加大，比较不同流量级别，在较大流量下的尾流处的湍流强度较大(见图9)，将造成的更大的水流能量损失，流速降低导致水位抬升。
Q = 4160 m³/s
Q = 3370 m³/s
Figure 9. Turbulence intensity distribution at bridge and downstream with different discharge
图9. 不同流量下桥梁处及下游湍流强度分布图
3.3.3. 河道开槽行洪分析
经洪水调查及模型仿真分析，石龙桥及下游河道不满足行洪标准，根据石龙桥洪水影响评价报告措施，对阻水河道新开43 m × 4 m河槽，槽低高程为277 m，可有效提高仕阳溪石龙桥断面行洪能力。开槽前后的三维模型见图10，计算20年一遇洪水流量下不同下边界水位条件对桥梁的壅水。
对于不同下游边界，分别计算开河槽前后桥梁附近水面线见图11，当下游边界水位为高水位时，河道开槽对桥梁壅水的消减作用不明显，桥梁上游水位仅下降0.3 m，Fr数最大值由0.46降至0.36。当下游边界水位为低水位时，河道开槽对桥梁壅水的消减作用非常显著，桥梁上游水位下降1.2 m，Fr数最大值由0.86降至0.59。同时不论是开槽前还是开槽后，下边界为高水位的桥梁壅水均小于下边界为低水位的桥梁壅水。
开槽前 开槽后
Figure 10. 3D model before and after dredging
图10. 新开河槽前后三维模型
Figure 11. Comparison of water level and Fr number with different downstream water level
图11. 开槽前后不同下边界水位沿程水位和Fr数比较分析图
4. 结果和讨论
4.1. 不同流量和Fr数壅水结果比较
由图12可知，不同重现期洪水流量下，桥梁壅水与流量呈正相关，相关系数R2值为0.96；桥梁壅水与Fr值也呈正相关，相关系数R2值为0.9548。
在河道中心线上，Fr最大值出现在桥梁下游处，洪水流量、Fr数及壅水高度三个参数呈显著正相关。
4.2. 三维与一维模型结合拓展了复杂形态桥梁壅水成果
像多孔拱桥这类过水断面较为复杂的桥梁可以通过三维模型壅水计算获得不同流量下的桥梁阻水特性，进而可对一维模型构筑物的水头损失系统进行修正，可为防洪规划设计、山区洪水预警预报提供技术支撑。同时，当采取河道开槽、桥梁拓建等措施时，也可以通过三维模型直观地模拟局部断面的行洪能力，与一维模型相比不依赖于经验参数，可为防洪工程方案比选提供量化依据。
Figure 12. Correlation analysis diagram of the cause of bridge backwater
图12. 桥梁壅水成因相关分析图
5. 结论
本文研究山区多孔拱桥在不同重现期洪峰流量下的阻水特性，提出了一维河道模型和局部三维桥梁阻水精细模型结合的计算方法，通过有效网格数量判别，提高计算效率，同时保证计算精度。计算结果表明：在不同流量下桥梁壅水与Fr数有很强的相关性。下游边界水位较低时，疏浚河床可以有效降低Fr数，进而减少桥梁的壅水，降低整体河道水面线；但是如果下游边界水位较高，疏浚河床对降低Fr数较小，桥梁壅水减小不显著，降低整体河道水位作用也不明显。
本研究中在不同重现期场景下，进行的桥梁壅水模拟仅考虑现状为特定形态的桥梁，在今后研究中，可开展桥梁形态流线型优化设计，减少桥下尾流的湍流强度，以达到较小阻力的目的；此外本研究采用恒定流边界，而实际洪水过程为非恒定流，在今后研究中可采用实际洪水过程对桥梁壅水进行非恒定流模拟分析，从而更加符合天然河道的槽蓄情况。