基础结构的冲刷病害会直接导致桥梁结构功能失效、丧失其安全性能.美国交通部已将桥梁基础冲刷看作是高速公路桥梁结构功能及安全性能失效的最常见原因之一^[1-3].长期以来，桥梁基础冲刷深度预测主要基于日常检查的主观经验判断，准确性不高.尽管针对个别特大型桥梁进行冲刷模型实验，但模型相似比难以确定，且实验人力物力花费较高.大多规范所采用的冲刷深度计算公式可快速对桥梁基础冲刷病害进行预判，但是计算假设条件苛刻，参数单一，难以准确模拟冲刷三维整体形态，甚至不同计算公式得出的结果差异性很大.

为保证桥墩冲刷分析的精细化、准确性与经济性，计算流体动力学(CFD)数值方法在冲刷分析中得到了越来越多的应用.近年来国内外学者提出了一些数值模型，郭辉等^[4]根据水、沙运动及河床冲淤变形方程，考虑压缩断面上下游出现回流的特点，建立了跨河桥梁基础压缩冲刷的一维数值计算模型.董天乐等^[5]从平面二维非恒定水流运动方程、不平衡泥沙输移方程及河床变形方程出发，建立了钱塘江河口平面二维河床冲刷数值模型.显然，一维、二维的冲刷模型无法对桥墩冲刷深度以及平面形态发展这一三维事件进行正确跟踪，计算结果也无法得到最深冲刷深度对应的平面位置.为解决这一问题，三维数值模型随着计算机能力的迅速提升也逐渐发展起来，但随着模型的不断繁杂，需要不断更新边界条件的动态处理以及流体结构计算算法的收敛性又成为难题，并限制了桥墩冲刷三维精细化分析技术的发展.

事实上，冲刷三维精细化分析的难题大大限制了桥墩抗冲刷优化设计的发展.如今对冲刷的防治主要针对桥墩进行长期监测与发生冲刷之后的桥墩加固设计；或者基于规范计算得到的冲刷深度而调整墩身高度来被动预防冲刷对基础的影响.而很少或没有通过针对不同冲刷环境的桥墩自身型式的优化设计来最大化主动降低冲刷影响；而事实上，这种方法抗冲刷效率优，经济性能好，是未来抗冲刷桥墩设计的重要发展方向之一.显然，这种设计方法需要准确的桥墩冲刷三维精细化模型以及基于该模型的冲刷环境参数精细化分析结果来提供理论支撑.

本文首先基于K-ε湍流模型建立桥墩冲刷及冲刷环境数值模型，对其三维仿真建模过程中的4个关键问题分别进行深入研究并提出优选的解决方法，特别关注三维边界条件的动态更新以及计算收敛性的改进.进而利用B. W. Melville经典实验对本文提出的桥墩冲刷仿真方法从三维形态的角度进行多方面准确性验证.最后基于该仿真方法，对CFD计算软件Fluent二次开发，以进行一系列参数分析.根据计算结果，确定不同冲刷环境参数，包括流场中结构物几何参数(桥墩形式以及桥墩尺寸)以及流场形态参数(河床沙粒半径，水流平均流速以及水流深度)，对桥墩冲刷空间形态的影响.该影响可对桥墩主动抗冲刷选型与设计方法提供理论依据.

1 仿真建模关键问题

仿真建模主要指的是结构本身描述的精细化程度以及与之相对应的计算方法的收敛性.对此，完成桥墩冲刷三维仿真建模至少需要解决以下4个关键问题.

1.1 网格选择

采用有限体积法对水流的三维空间形态进行数值模拟，流场的正确网格划分是关键.虽然非结构化的四面体网格对于复杂模型有着较强的适应性与变形能力；但是基于以下几点原因，本文认为结构化六面体才是精细化冲刷建模网格划分的首选：1)冲刷模型涉及到动网格的动态更新，更新过程非常复杂，对网格质量的要求比较高，因而采用质量较高的结构化六面体网格更容易收敛；2)结构化六面体网格和非结构化四面体网格相比，同样网格尺寸，单元数量较少，因此计算时间相对缩短；3)结构化六面体网格方向更能迎合流场方向，在边界层处六面体网格比四面体网格离散误差更小些. 4)由于桥墩附近水流处于复杂流动状态，水流流动梯度大，需要对网格进行局部加密.为了适应动网格的需要，非结构化四面体网格需要采用ICEM CFD中的“密度核”进行加密，其效果和便利性远远不如结构化六面体网格.因此本文选用结构化六面体网格来完成相应冲刷建模工作.

1.2 网格局部加密方法及加密区域确定

数值计算的结果不仅与网格质量有关，而且与网格划分过程中网格的拓扑结构有关.事实上，网格拓扑往往会导致伪扩散等数值问题；此外在结构网格中网格格线与水流流动方向存在偏角也一定程度上导致伪扩散出现.本文提出3种桥墩周围网格局部加密方式作为候选，如图 1所示.

方法1采用基于桥墩圆柱体位置的O型块直接剖分方式对桥墩圆柱体周围进行局部网格加密，网格畸变情况严重，网格质量最低.方法2将桥墩圆柱部分独立单次竖向切割，然后进行O型块剖分局部加密，有助于提高网格质量；但是圆柱四周网格线显然与水平方向夹角过大.方法3将圆柱周围横竖各切割两次，然后进行O型剖分对桥墩圆柱周边局部加密，从而减少外侧网格倾斜程度及数量，有助于将倾斜网格线的影响降低.其中，切割线距圆柱中心线的距离可取为2.5D，其中D为桥墩直径.

ICEM的determine函数可以用来判断网格质量好坏，即通过数学方法检测所划网格的畸角，最大、最小网格尺寸等指标来对网格质量进行客观评价.检测结果如图 2所示，绿色区域越靠近1表明网格质量越好，可以看出按方法3所划出的网格质量最高，完全符合数值计算的网格要求，并且方法3划分的网格质量远优于其他两种划分方法.据此，本文推荐方法3对桥墩处网格进行局部加密，加密范围可取以桥墩为中心的2.5D边长的正方形区域，经过计算Y+(垂直壁面的无量纲距离)均满足所使用的湍流模型要求.

1.3 床面变形求解及网格更新策略

根据输沙平衡，由推移质输沙率得到河床高程的瞬时变化率求解公式为^[6]

(1)

式中：h为河床高程；t为时间；n为河床泥沙孔隙率，本研究参考实验资料取0.41^[7]；q_b为水平床面推移质单宽体积输沙率；q_bx为x方向的推移质单宽体积输沙率；q_by为y方向的推移质单宽体积输沙率.

对于某一特定面单元，式(1)中面单元中心高程的变化梯度∂h/∂t可以利用下式进行离散化处理，即

(2)

式中：Δh, Δq_bx, Δq_by为离散化后的各参数在Δt内的增量；Δx, Δy为两个距离最近单元之间的坐标差值.

通过对网格面区域循环遍历所有的面单元，找出与该特定面单元距离最小的面单元，以两者之间的差值代入式(2)近似代表式(1)中各梯度变化项. Δt的取值根据计算精度及效率确定，本文取为0.075 s.

当采用CFD计算软件Fluent进行动态网格更新时，其需要通过改变各网格节点的坐标来实现.因此上述河床面网格中心点的高程变化值需要转化为相应节点坐标的变化值.本文通过计算节点周围网格面单元中心点高程变化的平均值来得到该相应节点坐标变化值.

以图 3中的节点A为例，其坐标变化值Δz_A为

(3)

式中：δz_i为网格单元i中心点的高程变化值；A_i为网格单元i在xy面内的投影面积.

1.4 河床面坍塌程序化处理

许多水槽试验表明，由于泥沙运动的复杂性和与水流之间的非线性作用、运动关系，会出现局部区域坡度超过泥沙休止角的情况，此时河床面会自动崩塌来调整坡度，从而逐渐形成冲刷坑.数值模拟时，若不考虑床面崩塌情况，计算结果显然与实际不符，造成数值模拟的失真；同时该区域网格将严重畸变，计算过程很容易发散.

解决这一问题，本文采用沙滑模型对河床面坍塌进行程序化处理.首先根据输沙平衡及床面变形方程对床面网格进行动态更新，进而对床面网格进行全面扫描，若发现某单元面倾角大于临界坡度(本文取泥沙休止角)，则对该单元节点进行循环操作来完成位置修正.以图 4为例，一旦发现A、B节点的倾角θ大于临界坡度，则将A点纵坐标调整到A′，将B点纵坐标调整到B′，从而使坡度降到，以满足床面变形稳定的要求，该变化可表示为

(4)

式中：z_A、z_B分别为A、B两点的纵坐标；δz_A、δz_B分别为A、B两点需要调整的纵坐标值.

同时，根据输沙平衡又可得到

(5)

式中：∑A_A为与节点A相邻的各面单元在xy平面内的投影面积之和；∑A_B为与节点B相邻的各面单元在xy平面内的投影面积之和.

进而联立式(4)、(5)，即可求得A、B两点的调整高度δz_A、δz_B.当采用CFD计算软件Fluent进行程序编译时，通过循环宏begin_f_loop(f, tf)对底面所有单元网格面进行循环，并采用循环宏f_node_loop(f, tf, n)对所循环到的面内的所有网格节点进行循环，依次判断面内相邻节点间坡度是否满足要求，并对大于休止角(临界坡度)的节点纵坐标进行调整，即可实现泥沙坍塌的过程模拟.但是，由于调整一个河床单元面坡度显然会改变其相邻的其他网格面坡度，所以必须通过多次迭代才能使所有河床单元面的倾角均降至休止角或以下.在实际操作中，逐步调整每个面单元倾角到休止角以下耗时太多；本文对整个河床面循环迭代操作次数限制在100次以内，即可同时满足计算精度与计算时间的要求.

2 数值模型及准确性验证

为验证上节所提出关键问题及相应解决方法的准确性，基于Melville经典冲刷试验建立相应的数值计算模型.由于Melville冲刷试验参数明确，试验数据完备，同时利用该试验数据对本文所提出的冲刷数值模型计算精度进行验证. Melville试验水槽长19 m，宽45.6 cm，水深0.15 m，在水槽中放置直径为5.08 cm的圆柱作为桥墩模型，圆柱型桥墩中心距水槽两侧的距离为22.8 cm，床底泥沙平均粒径d₅₀为0.385 cm，水流平均速度为0.25 m/s，泥沙休止角为32°^[8].事实上，根据Sarker的研究成果，圆柱下游12倍直径距离以外流动并不受圆柱影响^[9].故本文建立的数值计算模型所选取的试验范围总长度为20D=101.6 cm，其中D为圆柱桥墩直径，桥墩中心距下游出流面14D，满足大于12D的要求，桥墩中心距上游进口面6D，如图 5所示.河床面采用粗糙壁面边界，两侧以及桥墩采用光滑壁面边界，顶部采用对称边界.相应数值计算模型见图 6.基于1.2节所提出的加密方法，该模型中桥墩周围一定区域范围内已进行网格加密.

2.1 冲刷坑形状验证

图 7(a)为模型计算30 min后桥墩周围局部冲刷坑的地形等高线图，图 7(b)为Melville试验30 min后局部冲刷坑的地形等高线图.比较两者可以发现，数值模拟结果除了桥墩迎水面冲刷坑形态发展不够充分外，其他包括桥墩后方沙脊的出现，总体轮廓线的形式均能够与实验得出的冲刷坑地形图较好吻合.

事实上，由于本次计算选用了标准k-ε模型，对于流场的变化主要基于一定时间空间内的平均化结果，因此可能会引起计算结果在某些局部区域的偏差，如在主要冲刷坑外同时也出现了独立的小幅度冲刷坑，如图 8所示.

2.2 桥墩周围流场形态验证

随冲刷不断发展，流场底部的边界部分(河床面)会相应不断下移，从而引起边界条件的不断变化.该边界条件的变化反过来又会引起桥墩周围流场的不断变化，两者之间相互耦合.显然，局部冲刷是一个流体与边界条件不断相互作用的复杂变化过程.所以，能否准确模拟桥墩周围流场变化对保证桥墩局部冲刷计算精度有着至关重要的关系.

图 9(a)为冲刷仿真开始时刻距离河床面1 mm的截面的流场流线平面图，图 9(b)为Melville试验中的实际量测结果.通过比较可以发现，两者的流线轨迹符合较好，且仿真计算的尾流与桥墩的分离点位置(偏离水流方向75°左右)与试验结果完全一致.

2.3 冲刷坑深度验证

Melville试验结果表明冲刷达到最终平衡需要时间为220 min，此时平衡冲刷最大深度为6 cm；而事实上，试验进行到30 min时冲刷深度为4 cm，已达到最大深度的2/3.本次研究基于计算机硬件条件限制以及计算时间的考虑，仅模拟前30 min的冲刷状态，并将仿真计算的冲刷深度与试验结果进行比较.

图 10给出本次仿真计算中冲刷坑最大深度随时间发展的变化趋势.从图中可以看出在冲刷发展初期，冲刷速度很快，最大冲刷深度随时间的推移急剧增大.其原因在于，冲刷初期河床面最大切应力比临界起动切应力大得多，导致冲刷发展剧烈.而后期由于底部边界下切，流场运动空间增大，导致流动逐渐减弱，水流携沙能力降低，最终冲刷变缓.根据图中曲线斜率可以看出，如果继续计算，该曲线将趋近于某一稳定数值，即平衡冲刷深度，这与Melville试验所描述的现象相符.另外，仿真计算30 min时得到的最大冲刷深度为4.704 cm，仅仅比试验结果4 cm稍大而已，显然具有较好的计算精度.

2.4 冲刷坑坡度验证

30 min仿真计算后的桥墩前方冲刷坑坡度(图 11中直线AB倾角)为31.9°，此值与Melville试验的坡度值32°(同时也是泥沙休止角)非常接近.事实上，由于本次仿真模型计算并未达到最终冲刷平衡状态，比泥沙休止角偏小也属正常.另外，仿真计算中的桥墩后方冲刷坑坡度仅为24.9°(图 11中直线CD倾角)，这是由于泥沙在桥墩后方存在少量淤积，导致冲刷程度较弱，最终形成的冲刷坑坡度也较小.

通过与经典Melville试验的多方面验证可以看出，本文所提出的桥墩冲刷精细化建模方法具有较高的计算精度，计算结果与试验数据在各方面均具有一致性，完全可以利用该数值模型对局部冲刷行为进行准确预测.

3 环境参数对桥墩冲刷空间形态影响分析 3.1 环境参数选取

进一步利用所提出的仿真模型计算方法探讨环境参数对桥墩冲刷空间形态的影响.选取的主要环境参数可分为流场中结构物几何参数以及流场形态参数.其中，流场结构物几何参数包括桥墩形式以及桥墩尺寸；流场形态参数包括河床沙粒半径，水流平均流速以及水流深度.

3.2 流场中结构物几何参数分析 3.2.1 桥墩形式

桥墩形式对局部冲刷的影响，考虑到实际情况，主要采用应用较多的单柱墩、双柱墩以及排墩这3种形式.此处仍然基于Melville试验的尺寸与条件，计算数值模型如图 12所示(其中单柱墩计算数值模型见图 5).

经过计算得到30 min后单柱墩的最大冲刷深度为4.704 cm，双柱墩的最大冲刷深度为5.297 cm，排墩的最大冲刷深度为5.101 cm，三者冲刷深度随时间的变化趋势如图 13所示.

图 13可以看出双柱墩和排墩最大冲刷深度随时间发展的变化趋势与单柱墩类似，初期冲刷发展剧烈，后期逐渐变缓，最终达到冲刷平衡状态.但是有所不同的是，冲刷初期(图 13中5 min之前)单柱墩的冲刷深度比其他两者均稍大，但是到了中后期(图 13中5 min之后)却明显小于其他两者.这主要是因为，冲刷初期发展并不充分，单柱墩在平行水流方向的长度较小，导致墩周围平均流速较大，冲刷深度较深；但到中后期冲刷发展已较为充分，墩的阻水效应逐步显现，显然较大尺寸的双柱墩与排墩对水流的扰动更强.因此，此时平衡冲刷深度反而比单柱墩深.

根据单柱墩、排墩以及双柱墩30 min时冲刷坑的空间形态计算结果，可以看出单柱墩冲刷坑最大位置位于桥墩中心线偏后的位置，排墩冲刷坑最大深度位置位于排墩中心线两侧位置，而双柱墩最大冲刷深度位置(两个峰值)分别位于两墩中心线偏后的位置.双柱墩模型中其后方桥墩位置处的冲刷深度比前方小，主要是因为前方桥墩对水流的阻挡作用降低了水流强度而导致的.

所以，通过选择合适的桥墩形式可以达到降低冲刷深度，调整最大冲刷深度所对应平面位置的优化目的.

3.2.2 桥墩尺寸

桥墩尺寸大小显然会影响其对水流的阻挡作用，进而影响冲刷程度.本文仅针对单柱墩模型，分别计算桥墩直径为0.5D、0.8D、1.0D、1.2D、1.5D以及2.0D共6个工况的冲刷过程.

通过计算可以发现，冲刷深度h随桥墩尺寸的增大而逐渐变大，每一种工况下冲刷随时间的变化趋势也均相同，进一步验证了本文提出的仿真计算方法的稳定性. 30 min后各工况下最大冲刷深度分别为2.968 cm(0.5D)、4.382 cm(0.8D)、4.704 cm(1.0D)、5.543 cm(1.2D)、6.287 cm(1.5D)以及6.916 cm(2.0D).

为更好描述冲刷深度与桥墩尺寸的关系，利用MATLAB软件对计算数据进行公式拟合，其中横坐标为归一化后的桥墩尺寸α(D的倍数)，纵坐标为30 min时的最大冲刷深度h(cm)，得到两者近似关系为

(6)

可见，随桥墩尺寸的不断增大，桥墩局部冲刷深度将以近似线性的关系逐步增大.因此，实际工程中降低桥墩有效宽度是一种较为明显的减小桥墩局部冲刷深度的方法.

3.3 流场形态参数分析 3.3.1 水流速度

水流流速是影响桥墩局部冲刷的一个重要因素，因为水流流速直接影响到水流动能大小.一般认为，水流流速越大，桥墩周围局部冲刷程度越剧烈，这也是洪水季节容易发生桥墩周围泥沙掏空，桥墩倾覆的重要原因.

本文在保证其他各因素不变的情况下，仅仅改变水流流速，分别计算了0.200、0.225、0.250、0.275、0.300、0.325 m/s共6个工况的冲刷过程，得到30 min时各工况下最大局部冲刷深度分别为2.383、3.304、4.704、5.123、6.399、6.672 cm.值得注意的是，在计算过程中发现当水流流速较小时(小于0.15 m/s)，桥墩周围将几乎不产生局部冲刷现象，这是由于当流速小于起冲流速时，水流将不具备足够的动能使泥沙发生运动.

利用MATLAB软件进行公式拟合可以定量分析流速v(m/s)对最大冲刷深度h(cm)的影响，得

(7)

可见流速与最大冲刷深度之间存在近似的线性关系.洪水季节，水流流速急剧增大，此时可采取措施降低墩前水流流速来显著降低桥墩局部冲刷程度.

3.3.2 水流深度

本文针对0.02、0.04、0.06、0.08、0.10、0.15、0.18、0.20、0.22、0.25 m共10组水深分别进行计算，其最大冲刷深度分别为3.518、4.506、4.826、4.971、5.123、4.704、4.790、4.719、4.783、4.686 cm.

由以上计算结果可明显看出，当d₀/D>2.6(d₀为墩前水深，D为桥墩有效宽度或直径)的一组水深时(水深大于0.15 m)，冲刷深度均在4.7 cm左右，近似可认为冲刷深度与墩前水深关系不大.而对于d₀/D < 2.6的一组水深(水深小于0.1 m)，由MATLAB软件进行公式拟合得到最大冲刷深度h(cm)与墩前水深d₀(m)之间的关系为

(8)

由此可见，当d₀/D < 2.6时，局部冲刷深度随墩前水深的增大而逐渐增大，且两者之间满足指数增长关系.

根据Melville和Sutherland研究，存在一个水深数值，当大于该数值时，局部冲刷深度基本与水深无关^[9].这一结论也与本文的仿真计算结果相吻合.

所以，在桥墩设计过程中，可以通过选择合适的桥墩位置，采取合理的d₀/D数值，来降低桥墩局部冲刷程度.

3.3.3 河床沙粒半径

事实上，已有国内外学者关于泥沙粒径对局部冲刷深度的影响进行了一系列试验研究. Laursen曾经用中值粒径d₅₀=0.46~2.2 mm的床沙进行动床实验，结果表明泥沙粒径对冲刷深度无明显影响^[10]. Ettem也指出，当泥沙中值粒径满足d₅₀≤D/50(其中D为桥墩直径或有效宽度)时，泥沙粒径大小对局部冲刷深度无影响^[11].鉴于本文提出模型中，若考虑泥沙粒径对局部冲刷深度的影响，需要控制泥沙的均匀性不变^[12].此外，改变泥沙粒径后，相应的泥沙水下休止角及容重均将改变，而这些数据在数值计算模型中的相应改变缺少必要的试验依据和理论指导.因此，依据相关试验，可以认为粒径对局部冲刷深度无影响.

4 结论

1) 通过选择合适的桥墩形式，可以达到降低冲刷深度，调整最大冲刷深度所对应平面位置的优化目的.

2) 桥墩尺寸与最大冲刷深度存在近似线性关系.实际工程中降低桥墩有效宽度是可以明显减小桥墩局部冲刷深度.

3) 流速与最大冲刷深度存在近似线性关系.洪水季节，水流流速急剧增大，需采取措施降低墩前水流流速来降低桥墩局部冲刷程度.

4) 桥墩设计过程中，可以通过选择合适的桥墩位置，采取合理的d₀/D数值，达到降低桥墩局部冲刷程度的设计目的.

5) 后续将集中研究流固两相流模型在桥墩动床冲刷环境中的应用方法，并建立仿真模型与本文清水冲刷计算结果进行比较.