桥梁下部结构对水流存在束狭作用，促使河床发生相应冲淤变化。建桥引起的河床变化可分为压缩冲刷和局部冲刷，压缩冲刷是指桥墩压缩断面后发生的普遍冲刷，包括天然河床冲刷和桥墩压缩水流引起的冲刷；局部冲刷是指在水流漩涡作用下桥墩周围河床形成的冲刷坑。目前，大量研究集中在桥墩局部冲刷^[1-5]，对压缩冲刷研究较少^[6-8]。桥墩压缩冲刷使得桥址断面出现整体冲刷下切，对建筑物基础安全不利^[9]。就长江中下游而言，三峡水库蓄水后，下泄沙量急剧减小，加剧了下游河床冲刷，因此在过江桥梁设计中，桥墩压缩冲刷预测是关键技术问题之一。

以往学者通过模型试验或野外实测分析，得到桥墩压缩冲刷预测经验公式，其中较为著名的有：铁路及公路规范64-1式、64-1修正式和64-2简化式^[10-11]，美国HEC-18式^[12]，但此类公式只能反映某一水沙条件持续作用下压缩冲刷平衡值，且系数存在不确定性。近年来，国内外学者开始采用水沙数模预测桥墩压缩冲刷。刘有录等^[13]通过调整桥墩附近网格局部水头损失系数，模拟了百年一遇特大洪水期间黄河中游某铁路桥墩冲刷，床面冲淤趋势与试验基本吻合，但具体量值存在差异；Lai等^{[6-9, 14]}对长压缩段矩形水槽试验进行了验证，得出模型可较好地模拟沿程水位变化及压缩段冲刷分布，但未考虑上游泥沙过程及悬沙造床作用。对于长江下游而言，三峡水库蓄水拦沙使得下泄的水沙过程明显有别于天然情况。下游河床冲淤复杂，相比于过江隧道^[15-16]，桥墩压缩冲刷影响因素较多，因此要想较好地预测桥墩压缩冲刷，模型需具备以下几方面条件：①选取合理水沙过程；②模型参数需进行多次率定和验证，能复验河道年内和年际冲淤变化；③能较好反映建桥前后水流泥沙运动。

本文以新建连云港至镇江铁路世业洲大桥桥位方案为背景，对桥墩压缩冲刷模型的建立、参数验证、桥墩边界处理、水沙过程选取、冲刷深度预测等进行研究，为类似问题的研究提供参考。

1 工程概况

新建的连云港至镇江铁路全长249 km，北承新长铁路，中联宁启铁路，南与京沪高铁、沪宁城际铁路相接，是连接苏北、苏中、苏南地区的南北纵向铁路主通道。铁路拟在长江世业洲至五峰山河段处跨越长江，过江通道拟采用跨江桥梁，根据桥梁位置的不同，分为世业洲桥位、五峰山桥位方案。其中，世业洲桥位方案为：桥址位于世业洲中上部，采用悬索桥方案，左汊主跨为406 m，右汊主跨为910 m。平面、立面布置及桥墩结构布置见图 1~2。

世业洲河段为微弯分汊河型，江中有世业洲将河道分为左右两汊，见图 1。左汊为支汊，2012年12月实测分流比为38.8%，右汊为主汊。由于河道有良好的进出口节点控制及岸线守护作用，总体河势变化不大，但局部区域演变仍较为明显。主要表现为汊道分流格局变化及洲滩冲淤调整。

(1) 汊道分流格局变化主要是左汊分流比增加，右汊分流减小。20世纪70年代中期至90年代初期，左汊分流比缓慢增加，年均增幅为0.1%；1995年后，受到连续大洪水作用，增幅速度明显加大，其中1997—1999年，年均增大达2.8%；2000年以后增速趋缓。世业洲左汊分流比由20世纪70年代的20%增至40%左右。

(2) 洲滩冲淤调整体现为世业洲头低滩持续冲刷后退；左汊整体冲刷，局部位置冲深幅度4 m；世业洲右汊进口段左侧边滩冲刷、右侧深槽淤积，中下段则表现为世业洲右缘边滩淤长，深槽冲刷。桥址断面冲淤变化见图 3。

2 模型建立和验证 2.1 控制方程组

正交曲线坐标系下的平面二维水沙数模控制方程组如下：

(1) 水流连续方程：

$ \frac{\partial Z}{\partial t}+\frac{1}{C_{\xi} C_{\eta}}\left[\frac{\partial\left(C_{\eta} H u\right)}{\partial \xi}+\frac{\partial\left(C_{\xi} H v\right)}{\partial \eta}\right]=0 $

(1)

(2) 水流运动方程：

$ \begin{array}{l} \frac{{\partial (Hu)}}{{\partial t}} + \frac{1}{{{C_\xi }{C_\eta }}}\left[ {\frac{\partial }{{\partial \xi }}\left( {{C_\eta }Huu} \right) + \frac{\partial }{{\partial \eta }}\left( {{C_\xi }Hvu} \right) + Hvu\frac{{\partial {C_\xi }}}{{\partial \eta }} - H{v^2}\frac{{\partial {C_\eta }}}{{\partial \xi }}} \right] = - \frac{{gu/{u^2} + {v^2}}}{{{C^2}}} - \frac{{gH}}{{{C_\xi }}}\frac{{\partial Z}}{{\partial \xi }} + \\ \frac{1}{{{C_\xi }{C_\eta }}}\left[ {\frac{\partial }{{\partial \xi }}\left( {{C_\eta }H{\sigma _{\xi \xi }}} \right) + \frac{\partial }{{\partial \eta }}\left( {{C_\xi }H{\sigma _{n\xi }}} \right) + H{\sigma _{\xi \eta }}\frac{{\partial {C_\xi }}}{{\partial \eta }} - H{\sigma _{\eta \eta }}\frac{{\partial {C_\eta }}}{{\partial \xi }}} \right] \end{array} $

(2)

$ \begin{array}{l} \frac{{\partial (Hv)}}{{\partial t}} + \frac{1}{{{C_\xi }{C_\eta }}}\left[ {\frac{\partial }{{\partial \xi }}\left( {{C_\eta }Huv} \right) + \frac{\partial }{{\partial \eta }}\left( {{C_\xi }Hvv} \right) + Huv\frac{{\partial {C_\eta }}}{{\partial \xi }} - H{u^2}\frac{{\partial {C_\xi }}}{{\partial \eta }}} \right] = - \frac{{gv/{u^2} + {v^2}}}{{{C^2}}} - \frac{{gH}}{{{C_\eta }}}\frac{{\partial Z}}{{\partial \eta }} + \\ \frac{1}{{{C_\xi }{C_\eta }}}\left[ {\frac{\partial }{{\partial \xi }}\left( {{C_\eta }H{\sigma _{\xi \eta }}} \right) + \frac{\partial }{{\partial \eta }}\left( {{C_\xi }H{\sigma _{\eta \eta }}} \right) + H{\sigma _{n\xi }}\frac{{\partial {C_\eta }}}{{\partial \xi }} - H{\sigma _{\xi \xi }}\frac{{\partial {C_\xi }}}{{\partial \eta }}} \right] \end{array} $

(3)

式中：(ξ, η)为曲线坐标系坐标；(x, y)为物理坐标；C_ξ，C_η分别为ξ, η方向拉梅系数；u，v分别为ξ，η方向流速分量；Z为水位；t为时间；H为水深；C为谢才系数，$C=\frac{1}{n} H^{1 / 6}$，n为糙率；σ_ξξ，σ_ηη，σ_ξη，σ_ηξ为紊动切应力。

(3) 悬移质不平衡输移方程：

$ \begin{array}{l} \frac{{\partial \left( {H{S_i}} \right)}}{{\partial t}} + \frac{1}{{{C_\xi }{C_\eta }}}\left[ {\frac{\partial }{{\partial \xi }}\left( {{C_\eta }Hu{S_i}} \right) + \frac{\partial }{{\partial \eta }}\left( {{C_\xi }Hv{S_i}} \right)} \right] = \frac{1}{{{C_\xi }{C_\eta }}}\left[ {\frac{{{\varepsilon _\xi }}}{{{\sigma _s}}}\frac{{{C_\eta }}}{{{C_\xi }}}\frac{{{C_\eta }}}{{{C_\xi }}}H\frac{{\partial {S_i}}}{{\partial \xi }}} \right) + \\ \frac{\partial }{{\partial \eta }}\left( {\frac{{{\varepsilon _\eta }}}{{{\sigma _s}}}\frac{{{C_\xi }}}{{{C_\eta }}}H\frac{{\partial {S_i}}}{{\partial \eta }}} \right)] + {\alpha _i}{\omega _i}\left( {S_i^* - {S_i}} \right) \end{array} $

(4)

式中：α_i为泥沙的含沙量恢复饱和系数；ω_i为第i组泥沙沉速；S_i为第i组悬移质含沙量；S_i^*为第i组悬移质挟沙力，ε_ξ，ε_η分别为ξ，η方向泥沙扩散系数；σ_s取常数1.0。

(4) 推移质不平衡输移方程：

$ \frac{{\partial H{S_{{\rm{b}}L}}}}{{\partial t}} + \frac{1}{{{C_\xi }{C_\eta }}}\left[ {\frac{\partial }{{\partial \xi }}\left( {{C_\eta }Hu{S_{{\rm{b}}L}}} \right) + \frac{\partial }{{\partial \eta }}\left( {{C_\xi }Hv{S_{{\rm{b}}L}}} \right)} \right] = {\alpha _{{\rm{b}}L}}{\omega _{{\rm{b}}\mathit{L}}}\left( {S_{{\rm{b}}L}^* - {S_{{\rm{b}}L}}} \right) $

(5)

式中：S_bL^*为第L组推移质挟沙能力；S_bL为第L组推移质含沙浓度；α_bL为第L组推移质恢复饱和系数；ω_bL为第L组推移质的沉速。

(5) 河床变形方程：

$ 悬移质冲淤引起的河床变形:\;\;\;\;\;\;\;{\gamma _{0i}}\frac{{\partial {Z_i}}}{{\partial t}} = {\alpha _i}{\omega _i}\left( {{S_i} - S_i^*} \right) $

$ 推移质冲淤引起的河床变形:\;\;\;\;\;\;\;{\gamma _{0{\rm{b}}L}}\frac{{\partial {Z_{{\rm{b}}L}}}}{{\partial t}} = {\alpha _{{\rm{b}}\mathit{L}}}{\omega _{{\rm{b}}L}}\left( {{S_{{\rm{b}}L}} - S_{{\rm{b}}L}^*} \right) $

$ 河床总的冲淤厚度为:\;\;\;\;\;\;\;Z = \sum\limits_{i = 1}^M {{Z_i}} + \sum\limits_{L = 1}^N {{Z_L}} $

(6)

式中：M为悬移质分组数；N为推移质分组数；γ_0i为i组悬移质泥沙干重度；Z_i为第i组悬移质引起的河床冲淤厚度；γ_0bL为L组推移质泥沙干重度；Z_L为第L组推移质引起的河床冲淤厚度；Z为河床总冲淤厚度。

2.2 模型设置

为数值求解水流运动方程组(1)~(3)，采用积分守恒的控制体积法离散方程，并利用SIMPLER方程式求解耦合方程，反复迭代直到流场收敛；对悬移质输移方程(4)和推移质输移方程(5)在离散后采用欠松驰技术及逐行扫描的TDMA技术隐式求解；对河床变形方程(6)采用有限差分离散，显式求解。

水流运动方程主要参数有：糙率n和紊动黏性系数ν_t。糙率反映了天然河流的阻力，可通过与实测水位及流速验证确定；紊动黏性系数${\nu _{\rm{t}}} = k{u_*}H$, u_*为摩阻流速，k为常数。

泥沙运动方程主要参数有：悬移质分组挟沙能力S_i^*、推移质分组输沙率g_bL^*，其中悬移质分组挟沙能力S_i^*=P_*S_*，P_*为分组挟沙力级配(P_*可根据文献[15]附录B提供的方法确定)，S_*为水流总的挟沙力，一般采用张瑞瑾公式确定：

$ S_{*}=k\left(\frac{\sqrt{u^{2}+v^{2}}}{g h \omega}\right)^{m} $

(7)

式中：ω为泥沙平均沉速；k，m为系数。推移质分组输沙率$g_{{\rm{b}}L}^* = {P_{{\rm{b}}L}}{\eta _L}g_{\rm{b}}^*$，P_bL为第L组泥沙在床沙中的百分比，η_L为隐蔽系数，g_b^*为推移质输沙率，可根据文献[15]附录A列出的公式计算确定。

2.3 计算中对桥墩的考虑

以往计算中主要采用局部地形修正、附加糙率等方法概化桥墩，但此类方法仅能模拟桥墩阻水效应，无法给出桥墩附近细部流场变化，因而也无法计算桥墩压缩冲刷。随着计算机性能和数值方法的提高，可采用加密网格法处理桥墩，使得网格与桥墩尺度一致。如何处理桥墩边界是准确模拟桥墩附近水沙运动的关键，根据模拟空间尺度不同，可分为小尺度和大尺度两类模型。

小尺度模型通常用于实验室桥墩模拟，网格尺度较小，为厘米级，桥墩边界可采用大系数法与壁面率相结合的方法进行处理，即桥墩内的网格(图 5中网格节点1~16所围的内部网格)采用大系数法处理，亦即在计算时，令其源项为一大系数，则其计算结果为0；而与桥墩壁面相邻的网格节点(网格节点1~13)则采用壁面率进行处理。

大尺度模型通常用于天然河道中桥墩模拟。河道边界复杂，网格尺度一般为米级别，采用壁面率时，壁面第一个计算节点与壁面间的距离为$z^{+}=z_{1} u_{*} / v$，30 < z⁺ < 100；对于长江下游而言，当流速取为2 m/s，水深取10 m，糙率为0.025时，z₁为0.001 m。可见，计算所需尺度远远小于桥墩尺度，故实际计算时，很少采用壁面率，而是直接在桥墩壁面上规定不可入条件及无滑动条件。实践证明，此方法可较好地模拟桥墩附近绕流现象。

2.4 模型验证

建立了世业洲河段平面二维水沙模型，模型范围上迄三江口，下至六圩河口，模拟河段全长约40 km。为提高模型预测精度，对模型参数进行了反复率定及验证。图 5为沿程水位和典型断面流速的验证结果。限于篇幅，这里仅给出中长期河床冲淤验证。模型采用2011年10月地形作为初始地形，2015年11月作为验证地形。进口流量和含沙量过程采用大通站实测资料，出口水位由镇江站水位流量关系曲线得到。模型验证了河段经过4个完整水文年后的河床冲淤。由图 6可见，计算与实测冲淤部位和冲淤分布符合较好。

3 桥墩压缩冲刷预测 3.1 水沙过程的选取

预测成果合理性与水沙过程的选取密切相关。针对长江下游河段来水来沙特点，考虑三峡水库蓄水后长江下游水文控制站(大通)水沙变化，从偏安全的角度出发，提出水沙过程的确定方法。

(1) 特大洪水年水沙过程：根据世业洲桥位方案的设计标准(100年一遇洪水设计，300年一遇洪水校核)，选取300年一遇洪水作为特大洪水年。其中，流量过程是以1998年(100年一遇)为基础，根据频率对洪峰流量进行放大；沙量过程则根据三峡水库蓄水后大通站流量和输沙率关系曲线，从工程安全角度出发，采用下包线推求特大洪水年所对应的沙量过程(见图 7)。

(2) 系列水文年水沙过程：采用典型系列年，并考虑特大洪水年的组合。根据大通站1950—2014年水沙资料，三峡水库蓄水以来大通站流量变化不大，且多为中小水年，但输沙量明显小于蓄水前，这一变化反映了三峡水库的蓄水拦沙作用。三峡水库蓄水前，大通站多年输沙量为4.27亿t，蓄水后，2004—2014年均输沙量为1.43亿t，较蓄水前减少了66%(图 8)。2007—2010年能较好地代表三峡建库后2004—2014年水沙特点；据统计2007—2010年径流量为8 372亿m³，与三峡建库后平均径流量8 321亿m³基本相当，年最大流量为64 600 m³/s，最小流量为10 000 m³/s，见图 9。考虑到1998年大洪水的影响，最终确定系列水文年为2007—2010年+1998年(沙量过程用图 7下包线)。

3.2 桥墩压缩后桥址断面水动力变化

桥梁建设引起桥区一定范围内水位与流速分布变化，二者将直接影响到桥址断面单宽流量分布。一般而言，建桥后桥墩间单宽流量有所增加(见图 10)，桥墩上下游受墩柱壅水及掩护影响，单宽流量有所减小。在世业洲右汊左侧，由于桥墩间距较小，阻力较大，使得桥墩间单宽流量减小，同时右侧河道单宽流量增幅较大。

3.3 桥墩压缩冲刷预测

遭遇特大洪水年后，桥位上游由于桥墩阻水作用，流速减小，河床以淤积为主；桥位下游受桥墩挤压水流和掩护作用，河床有冲有淤；桥位断面，因桥墩存在，过水面积缩窄，单宽流量增加，出现压缩冲刷，其中冲刷幅度较大的位置是左汊N3#桥墩(9.46 m)和N2#桥墩(5.92 m)，右汊桥址断面冲刷相对较小，2#辅助墩桥墩2.81 m，3#主墩3.70 m。

文献[17]采用1998年和2013年长江下游大通—吴淞口航行图，对长江下游典型桥梁桥址断面压缩冲刷进行统计(见表 1)，可见，因桥址所处位置和断面形态不一，各桥压缩冲刷深度有所差异，幅度在2.6~8.0 m之间，与世业洲桥位模型预测结果相当。

4 结语

(1) 建立了桥墩压缩冲刷预测模型，对小尺度模型，桥墩边界可采用大系数法与壁面律相结合方法; 对天然大尺度模型，壁面率无法适用，可直接在桥墩壁面上规定不可入条件及无滑动条件。模型参数进行多次率定和验证，具备复验河道年内和年际冲淤变形的能力。

(2) 针对长江下游河段水沙特点，考虑三峡水库蓄水以来大通站水沙条件变化，从工程安全角度出发，提出水沙条件的确定方法，并确定不利水沙过程。

(3) 受桥墩影响，桥位断面出现压缩冲刷，其中冲刷幅度较大的位置是左汊N3#桥墩(9.46 m)和N2#桥墩(5.92 m)，右汊桥址断面冲刷相对较小。通过对长江下游典型桥梁桥址断面压缩冲刷深度统计可见，冲刷幅度在2.6~8.0 m，与模型预测结果基本相当。