水流作用下圆柱局部冲刷三维数值模拟

李绍武，杨 航

(天津大学 水利工程仿真与安全国家重点实验室，天津 300072)

墩柱绕流现象所带来的墩柱周边冲刷问题是实际工程所不容忽视的一大问题，国内外有许多研究成果，采用数值模拟方法的研究工作起源于墩柱周围的流场模拟。1992年，Kobayashi[1]通过将二维离散涡模型进行扩展，提出了三维涡段模型，尝试了在振荡流的情形下模拟圆柱周围的三维流场分布，模拟出了尾涡的不稳定特征，并和实测的流速场分布吻合较好。1993年，Olsen 和 Melaaen[2]首次将数值模拟运用于计算三维圆柱周围局部冲刷问题，他们通过引入泥沙输运方程，模拟了清水冲刷条件下圆柱周围无粘性沙局部冲刷过程的初期发展阶段，计算结果与实验冲刷型态能够较好地吻合。1998年，Olsen和Kjellesvig[3]在1993年模型的基础上进一步开展研究，模拟的冲刷过程不再仅局限于冲刷的初始阶段，而是扩展到了冲刷的全过程，最终计算所得的清水条件下圆柱周围的冲刷深度与经验公式计算结果吻合良好。2003年，Catalano和Wang[4]在高雷诺数的条件下运用大涡模拟研究了墩柱周围的绕流问题，通过与RANS紊流模型所得计算结果进行比对后发现，大涡模拟所得的计算结果能够更好地反应高雷诺数情况下拖曳系数减小和边界层分离推后的现象，但是由于大涡模拟对计算资源的要求较高，并未得到广泛运用。2005年，Roulund和Sumer[5]基于SST k-ω紊流模型模拟了圆柱周围的三维冲刷问题，较好地模拟出了床面沙纹的产生以及冲刷坑的形成过程，模拟结果中的平衡冲刷深度与实验值吻合良好。

由上述国内外研究进展可以看出，前人在水流作用下圆柱周围绕流流场和局部冲刷上已取得了显著进展，但仍存在着一些问题。一方面，相对于其他紊流模型而言， LES可以更好地模拟圆柱绕流流场，但是，由于LES模拟的流场结果依赖于网格尺度，其对于局部冲刷结果产生的影响尚不明确。另一方面，前人的研究成果大多局限于小尺度的圆柱，不同尺度圆柱周围的冲刷特性存在何种差异也需要进一步探讨。本文运用FLOW-3D软件基于大涡模拟模型，对不同直径的圆柱周边的流态及冲刷过程进行数值模拟探讨，研究圆柱尺度变化及网格大小对冲刷坑深度及形态的影响。

1 计算模型及其数值理论

1.1 基本控制方程及模型理论

LES过滤后的连续性方程

(1)

LES过滤后的N-S方程

(2)

(3)

涡旋粘性模型

(4)

(5)

Smagorinsky模型

(6)

LS=CS(δxδyδz)1/3

(7)

(8)

在冲刷模型采用基于Mastbergen和Van den Berg[6]的经验输沙模型。首先计算无量纲数d*,i

(9)

式中：di为第i种泥沙的颗粒直径；ρf为流体的密度；ρi为第i种泥沙的密度；||g||为重力加速度的范数；μf为流体的动力粘滞系数。

根据Soulsby-Whitehouse公式[7]计算θcr,i

(10)

根据泥沙休止角对临界谢尔兹数进行坡度修正[7]

(11)

式中：ψ为水流方向与逆坡方向的夹角；β为床面的坡度；φi为第i种泥沙的休止角。

本地谢尔兹数通过底部切应力τ求解

(12)

泥沙颗粒的挟带抬升速度ulift,i由下式计算[6]

(13)

式中：αi为挟带系数，本文中取0.018[6]；ns为床面的外法线方向。

本文使用的泥沙沉降速度usettling,i计算公式来自Soulsby[7]

(14)

式中：vf是流体的运动粘滞系数。

推移质输沙率采用Meyer-Peter的计算模型[8]，公式为

(15)

(16)

式中：Φi为无量纲的推移质输送率；βi为推移质输沙系数，取为8[9]；cb,i是第i种泥沙在床面泥沙中所占的体积分数。

悬移质泥沙浓度采用输运方程来求解，方程为

(17)

式中：Cs,i是第i种泥沙的悬移质质量浓度，即单位体积的水沙混合物中泥沙的质量；us,i为悬移质速度；D为扩散系数。

悬移质速度us,i即可按下式计算

(18)

1.2 计算域网格划分

由于模型计算过程颇为耗时，因此在计算域的选择上，在不影响正常计算结果的前提下，尽量选取较小的计算区域。在正式计算前通过一段40 m长的数值水槽获得充分发展的紊动水流，并通过热起动将充分发展的紊流条件作为计算域的入流边界条件。冲刷模型中，圆柱中心位于距入流边界8D处，距离出流边界10D，左右边界和圆柱中心的距离均为5D。计算区域设置及网格划分如图1所示。

1-a 平面图 1-b X-Z剖面图图1 计算域设置及网格划分Fig.1 Setup of computational domain and mesh discretization

1.3 模型参数及边界条件设置

分别选用直径为0.1 m，0.2 m，0.3 m及0.5 m的圆柱进行模拟，模型底部铺设均匀床砂，粒径d为0.26 mm，密度ρs为2 650 kg/m3，休止角α为32°，床面粗糙度取2.5倍中值粒径。入口处和出口处均布置一个与床砂等高的挡水护板，为保证其与床砂的交界处平稳过渡，设置其表面粗糙度为0.65 mm。床砂区上层为水，水体温度取20°C，水深为0.4 m，断面平均流速为0.46 m/s。重力加速度取9.81 m/s2。

为在保证精度的情况下减少网格总数以减少计算量，模型采用嵌套网格，共分为三层，相邻层网格的尺寸比例为1∶2。计算域的底部边界采用壁面边界条件，法向流速为零；侧边界和顶部边界采用对称边界条件，此边界处流量为零，流体剪切应力为零；入流边界采用网格重叠边界条件，将充分发展的断面紊动流速作为入流条件；出流边界采用压力边界，控制出流边界的水深。

2 模型验证

2.1 水流模型验证

参照Roulund[5]的物理模型实验，建立光滑床面下的水流数值模型。模型中，圆柱直径D为0.536 m，水深h为0.54 m，平均入流速度V取0.326 m/s，并经一段40 m长的无墩柱水槽的模拟，获得满足对数流速率的入流边界水流流速分布。

图2和图3分别展示了墩柱周围水流的x向流速及z向流速模拟结果与实测数据的对比，选取的水深位置分别为距床面0.005 m，0.01 m，0.02 m，0.05 m，0.1 m及0.2 m处，可以看出在不同的水深上，模拟的x方向流速和z方向流速与实测数据均吻合较好。

2-a z=0.005 m 2-b z=0.01 m

2-c z=0.02 m 2-d z=0.05 m

2-e z=0.1 m 2-f z=0.2 m图2 墩柱周围x向流速模拟值与实验值对比Fig.2 Comparison of modeled velocity in x direction around the pile with measured results

3-a z=0.005 m 3-b z=0.01 m

3-c z=0.02 m 3-d z=0.05 m

3-e z=0.1 m 3-f z=0.2 m图3 墩柱周围z向流速模拟值与实验值对比图Fig.3 Comparison of modeled velocity in z direction around the pile with measured results

2.2 泥沙冲刷模型验证

参照Roulund[5]的物理模型实验，建立泥沙冲刷数值模型。模型中，圆柱直径D为0.1 m，水深h为0. 4 m，平均入流速度V取0.46 m/s，模拟时长为60 min。

图4-a是不同时刻圆柱迎水侧冲刷深度模拟值与实测值的对比。由图可以看出，在冲刷的初始阶段，冲刷坑深度的模拟值与实验值相比略低，但随着时间的推移，在冲刷过程的中段直至接近平衡，模拟值和实验值基本保持一致。

图4-b是不同时刻圆柱背水侧冲刷深度模拟值与实测值的对比。根据图中结果，背水侧和迎水侧类似，在冲刷初始阶段模拟值低于实验值，但在冲刷过程的中后期，二者大体吻合。

4-a 迎水侧 4-b 背水侧图4 墩柱迎水侧和背水侧冲刷深度模拟值与实验值对比图Fig.4 Comparison of modeled scour depth on the upstream side and rear side of the pile with measured results

图5 冲刷坑纵向剖面模拟值与实验值对比图Fig.5 Comparison between modeled results of longitudinal scour profile with experimental results

图5是沿圆柱对称面的冲刷坑纵剖面图。从图中冲刷坑剖面的轮廓线中可以看出：在圆柱的迎水侧，冲刷坑的最大深度以及其总体轮廓都能与实验数据良好契合；而在圆柱的背水侧，冲刷坑的最大深度和实验数据吻合良好，但离开圆柱后，模拟所得的冲刷坑较实验测得的冲刷坑略深一些。

综上所述，在时间上，模拟值在冲刷过程的初始阶段与实验值相比较低，而在冲刷的中后期可以与实验值良好契合；在空间上，模拟的冲刷坑在圆柱迎水侧与冲刷坑实测结果基本吻合，而在背水侧圆柱根部吻合较好，离开圆柱模拟值偏大。总体而言，所建立的泥沙冲刷模型能够较好地与实际情况相吻合。

3 冲刷模拟结果及分析

保持模型的水深h为0.4 m，平均入流速度V为0.46 m/s，均质床砂粒径d为0.26 mm，砂密度ρs为2 650 kg/m3，休止角α为32°，床面粗糙度为0.65 mm不变，圆柱表面设为光滑，直径分别取为0.1 m，0.2 m,0.3 m及0.5 m，计算域的大小及沙床厚度进行相应调整。图6是各组模型在t=1h时冲刷坑沿水流方向对称面的剖面图。

6-a D=0.1 m 6-b D=0.2 m

6-c D=0.3 m 6-d D=0.5 m图6 t=1 h时冲刷坑纵向剖面图Fig.6 Longitudinal scour profile after 1 hour scour

由于冲刷模型的计算量巨大，尤其是对于大墩柱冲刷模型，因此D=0.1 m,0.2 m及0.3 m时，模型的模拟时长为1.5 h，D=0.5 m时，模拟时长取为1 h。之后，统一运用指数函数曲线进行非线性拟合并推算至t=5 h处，再进行比对分析。

在冲刷的初始阶段，无论是迎水侧或是背水侧，冲刷坑的发展过程都比较迅速，而随着时间的推移，迎水侧和背水侧的冲刷速率都逐渐趋于平缓，逐渐达到平衡。在马蹄涡的影响下，圆柱的迎水侧产生持续且大范围的泥沙起动，进而产生冲刷坑并不断扩大。从圆柱迎水侧带走的泥沙起初堆积在圆柱的后侧，但随着冲刷过程的持续进行，圆柱迎水侧的冲刷坑深度不断增加，与此同时，冲刷坑形态围绕圆柱向着背水侧不断发展，并使得原本堆积于圆柱后侧的泥沙向着下游继续推移。

对于直径为0.5 m的圆柱，当x和y方向网格尺寸分别为0.04、0.03和0.02 m时，计算所得的迎水侧和背水侧的冲刷深度出现了一定程度上的差异，如图7所示。

7-a 迎水侧 7-b 背水侧图7 网格尺寸造成的冲刷深度差异Fig.7 Discrepancies of scour depth caused by grid size

图7结果显示网格尺寸对计算结果有一定的影响，但随着网格的细化，计算结果的差异呈收敛趋势。考虑到再细化网格模型计算量过大(采用0.02 m的网格尺寸时，全计算域网格数达300余万，32核工作站模拟1 h花费35 d)，以0.02 m网格尺寸所得计算结果近似代表实际值。

图8和图9分别是各组直径的墩柱迎水侧和背水侧处的冲刷深度随时间的变化趋势。

图8 各组直径的圆柱迎水侧冲刷深度图Fig.8 Scour depths at the upstream side of the pile for different pile diameter图9 各组直径的圆柱背水侧冲刷深度图Fig.9 Scour depths at the rear side of the pile for different pile diameter

从图中可以看出，在冲刷速率随时间的变化特性上，迎水侧和背水侧大体一致：在冲刷的初始阶段，圆柱的两侧均有较大的冲刷速率，随着冲刷过程的进行，两侧的冲刷速率均逐渐减缓并最终趋于零，此时圆柱周围的冲刷坑趋于平衡状态，从图中可以直观看出，各组模型在经过不同冲刷时间后，均呈现冲刷速率趋缓，达到平衡的趋势。

对于同一直径的圆柱而言，在冲刷深度上，整个冲刷过程中圆柱迎水侧的冲刷深度始终大于背水侧；在平衡时间上，圆柱上下游两侧冲刷到达平衡状态所需的时间大体一致。对于不同直径的圆柱而言，随着圆柱直径的增大，圆柱周围的冲刷坑达到平衡状态的时间不断延长，相应地，圆柱迎水侧和背水侧的冲刷深度也随着圆柱直径的增加而增加，然而二者之间并非简单的线性关系。由图9可以看出，虽然在最终平衡状态下，直径0.5 m圆柱的背水侧冲刷深度大于直径0.3 m圆柱，但在冲刷过程的前期阶段，前前者的冲刷深度却始终小后者。造成这种现象的原因在于，在圆柱背水侧发生冲刷作用的同时，从圆柱迎水侧冲刷下来的泥沙在水流的挟带作用下，在背水侧产生了堆积，导致较大直径的圆柱在冲刷过程的初始阶段，其背水侧的冲刷深度可能小于较小直径圆柱。

图10 各组直径的圆柱对称侧冲刷深度图Fig.10 Scour depths at lateral side of the pile for different pile diameter

图10是各组直径的墩柱横向两侧的冲刷深度随时间的变化趋势。可以看出，不同直径的圆柱，其横向两侧的冲刷过程与上下游两侧颇为相似，达到平衡状态时的冲刷深度也大体相同，但冲刷过程与背水侧有所不同。

表1给出了按照各组模型的条件参数，根据Breusers[10]公式、中国公路工程水文勘察设计规范[11]、Sheppard[12]公式计算所得的最大冲刷深度，并与拟合拟合平衡深度数据进行比对。

从表中可以看出，拟合平衡冲刷深度结果与Breusers给定的预测最大冲刷深度较为接近。

表1 各计算公式所估算的最大冲刷深度Tab. 1 Maximum scour depth estimated using different formula

注：h=0.4 m,U=0.46 m/s

表2记录了各组模型经拟合到达平衡状态时，迎水侧，背水侧及对称侧的冲刷深度。

从表2中可以看出，对于不同直径的圆柱，其到达平衡状态时，背水侧冲刷深度与迎水侧冲刷深度之比以及横向两侧冲刷深度与迎水侧冲刷深度之比均保持为一常数，前者约为0.7，而后者在0.95～1.0之间。可以得出对于不同的圆柱直径，圆柱背水侧及横向两侧的平衡冲刷深度与迎水侧的平衡冲刷深度始终保持着特定的比例关系。

表2 各组模型到达平衡状态时的冲刷深度Tab.2 Modeled results of scour depth at equilibrium status

图11 各组直径圆柱迎水侧冲刷深度与直径比值S1/D对比图Fig.11 Variation of S1/D with diameter D

另，从表2和图11可以发现，虽然圆柱迎水侧的冲刷深度随着圆柱直径的增大而增大，然其增长率随着圆柱直径的增大而逐渐减小，从各组模型迎水侧的平衡冲刷深度与圆柱直径的比值可以看出，相对于直径为0.1 m的圆柱，直径为0.2 m的圆柱的S1/D值下降了24.6%，而直径为0.3 m和0.5 m的圆柱，其S1/D则分别下降了32.6%和44.7%。

4 结论

本文利用FLOW-3D三维模拟软件中的紊流模型和泥沙冲刷模型对不同尺度的圆柱局部冲刷过程进行了数值模拟。计算结果表明，对于不同直径的圆柱，其迎水侧、背水侧以及横向两侧的平衡冲刷深度始终保持着固定的比例关系。通过建立不同尺度的圆柱冲刷模型，开展了圆柱局部冲刷数值进行模拟，其中最大的圆柱直径达到了0.5 m，计算结果表明，圆柱周围的局部冲刷深度随圆柱直径的增大而增大，但其与圆柱直径的比值则随着圆柱直径的增大而减小。