三维有限长圆柱绕流数值模拟

王晓聪，桂洪斌，刘洋

哈尔滨工业大学（威海）船舶与海洋工程学院，山东威海264209

0 引 言

绕流现象广泛存在于自然界中，无论是陆上结构物（如桥梁上的缆索），还是海洋中各类柱状结构（如石油平台的桩腿），都会受到绕流的影响。圆柱绕流一般会产生漩涡脱落并诱发流体的脉动载荷与结构振动［1］，从而发生涡激振动现象。当振动严重时，会对结构物本身造成疲劳等结构破坏［2］，所以对经典的圆柱绕流现象进行研究具有重要意义。截至目前，尽管学术界对于经典的圆柱绕流流体力学问题进行了诸多研究，但对于有限长圆柱（Finite Cylinder，FC）的研究却并不多见。鉴于圆柱上端端面的下洗作用，有限长圆柱的绕流现象与无限长圆柱（Infinite Cylinder，IFC）的绕流现象相比区别较大，其三维特性显著［3］，故非常有必要对其开展研究。

首先，在实验方面，国内外学者主要从圆柱长径比AR、自由端斯托罗哈尔数St、自由端阻力系数的变化规律等方面进行了研究。Park等［4］利用粒子图像测速（Particle Imagine Velocimetry，PIV）技术对雷诺数Re=7.5×103的有限长圆柱进行了实验分析，发现其自由端成对出现梢涡，经分析后发现梢涡是由于端面的下洗流动和圆柱的绕流流动相互作用产生。Benitz等［5］在雷诺数 Re=2.9×103和傅汝德数Fr=0.65的条件下对有限长圆柱同时进行了一系列长径比的实验和数值模拟研究。研究发现，相比于无限长圆柱，有限长圆柱绕流的阻力系数更小，且在AR=2时阻力系数减小得最为明显。此外，当AR＜3时，由于自由端效应，下洗的流动会将卡门涡街完全抑制，如图1所示。

其次，在数值模拟方面，Fröhlich 等［6］使用LESOCC2代码，通过大涡模拟（Large Eddy Simula⁃tion，LES）方法求解不可压缩的N-S方程，并计算了AR=2.5的有限长圆柱。模拟中，分别使用了SM（Smagorinsky model）和 DSM（Dynamic Smago⁃rinsky model）亚网格尺度（Subgrid Scale，SGS）模型。结果发现，DSM对于模拟当前长径比下的有限长圆柱缺陷明显，它会导致计算中涡粘性偏大，致使下壁面边界条件从边界层向对称面转变，而SM模型则能更好地完成仿真模拟。Palau-Salvador等［7］基于浸入边界法，使用大涡模拟对2种有限长圆柱（AR=2.5，5）在 2种雷诺数（Re=4.3×104，2.2×104）工况下进行了实验和仿真分析。结果发现，较短圆柱的交替漩涡脱落仅发生在接近壁面处，且不稳定。

本文基于前人的研究结果，拟分别对具有相同几何尺寸且在相应计算域下为有限长和无限长的圆柱，选取前人研究较多的Re=3.9×103工况进行模拟，以验证本文所提方法的可行性。然后，分别对有限长圆柱流场的时均特性和瞬时特性进行分析，以模拟出有限长圆柱特有的“马蹄”涡、梢涡及“蘑菇”涡。最后，将得到的数据和流场特征与无限长圆柱的模拟结果进行对比，以了解有限长圆柱流场的三维特性和自由端面以及固定壁面对流场的影响，并期望从流场的涡泄角度来解释有限长圆柱阻力系数损失的原因，为半空间结构声振问题研究提供参考。

1 数学模型

1.1 控制方程

本文将采用大涡模拟的湍流模型进行数值模拟，该方法的基本原理是将漩涡在空间过滤，其中大漩涡采用直接模拟方法，而小漩涡则采用SGS模型来模拟其对大涡的影响。当考虑粘性、不可压缩性的水作为流体介质时，N-S方程可写为

式中：ρ为流体密度；p为压力；ν为流体的运动粘性系数；ui，uj为速度分量；xi，xj为位移分量；t为时间。从湍流瞬态运动方程中将尺度比滤波函数尺度小的涡滤掉，分解出描写大涡流场的运动方程。采用滤波函数处理下的N-S方程及连续性方程，即大涡模拟的控制方程为：

本文采用SM模型，定义亚网格应力为

式中：τkk为各向性部分的SGS应力；δij为Kronecker符号；为 SGS湍流粘度，其中Cs为Smagorinsky常数，在外流条件下通常取分别为沿x，y，z轴方向的网格尺寸，Δ为滤波尺寸，为可解尺度变形率的张量表达式。

1.2 几何模型与网格划分

为对比有限长圆柱绕流的流动特性，本文选取具有相同几何尺寸的有限长及无限长圆柱模型进行计算和分析，如图2所示。图2（a）中，无限长圆柱计算域的展向高度等于圆柱高度；图2（b）中，有限长圆柱计算域的展向高度大于圆柱高度。

上述几何模型采用笛卡尔坐标系作为计算域参考坐标，图中x，y，z分别表示顺流向、横流向和展向，坐标原点位于圆柱地面中心。其中，有限长圆柱计算域尺寸为10D×20D×6D，无限长圆柱计算域尺寸为10D×20D×πD，D为特征尺寸，即圆柱直径）。此圆柱展向高度为πD，圆柱中心距离上游入口为5D，距离下游出口为15D，模型阻塞比为5.23%。由于圆柱自由端面与流体域上表面必须留出足够的空间，以便精确捕捉自由端面附近的流动情况，考虑到自由端上部流场的发展，本文中自由端面距离流场域上表面约为3D。

图3所示为圆柱计算域网格尺寸划分，其中无限长圆柱划分俯视图与有限长圆柱的一致，采用O形网格对圆柱周边网格区进行局部加密。为有效监测有限长圆柱自由端的流动，该区域网格划分也进行了局部加密，以精确捕捉流动特性。

1.3 边界条件和时间步选择

边界条件设置为速度入口、压力出口、壁面、对称面，分别对应几何体的入口面、出口面、侧面、上下壁面。

对于n维问题的显式计算，时间步应满足Courant-Friedirchs-Lewy（CFL）条件，其表达式为

式中：C为库朗数（Courant number），Cmax=1；ui为i方向流速；Δi在此处为最小计算网格尺寸；Δt为计算时间步。

对于本文的计算，取n=1，因此只需满足u(Δt/Δx)≤1即可。此时，u即为顺流向速度U0。设置无因次时间Δt*=ΔtD/U0=5.13×10-4，则计算时间步Δt选为2.5×10-3s，C=0.975，故满足条件。

1.4 模型验证

为有效验证计算的准确性，在同等初始条件下，分别对有限长和无限长圆柱计算域进行计算。在Fluent计算流体力学软件下，设置压力速度耦合求解的SIMPLE格式，差分格式则采用中心差分。以远点流速和圆柱直径为特征参数的雷诺数Re=3.9×103。将计算结果与前人所做实验和数值模拟结果进行了对比，如表1所示。表中，Sim1指文献［8］得到的数值模拟结果，EXP1，EXP2和EXP3分别指文献［9］、［10］和文献［5］的实验结果。

表1 数值计算结果与文献实验结果的对比Table 1 Comparisonof numericalsimulation and experimental results provided by the literatures

表中列出了时均阻力系数d、升力系数均方根值Clrms、漩涡脱落频率相关参数斯托罗哈尔数St。斯托罗哈尔数St计算公式为

式中，f为漩涡脱落主频率。

由表1可知，在同等SGS模型下，无限长圆柱算例IFC1和有限长圆柱算例FC1与实验结果非常接近，尽管无限长圆柱算例IFC3与有限长圆柱算例FC3更加接近，但是相比计算资源消耗的增加，无限长圆柱算例IFC1和有限长圆柱算例FC1是更合理的方案。

为验证文献［6］中关于SGS应力模型的结论，本文针对有限和无限长圆柱，在SM及DSM模型中均计算了相关参数。经对比发现，对于当前工况，SM模型确实能更精确地模拟亚临界区圆柱绕流。在后续计算中，无论是有限长还是无限长圆柱绕流的计算，本文均选用SM亚网格尺度应力模型，且SM模型常数Cs=0.2［11］。同时，本文将计算结果与前人的实验结果进行了对比，发现在IFC1和FC1工况下，得到的结果与实验结果非常接近，且与作为衡量准确性最强的物理量d最接近［6］，故本文后续计算均取IFC1和FC1的工况进行运算。

为进一步验证算例的正确性，分别列出有限长和无限长圆柱沿中心线(y=0，z=πD/2)上的时均顺流向速度ux及圆周压力系数Cp分布，并与文献［12］的计算结果进行对比，如图4所示。其中，压力系数Cp的计算公式为

式中，p，p0分别为监测点压力和无穷远处压力。

由图4可知，本文计算结果与文献［12］结果误差在3%以内，这进一步说明了本文计算工况选取合理。此外，通过对比无限长圆柱的计算结果还可发现，有限长和无限长圆柱的顺流向速度在时变为正值，这说明有限长圆柱的回流区更小，自由端面引起的下洗对回流区的影响导致顺流向阻力系数Cd损失，这与表1的计算结果一致。

2 三维流动特性分析

2.1 时均流场特性分析

图5所示为有限长和无限长圆柱的升、阻力系数稳态时历曲线。由图可知，相比于无限长圆柱，有限长圆柱的时均阻力系数要小23.6%，这说明自由端的存在明显改变了圆柱受力，进一步说明了下洗作用对圆柱绕流流场的干扰。此外，相比于无限长圆柱，有限长圆柱的升力系数Cl的峰值更小，同时也表现得更无规律性。鉴于2种圆柱的实际几何尺寸一致，说明在有限长圆柱计算域下，下洗作用至少局部地扰乱了原有亚临界区表现出的卡门涡街。

为研究圆柱存在的自由端面对流体下洗作用的影响，将对圆柱展向居中处水平截圆，以其前缘点为初始监测点（θ=0°），每隔10°沿主流方向至180°均匀分布圆周监测点，分别计算各监测点处的压力系数时均曲线。为与无限长圆柱进行对比，在同等监测位置（即无限长圆柱的展向中间）处设置了同等数量的监测点，检测结果如图6所示。

由图6可知，有限长和无限长圆柱在中间位置处分离角相同，不同的是对于压力系数Cp，有限长圆柱与无限长圆柱的差值达到了35%。由此可见，对于较小长径比的有限长圆柱，端面下洗流可对流动造成巨大干扰，从而对流场带来极大的改变。此外，由于下洗的作用，圆柱尾流近壁面处流动受到了一定的抑制。

为进一步研究自由端面和固定壁面对圆柱近壁区尾流的影响，在y=0，x=D的垂向直线上分别为有限长和无限长圆柱均匀设置了60余个监测点，用于监测顺流向速度ux，无因次化处理后的结果如图7所示（图中，h为圆柱展向高度）。由图可知，对于无限长圆柱，在x=D处，其顺流向速度ux均为负值。结合对顺流向速度ux沿x方向的分布曲线，说明圆柱的回流区大于D，与前文对压力系数的分析一致。同时，在z≤0.8h的区域内，无论有限长圆柱还是无限长圆柱，ux同时为负，但当z≥0.8h时，有限长圆柱速度迅速跃升，从自由端面直到距离自由端面约1.2D的位置，速度变为无穷远处流速。从本文进一步分析z/h=0和z/h=π时的有限长及无限长圆柱的相对差值可以发现，自由端对顺流向速度ux的影响要大于固定壁面的影响。

2.2 瞬时流场特性分析

由于有限长圆柱的实际流动时，下壁面，即圆柱固定壁面并非光滑壁面，固定壁面的非滑移性质导致产生了“马蹄”涡。为能同时得到梢涡和“马蹄”形漩涡，本文将原有有限长圆柱计算中的上、下面对称条件改为下壁面非滑移，并重新计算得到了有限长圆柱绕流的流场。

为能清晰得到圆柱后漩涡脱落的流动，本文采用Q准则（Q-criterion）定义涡量如下：

式中：Sij为应变张量；Ωij为计算涡量。

在三维笛卡尔坐标系下，简化公式为

式中，u，v，w分别为3个方向的速度分量。

图8所示为以Q准则为基准的涡量等值线分布图。图8（a）为2个时刻自由端面的梢涡涡量等值线图。由图可看出不同时刻的涡量发展，梢涡是成对出现，另外顶面处漩涡分离角约为70°；图8（b）为圆柱展向中间位置处的平面涡量等值线图。由图可看出，此时的漩涡已不再是标准的卡门涡街漩涡分布，上、下漩涡也不再对称。可见自由端的确已对中间处的漩涡发展产生影响，这也印证了本文之前对升、阻力系数分析时得到的结论；图8（c）为固定圆柱壁面处的涡量等值线图，由图可看出，在底面处因同时考虑了底面的非滑移属性，底面漩涡流动非常复杂，但在圆柱前端仍可清晰地看到“马蹄”形漩涡的涡线。值得注意的是，底面处圆柱后因受到的下洗作用影响较小，故保留了一定的、类似于无限长圆柱的漩涡脱落属性。这也间接说明，对于AR=π的圆柱，下洗并未将整个圆柱的流场完全扰乱，自圆柱底面至其展向高度的一半处，卡门涡街仍是流动的主要成分。

图9所示为有限长圆柱自由端面的x方向速度分量云图及其自由端面流线图的局部放大图。由图可知，自由端面自前缘开始产生漩涡，漩涡逐渐发展，并覆盖自由端面，形成“蘑菇”形漩涡，漩涡发展区域约占自由端表面长度的2/3。这说明，从圆柱前缘发生分离的流体在2/3处出现回旋后，其余流体继续向前流动，未发生回贴；在自由端的最后向下，流体继续分为2个部分：一部分向圆柱壁面流去，形成回旋，即梢涡，正是这种梢涡的存在与卡门涡街一起构成了有限长圆柱的主要流动特性；另一部分继续向后流去，但其流动轨迹已不再平直，而是向着偏下的方向流动。

图10所示为圆柱壁面极限流线图。图10（a）中，极限流线交汇成沿展向的某一直线，此直线即为流动绕圆柱壁面的分离点构成的线。由图10（b）所示顶面自由端极限流线图可明显看到2个流线源点，该源点即为“蘑菇”形漩涡形成的源点，正是在此分离的流体分别沿中线向两侧回旋，形成漩涡，这与文献［2］的研究结果一致。

为更清晰地了解圆柱三维涡结构，分别得到有限长和无限长圆柱基于Q准则的涡量等值线图，如图11所示。从图中可清晰地看到，相比于无限长圆柱，有限长圆柱的漩涡结构更复杂，但却相对较弱，其绕流被限制在相对较小的区域内，但有限长圆柱的漩涡在趋近下游区域逐渐变得光滑。从图 11（a）和图11（b）中可以看到，由于存在下洗作用，有限长圆柱的回流区更短，且涡截面形状明显呈“三角形”，这与前文流线图中圆柱后的尾流流场中流线斜向下发展的现象一致；而无限长圆柱则不存在这种现象，圆柱后的流场比较均匀，截面呈标准的“矩形”。正是由于上述原因，造成了有限长圆柱阻力系数Cd的损失。除尾流涡结构外，从图11（c）和图11（d）可以看出，在圆柱柱体与地面交界处，可以清晰地观察到“马蹄”形漩涡，这正是由底面的非滑移属性造成的。

3 结 论

本文基于大涡模拟数值方法，分别对Re=3.9×103工况下的有限长和无限长圆柱进行了数值模拟。通过与文献研究结果的对比，证明了本文所提方法的正确性。在时域范围内，本文分别从有限长圆柱的流场速度、阻力系数、升力系数和瞬态流场云图等方面与无限长圆柱进行了对比分析，并得到如下结论：

1）相对于无限长圆柱，有限长圆柱的流动更为复杂，其阻力系数要小23.6%；而对于圆周压力系数的监测结果表明，有限长圆柱的压力系数小35%，其主要原因是自由端面的下洗流动改变了原有卡门涡街的流动模式，造成其尾流流动减弱。同时，有限长圆柱的升力系数也比无限长圆柱的更小，这说明有限长圆柱对圆柱受力起到了明显的抑制作用。在垂直方向上对有限长和无限长圆柱顺流向速度的分析，得出自由端面对有限长圆柱的影响大于固定壁面。

2）基于Q准则（Q=6）对有限长的流场进行了分析，在展向（z=πD，z=1/2πD，z=0）3个面分别进行了涡量等值线分析，发现自由端对圆柱后缘尾流影响的范围在AR=3的工况下超过了圆柱高度的一半。在圆柱底面，由于存在固定壁面，圆柱底部前缘会出现有限长圆柱特有的“马蹄”形漩涡，该漩涡对于圆柱后缘底部流动同样具有一定影响，导致底部流场同样非常复杂，但卡门涡街相对于顶面受到的影响相对较小。对圆柱中面(y=0)进行顺流向速度云图及流线分析可以发现，顶涡，即“蘑菇”形漩涡的发展区域约占自由端面长度的2/3，观察自由端流线分布可清晰得到2个流线源点。此外，绕过自由端的流体会分为2个部分：一部分回贴向壁面，形成梢涡；另一部分沿斜向下的方向流去。后面的涡量等值面则间接证明了这一点。

3）同样基于Q准则，分别对有限长和无限长圆柱作出涡量等值面，可以看到有限长圆柱的流动比无限长圆柱的弱，回流区较短，下洗流动的存在致使其界面流场并非像无限长圆柱那样呈“矩形”，而是呈“三角形”形状。

参考文献：

［1］王露，李天匀，朱翔，等.浸没边界格子Boltzmann方法的改进及转动圆柱绕流模拟［J］.中国舰船研究，2016，11（6）：97-103.WANG L，LI T Y，ZHU X，et al.Improved immersed boundary lattice Boltzmann method and simulation of flow over rotating cylinder［J］.Chinese Journal of Ship Research，2016，11（6）：97-103（in Chinese）.

［2］NAUDASCHER E，ROCKWELLD.Flow-induced vi⁃brations：an engineering guide：IAHR hydraulic struc⁃tures design manuals 7［M］.Rotterdam：Balkema Pub⁃lishers，1994.

［3］SUMNER D.Flow abovethefreeend of a sur⁃face-mounted finite-height circular cylinder：a review［J］.Journal of Fluids and Structures，2013，43：41-63.

［4］PARK C W，LEE S J.Effects of free-end corner shape on flow structure around a finite cylinder［J］.Journal of Fluids and Structures，2004，19（2）：141-158.

［5］BENITZ M A，CARLSON D W，SEYED-AGHAZA⁃DEH B，et al.CFD simulations and experimental mea⁃surements of flow past free-surface piercing，finite length cylinders with varying aspect ratios［J］.Comput⁃ers&Fluids，2016，136：247-259.

［6］FRÖHLICH J，RODI W.LES of the flow around a cir⁃cular cylinder of finite height［J］.International Journal of Heat and Fluid Flow，2004，25（3）：537-548.

［7］PALAU-SALVADOR G，STOESSER T，FRÖHLICH J，et al.Large eddy simulations and experiments of flow around finite-height cylinders［J］.Flow，Turbu⁃lence and Combustion，2010，84（2）：239-275.

［8］WORNOM S，OUVRARD H，SALVETTI M V，et al.Variational multiscale large-eddy simulations of the flow past a circular cylinder：Reynolds number effects［J］.Computers&Fluids，2011，47（1）：44-50.

［9］KAWAMURA T，HIWADA M，HIBINO T，et al.Flow around a finite circular cylinder on a flat plate：cylin⁃der height greater than turbulent boundary layer thick⁃ness［J］. Bulletin ofJSME，1984， 27（232） ：2142-2151.

［10］OKAMOTO S，SUNABASHIRI Y.Vortex shedding from a circular cylinder of finite length placed on a ground plane［J］.Journal of Fluids Engineering，1992，114（4）：512-521.

［11］WELLER H G，TABOR G，JASAK H，et al.A tenso⁃rial approach to computational continuum mechanics using object-oriented techniques［J］.Computers in Physics，1998，12（6）：620-631.

［12］ZHANG H，YANG J M，XIAO L F，et al.Large-ed⁃dy simulation of the flow past both finite and infinite circular cylinders at Re=3900［J］.Journal of Hydrody⁃namics（Ser.B），2015，27（2）：195-203.