基于大涡模拟的圆端形沉井下沉过程水力特性研究

肖苡辀 ，吴启和 ，高宁波

（1.中交第二航务工程局有限公司，湖北 武汉 430040；2.长大桥梁建设施工技术交通行业重点实验室，湖北 武汉 430040；3.交通运输行业交通基础设施智能制造技术研发中心，湖北 武汉 430040）

0 引言

沉井作为可靠的深水基础形式在水工、港航、桥梁建设中发挥着重要的作用，确保沉井基础安全平稳下沉是沉井施工的关键问题[1-2]。在沉井下沉前，需要结合施工水域的水文资料对沉井周边水动力特性进行详细的分析计算，计算结果将可视化呈现三维流场结构，直接为沉井下沉时机、助沉措施选择提供参考数据和意见。

常用的沉井基础形式包括圆形、矩形、圆端形以及其他形式。目前关于圆形沉井和矩形沉井的水动力研究较多，而对圆端形沉井的水动力特性论述较少。水动力特性的分析方法包括原型观测、模型试验、计算流体力学（Computational Fluid Dynamics，简称CFD）以及解析解等。前人对沉井水动力特性试验研究一般包括沉井下沉过程中底部垂线流速分布、侧壁摩阻力分布规律的分析[3-4]。近些年来，CFD方法广泛应用在基础性研究和工程实际中，能够较为准确高效地得到物理量数据，也能够可视化地展示研究成果[5-6]。张永涛等[7]通过CFD软件对矩形沉井基础浮运阻力进行了计算，探讨了波流、自由面以及流向角等因素对浮运阻力的影响。沉井在下沉过程中对水流产生阻滞作用，形成绕流，因此对沉井下沉的水动力研究也可参考绕流水体水力特性研究成果，从而更加系统地阐述沉井下沉过程中的水动力特征。槐文信等[8]采用大涡模拟方法（Large eddy simulation，简称LES）对圆柱绕流的尾流特性进行了研究，表明绕流的尾流特性会对物体自身结构稳定有着重要的影响。Prasad和Williamson[9]对阻流体绕流后的剪切层不稳定性进行了研究，揭示了剪切层不稳定性与雷诺数之间的关系。综上所述，现有沉井的水动力研究侧重于阻力、流速等方面，而对相干结构、尾流特性的关注较少。

常泰长江大桥5号墩沉井为目前世界上规模最大的水中沉井基础，下沉难度大、精度控制高、水流结构复杂，保证超大型沉井安全平稳下沉是施工的关键问题。本文采用大涡模拟方法，依托常泰长江大桥5号墩沉井基础设计，对圆端形沉井下沉过程中流场进行了模拟研究，将数值计算结果与试验数据进行了对比，验证该数值方法的可靠性；并且分析其时均流场特性以及相干结构，为沉井下沉施工方法的选择提供参考数据。

图1 三维模型及沉井尺寸Fig.1 Three-dimensional model and size of the caisson

1 数值计算方法

1.1 控制方程

在大涡模拟方法中，对于不可压缩流体，经过滤波处理的连续和动量方程如下所示[10]。

连续方程和动量方程：

式中：xi、x（ji=1，2，3；j=1，2，3；xi≠x）j均为空间坐标系中x、y和z三个坐标轴方向上的长度分量，m；=1，2，3；j=1，2，3；均为滤波后的速度分量，m/s；滤波后的压强，Pa；μ 为动力黏滞系数，kg（/m·s）；ρ为流体密度，kg/m3；t为时间，s；τij为亚格子尺度应力（Subgrid scale stress，简称SGS应力），Pa，其表达式为：

在本研究中解决SGS应力的数值计算模型选择 Smagorinsky-Lilly 模型[11-13]。

1.2 计算区域及网格划分

依托常泰长江大桥5号墩沉井基础设计，基于水工模型试验，本文对其进行了数值模拟研究。三维模型及沉井尺寸见图1。

图1（a）展示了数值模拟三维计算模型，计算域在x、y、z方向上分别为11.1 D×8.0 D×0.3 D（D为沉井基础半圆部分的直径）。坐标系为直角坐标系，原点位于沉井中心位置对应的床底面。沉井基础水平截面为圆端形，由2个半圆和1个矩形组成，如图1（b），其中D=57.8 m，L=37.2 m，R=D/2；由于考虑沉井与水流方向的关系，该项目设计沉井纵轴与水流方向成6°夹角。水深H为15.7 m，水流流速Ua分为1.0 m/s、1.5 m/s、2.0 m/s和2.5 m/s 4种情况，沉井底面距离床面高度h分为1.0 m、3.0 m、6.0 m和8.0 m 4种情况，具体工况信息均列于表1。网格划分方式为分块结构化网格，沉井外壁近区采用O形网格，而其他区域则采用H形网格（见图2）。同时对沉井近壁区域、近床面区域进行加密，使靠近壁面的无量纲距离y+，其中：y为第1层网格距壁面的距离，m；ρ为流体密度，kg/m3；u*为摩阻流速，m/s；μ为动力黏滞系数，kg（/m·s），小于1）。

表1 计算工况Table 1 Calculation conditions

图2 沉井近壁面网格（x-y平面）Fig.2 Near-wall grid of the caisson in the x-y plane

1.3 边界条件

在水流入口处采用速度进口边界条件，速度为Ua。在水流出口处设置为压力出口边界条件。自由水面采用刚盖假定，设定为对称边界。其余边界设定为无滑移壁面条件。为了消除初始瞬态流动的影响，模拟在数据统计前先进行52 D/Ua时间单位的初始计算时间，然后进行104 D/Ua时间单位的数据统计与平均计算。在统计平均之前，计算时间足够让流场得到充分发展，使统计量更符合试验情况，从而得到合理且收敛的解。

2 计算结果分析

2.1 算例验证

沉井在下沉过程中，由于作用在沉井上的水动力压强分布不对称形成的压强阻力，以及由于表面切应力造成的表面摩擦阻力，两者之和为绕流阻力，其方向与流体运动方向相反[14]。绕流阻力FD公式如下：

式中：ρ为水的密度，kg/m3；CD为绕流阻力系数；Ua为水流流速，m/s；A为沉井在垂直于水流方向平面上的投影面积，m2。表2为各工况阻力系数计算结果与试验数据的对比，试验在南京水利科学研究院进行。

表2 沉井基础阻力系数CD计算结果与试验数据对比Table 2 Comparison of resistance coefficient CDof the caisson foundation between calculated and experimental data

由表2结果计算可得计算结果与试验值较为一致，误差大致控制在11%以内，这个精度可以满足工程应用。当沉井底面距床面越近，阻力逐渐增大，而阻力系数则逐渐减小。这是因为沉井下沉过程中阻水面积逐渐增大，阻水面积增大的速度快于阻力增大的速度。总而言之，本次模拟计算结果得到了验证，同时也证实了数值方法的可靠性。

2.2 时均流场特性

分析时均流场特性能够有效反应沉井周围水流流动变化规律，对指导施工具有重要意义。图3表示来流流速为1.5 m/s时不同工况水平截面时均流速分布。由图3可见，沉井吃水深度的增加会导致沉井横向两侧流速增大，背水侧流速减小，这说明沉井吃水深度的加深会增强沉井对水流的阻力。在本研究中，由于沉井是其纵轴与水流方向存在6°偏斜角的设计，所以此设计下的流场分布特性也具有明显的非对称特点。例如沉井背水侧的回流区域，由于偏斜角的存在，负y轴后方的背水侧回流区域略大于正y轴后方的背水侧回流区。

图3 来流速度为1.5 m/s时各工况水平截面速度分布Fig.3 Velocity distributions of horizontal plane under different working conditions when the inflow velocity is 1.5 m/s

2.3 相干结构

为了更直观地观察相干结构，本研究采用了一种漩涡判别法则——Q准则。对于不可压缩流动，Q 定义为 Q=（||Ω||2-||S||2）/2，其中，||Ω||和||S||分别代表速度梯度张量的对称和反对称部分[15]。图4给出了沉井在不同吃水深度下尾流的相干结构图。从图4中可以看出，沉井对水流的阻滞作用在相干结构上表现为尾涡，尾涡与圆柱绕流产生的涡结构较为类似，也可以称之为“卡门涡街”。尾部的回水区也可称为漩涡生成区。在漩涡生成区，相干结构尺度较小；随着水流向下游发展，相干结构尺度逐渐拉长、增大。

图4 基于Q准则的速度场的相干结构等值面图（Q=3.4×10-5s-2）Fig.4 Coherent structure isosurface of the velocity field based on Q criterion(Q=3.4×10-5s-2)

由图4可见，漩涡生成区的长度与沉井吃水深度有关，吃水深度的增加导致沉井漩涡生成区长度增长，相干结构尺度增大，涡街的摆幅也相应增大。

3 结语

依托常泰长江大桥5号墩沉井基础设计，采用大涡模拟方法，结合水工模型试验数据系统地研究了沉井下沉过程中的水动力特性，主要研究结果如下：

1）模拟计算的沉井侧壁阻力系数与试验结果吻合度较好，验证了数值模拟方法的可靠性；阻力随着沉井吃水深度的增加而增大，阻力系数则反之。

2）基于时均流场，研究了沉井下沉过程中的流场特性，对各工况水平截面的速度分布进行了分析，结果表明沉井吃水深度的增加会导致沉井横向两侧流速增大，背水侧流速减小。

3）沉井对水流的阻滞作用在相干结构上表现为尾涡，尾涡形状类似于卡门涡街；沉井吃水深度的增加导致沉井漩涡生成区长度增长，相干结构尺度增大，涡街的摆幅也相应增大。