桥墩局部冲刷发展过程的三维动网格模拟

文章编号：1672-1683（2017）02-0132-06

冲积河道在遇到阻水构筑物（桥墩、丁坝等）时，构筑物周围河床的局部冲刷对阻水构筑物的稳定有很大的影响。对于桥墩而言，水流在遇到桥墩后，由于桥墩的阻水使得过水面积减小，墩周流速增大，河床剪应力增加，墩周河床沉积物被水流搬运，墩周床面高程逐步降低，并产生冲坑，导致桥墩基础的埋深减小，进而会导致桥梁的倒塌，甚至生命和财产的损失。

桥墩的局部冲刷是一个动态的发展过程，影响因素众多，空间分布具有很强的三维特性，这就使得冲刷模型试验成为以往研究冲刷问题的主要手段。但模型实验存在费用高，无法普遍应用，条件单一，存在模型尺寸效应等不确定因素，数值模拟方法的不断改进使得其作为一种研究手段越来越显示出其不可替代的作用。

近年来，国内外学者针对桥墩冲刷三维性态发展开展了一系列数值模拟研究。Ehteram运用SSIIM软件对桥台的冲刷过程进行了三维模拟，得到了冲刷坑深度和形状并与试验结果进行了比较。Khosronejad对不同横截面形状的桥墩进行了三维动床模拟，采用了流固耦合曲线浸入边界的技术。Kim采用大涡模拟的方法对相邻的两个圆柱形墩的局部冲刷坑进行了模拟，得到的最大冲深位置与试验结果较为一致。韦雁机等基于Open-FOAM开源软件的动网格技术，用输沙率计算床面地形随时间的变化，构建起桩周局部冲刷的动态三维数学模型。祝志文等根据床底泥沙的单宽体积输沙率得到河床高程坐标的瞬时变化，采用边界自适应网格技术修改动边界计算域网格，得到圆柱形桥墩周围局部冲刷坑的演化过程。以上研究对局部冲刷的数值模拟起到了很好的推动作用，在实际应用中多少都存在一些不足的地方，如采用虚拟的浸入边界、地形函数等很难与实际条件一致，大涡模拟或分离涡的模拟计算消耗极大，基于单宽体积输沙率来计算河床的变形，计算过程复杂，涉及到梯度的计算，会使误差增加等。

本文基于CFD计算软件FLUENT的动网格技术和用户白定义函数（UDF）功能，使用基于雷诺平均N-S模型的Realizable k-e湍流模型，将床面瞬时剪应力和临界剪应力带入Van Rijn沉积输运函数，得到床面坐标的变化，通过网格重构和弹性光顺结合的方法来不断修正变形较大的网格，实现了局部冲刷过程的动态模拟。

1数值模型

1.1物理模型的选取

在已有的沖刷试验中，Melville和Raudkivi（MR）对局部冲刷坑发展的三个不同的阶段进行了相对较详细的定量描述。所以本次研究选取MR的经典冲刷试验资料建立数值模型，并进行对比分析。MR试验水槽长19 m，宽0.456 m，水深0.15m，平均来流流速为0.25 m/s。模型布置见图1。MR分别选取初始定床阶段（测定了河床面附近的流速），中间发展阶段（冲刷30 min时，冲坑深度达到0.04 m）和冲刷终止的平衡阶段进行分析，给出了详细的试验结果。试验中前30 min发展较为剧烈，而之后冲刷发展开始缓慢。30 min时的冲刷深度达到总冲刷深度的75%。数值模拟取前30 min进行研究。

1.2计算域及网格划分

试验研究表明，圆柱形桥墩绕流流态以x轴基本呈对称分布，因此取其中一半作为本次模拟的计算域，以缩短计算时间。根据试验布置和模拟要求，将计算域高度设置为15 cm，宽度设置为3d即15.24cm，圆柱上游距桥墩中心为3d，出口处要求尾流充分发展，所以设置下游距桥墩中心为10d^[s]即50.8 cm。取圆柱竖向为z轴方向，床面为x-y面水流方向为x轴正方向。具体见图2。

本次模拟为动床模拟，桥墩周围由于冲刷作用使得局部变形较大，且变形不规则，所以选用适应变形能力较强的四面体非结构化网格。为了提高床面剪应力的模拟精度，在床面设置了0.2 cm（约0.5d50）边界层，且边界层会随着床面结点的移动而跟随移动，这样更大程度上保证了床面剪应力获取的精度。因此本次模拟在床面边界层内为三棱柱体网格，其余部分为四面体网格。由于墩周及靠近床面的部位各物理量梯度较大，采用尺寸函数功能（最小网格尺寸0.3 cm，比率1.2，最大网格尺寸1.5cm）对局部网格进行加密。动床面靠近墩周的部位考虑到湍流边界层及后续发生较大的局部变形，采用尺寸函数（最小网格尺寸0.2 mm，比率1.05，最大网格尺寸1.5 cm）进行加密。整个模型共划分网格单元数127 224个，见图2。

1.3湍流模型

针对湍流求解，大涡模拟（LES）和分离涡模拟（DES）方法对通过桥墩的大尺度涡的动力特性能够精确的预测，但是由于其计算消耗过大，应用到工程中有很大的挑战。冲刷达到平衡的时间尺度（小时或天）比湍流脉动的时间尺度（秒或更小）要大的多，如此大的悬殊使得用LES和DES方法进行冲刷的水动力耦合模拟不太实际。本次模拟采用的湍流模型为更经济实用的雷诺平均N-S模型。以往的研究表明，对于圆柱型墩，雷诺平均N-S模型的缺陷在于其不能有效捕获上游面桥墩与河床相接处的高能湍流涡，而这样的湍流脉动对于冲刷的发展是有影响的。因此可以预见，采用基于雷诺平均N-S模型的局部冲刷模拟在圆柱型桥墩前缘的河床冲刷深度会低于实际值。因此，本次采用动床模拟冲坑深度，可以比较得出基于雷诺平均N-S模型的数值模拟结果与试验结果的误差大小并分析误差的来源。另外，桥墩前缘的湍流脉动，对于桥墩前缘为非圆柱形状的情况，如尖角型，可能会产生不同的流动类型和冲刷动力。所以，桥墩形状，尤其是前缘的形状对数值模拟结果的精度在文末进行了探讨。

1.4边界条件

由于计算域选取的流场入口段距离较小，所以要经过计算给定一个稳定的，边界层充分发展的流速剖面，作为速度入口边界条件。因此，在三维计算之前，首先建立二维无圆柱流场并给定速度入口条件让其充分发展，模拟结果和Melville试验结果均显示出充分发展的速度剖面分布基本符合最广泛使用的karman-Prandtl对数流速分布公式，即：

（1）

底部河床指定为粗糙壁面，根据Melville试验结果，其有效粗糙高度取为2d50。桥墩面，水槽侧壁均设置为光滑壁面。顶面设置为对称边界来模拟自由水面。由于取一半流场进行数值模拟，沿x轴剖分出来的面均设置为对称边界。具体边界设置见图3。

1.5动态网格更新

CFD模拟中，流场的计算采用单相，瞬态求解。在FLUENT中激活动网格，河床面设置为动边界，用DEFINE_GRID_MOTION宏命令来控制边界各个节点的运动。在每个时间步开始计算时，比较床面（动边界）各结点实时剪应力T与床沙起动临界剪应力值τ^σ，若存在超临界剪应力（τ-τ^cr>0），该点表现为冲刷，结点下移，否则表现为静止，结点位置不变。局部冲刷问题，局部网格变形较大且不太规则，随着冲刷坑逐渐发展，局部网格必然变大或扭曲造成数值发散，所以为了保证各区域变形后网格尺寸不至于过大或者过于扭曲而使数值发散，选用局部网格重构与弹性光顺相结合的方式进行动态网格更新，网格变形前后对比图见图4。

在FLUENT软件中，通过用户白定义函数（UDF）获取床面实时剪应力值并存储。平床下床沙起动的临界剪应力值通过由希尔兹公式推导得出的临界剪应力的计算方法：

（2）

（3）

随着冲刷的发展，床面开始出现坡度并逐渐增加，由于床沙重力在水流方向产生分力，床沙的起动临界剪应力将不再等于平床下的临界剪应力值。本文采用Dey提出的经验方程来计算变临界剪应力：

（4）

除了通过剪应力确定各结点在每个时间步是否冲刷产生向下位移外，还需要确定各结点向下位移的大小。在每个时间步内，各结点向下的位移等于时间步长与网格移动速度的积。通过调整时间步长可以控制在每一时间步内网格移动增量处于一个较小的值。底边界网格移动的速度根据Van Rijn基于水槽试验提出的沉积输运函数来表示：

（5）

2模型验证

2.1流场对比

在进行动床模拟前先进行了定床条件下的流场模拟。从图5可以看出，墩前流线的分离和墩后尾涡的形态，墩前垂直剖分面上的下降水流都基本与试验结果吻合。Melville试验定床冲刷时，观测到在桥墩附近流速明显变大，最大流速在迎流面中轴线两侧±100°的位置，数值约为1.5倍的平均来流速度即0.37 m/s。图6是数值模拟得出的流速等值线分布图和床面剪应力分布图，无论是最大流速的位置还是大小都与试验结果基本一致。最大河床剪应力的位置与最大流速的位置基本一致，这也与试验结果基本吻合，说明流速的增加引起了床面剪应力的增加，当床面剪应力超过床沙起动应力时就表现为冲刷。

2.2冲刷过程对比

试验和数值模拟结果均表明桥墩的局部冲刷经历了一个先强后弱的过程（图7）。在冲刷开始阶段，试验和模拟结果高度吻合，冲刷速率很大。当冲刷约7 min后，冲刷速率减缓，试验深度开始大于数值模拟的深度。在30 min时，试验深度达到4 cm，数值模拟深度为3.5 cm，数值模拟最大冲深比试验结果小了13%。

从Melville试验及其他研究者的研究成果可以看出，冲刷首先出现在桥墩两侧（迎流面中轴线两侧45°～100°之间），冲刷的产生是由于桥墩两侧水流的加速引起了床面剪应力超过床沙起动剪应力，从而产生了冲刷。从图8看出数值模拟冲坑也首先出现在桥墩两侧45°～100°之间，与试验结果吻合。冲刷几分钟后，桥墩前端下降水流引起的湍流涡系开始显现并逐渐发挥对床沙冲刷的加速作用，冲坑开始由桥墩两侧往桥墩前端贯通，最终冲坑最深的位置处于桥墩迎流面正前方到两侧75°范围内，如图9试验冲坑等值线所示。数值模拟冲坑等值线整体形态与试验结果相近，主要区别在于桥墩前端冲深比试验结果小。正如1.3节所述，数值模拟采用的雷诺平均N-S模型对圆柱型墩正前端湍流脉动对冲刷的贡献不能有效的模拟出来，使得桥墩正前端的冲坑深度低于試验结果。

从以上分析可以看出，本文所采用的动态网格模拟冲坑的方法能够较好地反应冲刷发展的动态过程，模拟结论可信。

3桥墩横截面形状影响的探讨

从前面的分析可以看出，对于圆柱型的桥墩，产生冲刷的原因主要来源于以下两个方面：（1）由于桥墩的存在，局部流速加大，引起局部河床剪应力超过了临界剪应力，引起墩周局部冲刷；（2）湍流脉动速度引起的瞬时的河床剪应力增大引起的局部冲刷。本文所采用的数值模拟方法是基于雷诺平均N-S模型的，因此对于湍流脉动速度估计不足，从而使得墩头处模拟深度与实际相比偏小。以往的研究表明，墩头钝度不同，湍流脉动强度也不同。比如方形桥墩，钝度最大，圆端型次之，尖角型墩钝度因子（BF）为0。Olcmen和Simpson经过研究发现，湍流脉动和马蹄涡强度随障碍物前缘钝度的减小而逐渐降低，当钝度因子为O时，湍流脉动微弱到无法测得。因此本文所采用的数值模拟方法，对于尖角型墩可以达到准确预测，随着钝度的增大，此方法在墩头处的预测深度会比实际的偏小，已有的试验研究成果也已观测到，对于尖角型墩，其最大冲刷主要位于两侧角部而墩前端的冲刷深度很小。而圆端型和方型的墩最大冲深均处于墩前端及附近一个较小的范围内，这也间接地证明了以上分析的合理性。具体更详细的验证有待后续进一步的研究。

4结论

床面的冲刷类似于柔性变形，因此河床底面坐标的改变要对河床面的各个结点进行控制。各个结点在每个时间步是否下移和下移多少主要取决于两个关键参数：床面剪应力和临界剪应力。通过对模型底部网格设置边界层来提高床面剪应力的提取精度。临界起动剪应力会随着床面坡度的改变而变化，所以采用变临界剪应力的方法更符合实际。各结点运动的速度通过Van Rijn沉积输运函数来获得运动的速度，速度乘以时间步长就是在每个时间步内坐标移动的大小。与求解输沙率来求得床面坐标变化的方法相比，本方法不需求解标量的梯度和微分方程，计算简单，避免了梯度求解产生的误差，冲刷的速率更接近试验结果。本方法在应用时，钝形前端桥墩的最大冲深比实际情况会小约13%，本方法是在FLUENT软件基础上的二次开发，可塑性较高，在后续的研究中可通过尝试增强桥墩前端床面剪应力的方法来减少钝形桥墩前端冲深预测不足的缺点。