有限流道内低雷诺数二维圆柱绕流数值模拟

桑文慧，孙志强，周孑民

(中南大学 能源科学与工程学院，湖南 长沙，410083)

摘要：为揭示尾迹流场特性及其演化规律，采用多块辐射型网格，以计算流体力学软件FLUENT为平台，对有限流道内低雷诺数二维圆柱绕流进行数值模拟，讨论压力-速度耦合算法、压力离散格式和动量方程离散格式对模拟结果的影响。研究结果表明：动量方程离散格式对数值模拟误差的形成起主要作用，较佳数值模拟算法为SIMPLEC压力-速度耦合、二阶压力离散和QUICK动量方程离散的组合。尾迹流场中静压在沿轴向2.6~12倍流道宽度范围内波动较显著；并且随着流速增大，轴向压降最大点越靠近圆柱，径向静压幅度变化越大。

关键词：圆柱绕流；有限流道；压力-速度耦合；压力离散；动量方程离散

中图分类号：O357. 1         文献标志码：A         文章编号：1672-7207(2012)03-1166-05

Numerical simulation of two-dimensional flow around a circular cylinder at low Reynolds numbers in finite channel

SANG Wen-hui, SUN Zhi-qiang, ZHOU Jie-min

(School of Energy Science and Engineering, Central South University, Changsha 410083, China)

Abstract: To reveal the characteristics and the evolution of cylinder wakes, two-dimensional flow around a circular cylinder at low Reynolds numbers in finite channel was numerically simulated on the computational fluid dynamics code FLUENT using the multi-block radial girds. The effects of pressure-velocity coupling, pressure discretization, and momentum equation discretization on the computation accuracy were discussed. The results show that the discretization scheme of the momentum equation contributes principally to the formation of errors. The combination of SIMPLEC, second-order pressure discretization, and QUICK is found to be more appropriate for the simulation of the two-dimensional flow around a circular cylinder at low Reynolds numbers in a finite channel. It is also found that the static pressure fluctuates markedly along the channel axis within 2.6–12 times of the channel’s width. With the increase of the flow velocity, the position of the maximum pressure drop moves towards the circular cylinder, and the magnitude of the radial static pressure increases.

Key words: flow around circular cylinder; finite channel; pressure-velocity coupling; pressure discretization; momentum equation discretization

钝体绕流，即流体横向流过非流线型物体，广泛存在于自然界和工业生产之中。钝体绕流涉及流动分离、旋涡脱落及其相互作用等诸多问题，尾迹流场极为复杂，不少流动机理尚未完全认识清楚^[1-2]，因而对钝体绕流的研究具有十分重要的学术价值。同时，对于存在绕流的工业设备，流体还可能诱发钝体产生振动，使钝体疲劳损伤导致设备故障^[3]，因此，研究钝体绕流对提高现代工业设备长期运行的可靠性与安全性也具有特别重要的现实意义。钝体绕流是强烈的非线性时变湍流，难以用解析的方法求得其流场分布，所以，至今人们对钝体绕流的认识几乎全部依赖于经验和实验。计算机数值模拟技术特别是计算流体力学的飞速发展为钝体绕流研究提供了一种有效的新方法。在模拟对象方面，Kahawita等^[4]模拟了梯形钝体绕流，讨论了钝体形状对旋涡脱落的影响；Loc等^[5]模拟了圆柱绕流中的二次涡结构；Blackburn等^[6]模拟了二维振荡圆柱绕流的尾迹结构；Yu等^[7]模拟了多孔方柱绕流，讨论了雷诺数和达西数对旋涡尺寸和位置的影响；孙志强等^[8]通过对涡街流量计内钝体绕流的模拟，揭示了在钝体后部一定范围内壁面压力分布规律。在模拟方法方面，Breuer^[9]对雷诺数为140 000的圆柱绕流做了大涡模拟；Yang等^[10]采用动网格和Galerkin有限元法模拟了有限流道内振荡矩形柱绕流的尾迹结构；黄钰期等^[11]使用分块耦合求解方法模拟了圆柱绕流；王亚玲等^[12]采用SIMPLE算法和QUICK格式对圆柱绕流进行了计算，发现三维模拟的升力和阻力系数均小于二维模拟；Cao等^[13]对剪切流的圆柱绕流进行了直接数值模拟和大涡模拟，分析了雷诺数的影响；Rocchi等^[14]比较了圆柱绕流的实验与数值模拟结果。上述研究丰富了人们对钝体绕流的认识。工业上钝体绕流常发生在有限流道内，并且大多数钝体的长度远大于其宽度，因此，研究有限流道内二维钝体绕流具有重要的实际意义。然而，目前针对有限流道内钝体绕流的研究仍比较少。本文作者以计算流体力学软件FLUENT为平台，对有限流道内低雷诺数二维圆柱绕流开展数值模拟，着重讨论压力-速度耦合算法、压力离散格式和动量方程离散格式对模拟结果的影响，寻找较佳的数值格式组合，并根据模拟结果对圆柱绕流流场进行分析。

1  计算模型

当雷诺数较小时，流体紧贴圆柱体表面流动，尾迹中无旋涡；随着雷诺数增大，尾迹中逐渐形成周期的旋涡脱落现象。圆柱绕流过程受到流体流动状况和钝体结构参数等因素的影响。本文研究的二维圆柱绕流计算模型如图1所示。流体沿x方向流经流道内的圆柱钝体，圆柱直径d=14 mm，流道宽度D=50 mm，流道长20D，圆柱中心距流道入口和出口分别为3D和17D。由于流道壁面压力便于实验测量与验证，故设置P₁(160,1)和P₂(160,49)为监视点。流动介质为水，密度为998.2 kg/m³，动力黏度为1.003 g/(m·s)。

图1  计算模型

Fig.1  Calculation model

考虑到圆柱绕流流场的特点，采用多块结构化网格，将计算区域划分为圆柱周围区①、圆柱前区②、近壁面区③和尾迹区④等4个子区域，如图2所示。由于旋涡的产生和脱落发生在圆柱表面，故边长为0.8D的正方形圆柱周围区①采用正交性较好的辐射型网格。圆柱前和靠近壁面处的流场较平稳，为节省计算时间，圆柱前区②采用较稀疏的四边形网格；近壁面区③也采用四边形网格，但较圆柱前区②稍密集。在圆柱尾迹中旋涡强度沿流动方向由强变弱，故尾迹区④采用由密渐疏型网格。流道内壁和圆柱表面网格设置为边界层，并满足无滑移条件。计算区域内的网格总数为44 620，最小和最大网格的边长分别为0.1 mm和18.4 mm。各求解变量收敛标准设置为残差值小于10^-8。为了使每个旋涡脱落周期内的采样点不小于50，时间步长的设置随入口速度不同而相应改变。流道入口设置为速度入口，出口设置为压力出口。模拟的雷诺数为200~600，在该范围内流动为层流，故采用二维非定常不可压缩N-S方程组。

图2  多块辐射型网格

Fig.2  Multi-block radial gird

2  数值格式

FLUENT将控制方程转换为代数方程来求解，通过对每一个控制体积分，得到基于控制体的守恒离散方程，故合理选择离散格式至关重要。此外，由于压力是圆柱绕流求解的重要参数，数值计算中需要根据动量方程和连续性方程推导出压力方程，因此，需合理选择压力-速度耦合算法。

压力离散格式主要有Standard、二阶离散、PRESTO、线性和体积力加权格式。由于有限流道内低雷诺数圆柱绕流控制体中心的压力梯度变化不大，且计算时采用了四边形网格，因此，Standard、二阶离散和PRESTO格式等较为合适。动量方程离散格式主要有一阶迎风、二阶迎风、幂率和QUICK格式，根据圆柱绕流的特点，一阶迎风、二阶迎风和QUICK格式较为适合。压力-速度耦合算法主要有SIMPLE，SIMPLEC和PISO算法。由于进出流道的流体质量守恒，因此，3种方法均比较合适。但由于控制方程速度校正项的不同，使得三者在计算速度上存在一定的差异。

综上所述，适于有限流道内低雷诺数圆柱绕流的压力-速度耦合算法、压力离散格式和动量方程离散格式各有3种，通过组合可以有27种求解组合。本文首先确定较优的组合方法，实现对有限流道内低雷诺数二维圆柱绕流计算精度和速度俱佳的模拟。为了便于后续描述，不同数值算法采用如表1所示的代码表示。

3  斯特劳哈尔数

斯特劳哈尔数St是用于表征旋涡脱落特性的最重要的相似准则数之一，可视为旋涡脱落的无量纲频率。因此，本文以斯特劳哈尔数来判断数值模拟结果的准确性。对于有限流道内圆柱绕流，St的计算式为：

               (1)

式中：m为圆柱两侧弓形流通面积与流道横截面积之比；d为圆柱直径；f为旋涡脱落频率；U为来流速度。本文中m = 0.648。

人们对圆柱绕流斯特劳哈尔数St与雷诺数Re的关系做了大量研究。目前比较公认的St与Re经验关系式为^[15]：

        (2)

该式在47＜Re＜1 000时误差为0.000 5。

本文采用式(2)作为斯特劳哈尔数参考值的计算式，通过比较斯特劳哈尔数的仿真值与式(2)计算所得的参考值之间的相对误差来检验模拟结果的准确性。

表1  数值算法代码

Table 1  Code of numerical algorithms

4  计算结果与讨论

4.1  算法选择

本文模拟了表1所列的不同压力-速度耦合算法、压力离散格式和动量方程离散格式时有限流道内二维圆柱绕流流场。通过FLUENT求解控制方程，得到监视点P₁(160, 1)和P₂(160, 49)的静压，再对静压进行傅里叶分析，求得旋涡脱落频率，最后由式(1)和(2)计算得到斯特劳哈尔数的仿真值和参考值。表2给出了不同数值格式组合下斯特劳哈尔数仿真值与参考值的平均相对误差。

表2  不同数值格式组合下斯特劳哈尔数的平均相对误差

Table 2  Average relative errors of Strouhal numbers under different numerical scheme combinations

从表2可以看出：A₁，A₂和A₃在各种精度的结果中出现，三者对计算结果无影响，表明模拟时可在SIMPLE，SIMPLEC和PISO 3种压力-速度耦合算法中任选其一。B₁和B₃在Ⅰ，Ⅲ和Ⅴ类中出现，表明Standard和二阶压力离散也可任意选择。分别比较Ⅰ与Ⅱ，Ⅲ与Ⅳ，Ⅴ与Ⅵ可以看出：在A和C相同时，Standard(B₁)和二阶压力离散(B₃)的计算精度高于PRESTO(B₂)，这是由于PRESTO格式不太适用于低速流动，不同压力离散格式引入的计算误差相差不大。在A和B相同时，由于计算的是涡流且采用了四边形网格，因此QUICK(C₃)的计算精度最高。从Ⅴ和Ⅵ可知，较大计算误差主要是由一阶迎风(C₁)所引入。动量方程离散格式对于数值模拟误差的形成起主要作用，压力离散格式的影响次之，这是由于二维圆柱绕流模型的求解过程主要是二维N-S动量方程的求解。

根据以上分析，模拟精度最高的数值格式组合有6种。以下通过比较旋涡脱落稳定时间，从计算速度方面进一步确定较佳的数值格式组合。当监视点静压的某一峰值与前后两峰值相差小于0.1%时，认为旋涡脱落达到稳定，故将从开始计算到此峰值的前一峰值的时间间隔作为旋涡脱落稳定时间。表3给出了6种数值格式组合在不同雷诺数下的旋涡脱落稳定时间。比较各个雷诺数下的结果可以看出：压力-速度耦合算法中SIMPLEC(A₂)稳定时间最短，这是由于SIMPLEC与SIMPLE(A₁)速度校正项不同导致的计算速度不同，而PISO(A₃)在每个迭代中要花费较多的CPU时间，虽然能减少到达收敛所需要的迭代步数，但对于有限流道内圆柱绕流模拟，这一优势并不明显。比较同一雷诺数下不同数值格式组合的稳定时间可见，压力离散格式中二阶(B₃)基本上都优于Standard(B₁)，这是因为二阶离散格式适合于流过曲面边界的计算。

表3  A₁B₁C₃，A₂B₁C₃，A₃B₁C₃，A₁B₃C₃，A₂B₃C₃和A₃B₃C₃的稳定时间

Table 3  Convergence times of A₁B₁C₃, A₂B₁C₃, A₃B₁C₃, A₁B₃C₃, A₂B₃C₃ and A₃B₃C₃

综合考虑模拟精度和计算速度，对于有限流道内低雷诺数二维圆柱绕流的计算，建议采用SIMPLEC压力-速度耦合，二阶压力离散和QUICK动量方程离散格式的组合。

4.2  流场分析

采用SIMPLEC压力-速度耦合、二阶压力离散和QUICK动量方程离散格式，计算不同雷诺数下圆柱绕流的涡量，其等值线分布如图3所示。各个雷诺数下圆柱尾迹中都出现了规则的旋涡脱落，在向下游运动的同时，旋涡强度不断减弱。随着雷诺数的增大，相同位置处的旋涡强度逐渐增大。

图4所示为流道内y=0.02D直线上不同雷诺数下的轴向静压分布，圆柱中心位于x/D=3。不同雷诺数下的轴向静压变化趋势基本相同。当x/D＜1.6时，流动几乎不受圆柱影响，静压保持为常数；当1.6＜x/D＜2.6时，静压逐渐增大；当x/D＞2.6以后，随着圆柱边界层分离，静压出现波动，并且随流速增加，静压变化幅度增大，在x/D=3.7附近压降达到最大。流速越大，轴向压降最大点越靠近圆柱。在2.6＜x/D＜12范围内静压波动较显著，随着旋涡向下游运动，静压波动性逐渐减弱，直至稳定。

图3  不同雷诺数下的涡量等值线图

Fig.3  Contours of vorticity at different Reynolds numbers

图4  轴向静压分布

Fig.4  Axial static pressure distribution

图5所示为流道内x=3.6D直线上不同雷诺数下的径向静压分布，圆柱中心位于y/D=0.5。不同雷诺数下的径向静压变化趋势也基本相同。流速越大，径向静压幅度变化越大，这是因为流速增大时圆柱体后生成的旋涡强度越大，因而引起的尾迹波动也越显著。对于雷诺数为200，300和400的情形，静压幅度变化较小，压降最大值出现在圆柱后部上侧。对于雷诺数为500和600的情形，静压幅度变化较大，最大压降位于圆柱后部下侧y/D=0.47附近。

图5  径向静压分布

Fig.5  Radial static pressure distribution

5  结论

(1) 采用多块辐射型网格可较好地模拟圆柱绕流尾迹旋涡的产生与脱落。

(2) 动量方程离散格式对数值模拟误差的形成起主要作用，压力离散格式的影响次之。模拟精度和计算速度俱佳的有限流道内低雷诺数二维圆柱绕流数值模拟算法为SIMPLEC算法、二阶压力离散格式和QUICK动量方程离散格式的组合。

(3) 流速越大，轴向压降最大点越靠近圆柱，径向静压幅度变化越大。在2.6＜x/D＜12范围内静压波动较显著。

参考文献：

[1] Cai J S, Chng T L. On vortex shedding from bluff bodies with base cavities[J]. Physics of Fluids, 2009, 21(3): 034109.

[2] Williamson C H K. Vortex dynamics in the cylinder wake[J]. Annual Review of Fluid Mechanics, 1996, 28(1): 477-539.

[3] Williamson C H K, Govardhan R. Vortex-induced vibrations[J]. Annual Review of Fluid Mechanics, 2004, 36(1): 413-455.

[4] Kahawita R, Wang P. Numerical simulation of the wake flow behind trapezoidal bluff bodies[J]. Computers and Fluids, 2002, 31(1): 99-112.

[5] Loc T P, Bouard R. Numerical solution of the early stage of unsteady viscous flow around a circular cylinder: A comparison with experimental visualization and measurements[J]. Journal of Fluid Mechanics, 1985, 160: 93-117.

[6] Blackburn H M, Henderson R D. A study of two-dimensional ?ow past an oscillating cylinder[J]. Journal of Fluid Mechanics, 1999, 385: 255-286.

[7] Yu P, Zeng Y, Lee T S, et al. Wake structure for ?ow past and through a porous square cylinder[J]. International Journal of Heat and Fluid Flow, 2010, 31(2): 141-153.

[8] 孙志强, 张宏建, 周孑民. 涡街流量计内壁面压力分布的数值模拟[J]. 传感器与微系统, 2008, 27(10): 70-75.
SUN Zhi-qiang, ZHANG Hong-jian, ZHOU Jie-min. Numerical simulation of wall pressure distribution in vortex flowmeter[J]. Transducer and Microsystem Technologies, 2008, 27(10): 70-75.

[9] Breuer M. A challenging test case for large eddy simulation: high Reynolds number circular cylinder flow[J]. International Journal of Heat and Fluid Flow, 2000, 21(5): 648-654.

[10] Yang S J, Chang T R, Fu W S. Numerical simulation of flow structures around an oscillating rectangular cylinder in a channel flow[J]. Computational Mechanics, 2005, 35(5): 342-351.

[11] 黄钰期, 邓见, 任安禄. 黏性非定常圆柱绕流的升阻力研究[J]. 浙江大学学报: 工学版, 2003, 37(5): 596-601.
HUANG Yu-qi, DENG Jian, REN An-lu. Research on lift and drag in unsteady viscous flow around circular cylinders[J]. Journal of Zhejiang University: Engineering Science, 2003, 37(5): 596-601.

[12] 王亚玲, 刘应中, 缪国平. 圆柱绕流的三维数值模拟[J]. 上海交通大学学报, 2001, 35(10): 118-122.
WANG Ya-ling, LIU Ying-zhong, MIU Guo-ping. Three-dimensional numerical simulation of viscous flow around circular cylinder[J]. Journal of Shanghai Jiaotong University, 2001, 35(10): 118-122.

[13] Cao S Y, Ozono S, Tamura Y, et al. Numerical simulation of Reynolds number effects on velocity shear flow around a circular cylinder[J]. Journal of Fluids and Structures, 2010, 26(5): 685-702.

[14] Rocchi D, Zasso A. Vortex shedding from a circular cylinder in a smooth and wired con?guration: Comparison between 3D LES simulation and experimental analysis[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2002, 90(4/5): 475-489.

[15] Henderson R. Nonlinear dynamics and pattern formation in turbulent wake transition[J]. Journal of Fluid Mechanics, 1997, 352(1): 65-112.

(编辑 陈爱华)

收稿日期：2011-04-08；修回日期：2011-06-26

基金项目：国家自然科学基金资助项目(51006125)；中国博士后科学基金资助项目(200801346, 20100471227)

通信作者：孙志强(1980-)，男，河南武陟人，博士，副教授，从事复杂流体流动诊断与控制、新能源与高效节能研究；电话：0731-88879863；E-mail: zqsun@csu.edu.cn