雷诺数3 900时三维带齿圆柱流场特征

赵 通, 林海花, 孙承猛, 焦 波, 孙洪源

(山东交通学院 威海校区 船舶与港口工程学院, 山东 威海 264200)

0 引 言

一直以来，圆柱绕流问题是流体力学中的经典问题，其揭示涡的产生、脱落、消亡等复杂的流动现象。国内外学者均对此展开一系列研究。从二维流动开始：DOU等[1]在低雷诺数Re下研究二维平行壁面的圆柱绕流问题，揭示卡门涡街机理，发现在定常流动条件下，尾流区域的涡并不影响卡门涡街的产生；SINGHA等[2]对二维圆柱绕流流场中靠近壁面的旋涡脱落规律进行研究，发现在非定常旋涡脱落状态下，尾流模式与卡门涡街相反。在三维圆柱绕流研究上：KRAVCHENKO等[3]在Re为3 900时对三维圆柱绕流展开数值分析，采用大涡模拟方法，并采用基于B样条的高阶精确数值方法进行求解，验证了该方法在研究尾流中的准确性，但未说明圆柱绕流的三维效应；王亚玲等[4]用CFX-4研究当Re=1 000、10 000时圆柱绕流的三维效应，发现相对于二维圆柱绕流，当Re较小时三维效应并不明显，当Re较大时圆柱绕流有明显的三维效应，沿柱体不同截面的升阻力系数均不相同，与二维圆柱绕流结果也不相同；WORNOM等[5]用大涡模拟当Re=3 900、10 000、20 000时的流动特性，所得结果与以往试验结果有较好的一致性，发现网格密度对高Re的数值模拟结果有显著的影响；乔永亮等[6]用大涡模拟对Re=3 900时三维圆柱进行数值模拟，验证了该数值方法的正确，也发现在沿圆柱轴向上的三维效应；端木玉等[7]在Re=3 900时采用大涡模拟对三维圆柱流场进行分析，发现三维圆柱绕流有着很强的三维和湍流效应，并在靠近圆柱尾流区观察到U形速度分布。在有限长圆柱上：DOBROSELSKY[8]通过试验研究不同高Re下圆柱绕流流线变化，并采用粒子图像测速(Particle Image Velocimetry，PIV)法观察涡流长度和速度的变化，发现随着Re不断增大，最大回流速度、回流长度、涡流之间的距离都在减小；GAO等[9]对有限长圆柱进行数值模拟，观察其流场特性如压力系数等，并将其与无限长圆柱绕流进行比较，发现其与无限长圆柱绕流有很大差别，譬如有限长圆柱自由端生成的旋涡与尾流相互作用形成4个螺旋旋涡，尾流形状也受到长宽比的影响等，流场特性更复杂。

尽管对于在较高Re下的三维圆柱绕流已有一定的研究，但大多是围绕规则圆柱体展开的研究。在实际工程中，海洋结构物不全是规则圆柱体，还存在方形、椭圆形等更复杂的钝体结构。对此，丁九峰等[10]采用PIV技术研究凹坑表面圆柱尾流特性，发现凹坑结构削弱旋涡脱落过程。SOHANKAR等[11]采用层流模型研究方形柱体在Re=40～200时的流动特性，发现在不同入射角度(0°～45°)下，涡的分离情况不同，流场也有很大的差别。ZHU等[12]研究近壁面对椭圆周围流动的影响，发现在稳态流动过程中椭圆绕流的尾流由2个不对称涡流组成，且间隙比的降低抑制其下部的椭圆涡流脱落。李寿英等[13]在Re=8.28×103、8.28×104时对斜、直圆柱采用CFX层流模型直接对N-S方程进行离散求解，对平均压力系数等进行分析后发现，相对于直圆柱，斜圆柱的平均阻力系数更小。ZHAO等[14]研究不同Re下倾斜三维圆柱体的旋涡脱落，观察倾角为0°和45°时的流场，发现当倾角为45°、Re=700时，其均方根升力系数接近来流角为0°时的均方根升力系数。对于其他不规则三维钝体：XU等[15]在Re=3 900时研究波浪形三维柱体绕流，与规则圆柱相比，其阻力系数小，脉动力抑制效果明显，且其湍流波动出现明显下降；LUO等[16]在Re=3 900时对双曲圆柱进行仿真模拟，并对比其与圆柱绕流的特点，发现双曲圆柱涡脱落尺度比圆柱涡脱落尺度小，升阻力系数比圆柱升阻力系数小。

由此可见，在钝体绕流研究中针对圆柱绕流的研究较多，对于非圆柱形钝体也有一定的研究。然而对于自升式平台中带齿圆柱结构的绕流研究较少。在带齿圆柱数值模拟上，林海花等[17]研究不同Re和不同来流角度下的流场，发现二维带齿圆柱绕流流场与二维圆柱绕流流场有较大差别，当来流方向夹角不同时，尾流流场差异也较大。与圆柱相比，对于带齿圆柱的研究，因其结构复杂性，其流场更复杂。基于此，本文对带齿圆柱在Re为3 900时进行三维绕流分析，观察其尾流特征、瞬时速度和时均速度等特点。

1 研究方法

1.1 数学模型

带齿圆柱流场的流体为不可压缩的液体水，其连续性方程为

(1)

(2)

(3)

式中：μt为湍流黏度，取决于流动状态；k为湍动能；δij为克罗内克符号。

1.2 数值方法

湍流数值模拟最重要的就是对湍流黏度μt进行求解。采用k-ωSST湍流模型，k和ω的输运方程为

(4)

(5)

式(4)和式(5)中具体参数见参考文献[18]。

与LES和DNS相比，k-ωSST湍流模型不会对涡流黏度进行过度预测。根据以往的研究，当Re为3 900时研究绕流现象可很好地捕捉流动现象，因此当Re为3 900时对三维带齿圆柱进行数值模拟，分析该条件下带齿圆柱的流场特征。

2 分析模型

2.1 分析模型与网格划分

以最大作业水深为121.92 m(400英尺，1英尺≈0.304 8 m)[19]的自升式平台桩腿为例，在有限的计算资源条件下，为节省计算时间，取完全浸没在海水里的一段带齿圆柱，并对分析模型进行缩尺，缩尺比为1∶5。所创建的分析模型如图1所示。

图1 带齿圆柱模型

在图1中，带齿圆柱由2块半圆板和1块齿条板组合而成，中间长方形带斜纹部分为齿条板，两侧带斜纹部分为半圆板，齿条板两侧均有齿的存在。三维分析模型中沿圆柱展向共分布10对齿，令每对齿齿端之间的垂直距离为D，齿距为0.4D，齿高(齿顶到齿根的距离)为0.2D。

计算域如图2所示，坐标原点取在带齿圆柱底部的中心点位置，箭头方向是来流方向，以带齿圆柱特征距离D为特征尺寸，上游面距离带齿圆柱中心的距离为10D，下游面距离带齿圆柱中心的距离为20D，两侧面距离带齿圆柱中心的距离为10D，为体现绕流的三维效应，计算域中沿圆柱展向高度的尺寸为4D[20]，即计算域为30D×20D×4D。

图2 流场建模范围及边界条件

流场计算域的网格划分如图3所示。

图3 网格划分示例

计算域采用Fluent Meshing划分网格，采用正六面体进行网格划分，满足第1层网格高度Y+≤1，这里选取第1层网格高度为0.003 18D。为更好地捕捉流场的流动特性，对带齿圆柱周围网格进行加密：图3(a)为流场域整体网格划分图，结构附近网格较密，沿着结构表面法线方向向外网格尺寸逐渐变大；图3(b)为带齿圆柱附近边界层网格划分图。

表1 网格无关性验证

2.2 初始条件和边界条件

计算域入口设置为速度入口(Velocity Inlet)，定常速度U0=0.039 18 m/s；出口设置为压力出口(Pressure Outlet)，表压为0 Pa；上下左右等4个面均设置为对称面(Symmetry)，带齿圆柱表面设置为壁面(wall)。

流场的初始条件如下：流场无穷远处速度为0.039 18 m/s，流体设置为水，密度为998.2 kg/m3，动力黏性系数μ=0.001 003 kg/(m·s-1)，Re≈3 900，来流方向垂直于一对齿的方向。

2.3 数值方法验证

表2 数值方法验证

3 带齿圆柱绕流分析结果

3.1 带齿圆柱受力分析

图4为带齿圆柱升阻力系数的时历变化曲线，点线代表升力系数Cl，实线代表阻力系数Cd。由升阻力计算公式可知，升阻力系数变化代表升阻力的变化。由图4可知，带齿圆柱表面升力在不断地周期性变化，振动幅值也在0值上下不断变化，但计算稳定后，振动幅值变化不大，这符合湍流状态下流动的不规则性。

图4 带齿圆柱升阻力系数

由图4还可知：当Re≈3 900且来流方向垂直于一对齿的方向时，带齿圆柱的升力系数和阻力系数在一定范围内上下波动，这也符合亚临界状态下流动的特点，随着流动时间的增加升阻力系数变化趋缓；计算稳定后(200 s之后)，升力系数幅值稳定在0.62附近，阻力系数幅值稳定在1.59附近，升力的振荡变化周期为阻力振荡变化周期的2倍。

进一步对计算稳定后的升力进行功率谱密度分析，得到如图5所示的升力功率谱密度曲线，从而可以得到涡泄频率fs。

图5 升力功率谱密度

由图5可知，所计算的带齿圆柱在Re为3 900时的涡泄频率在0.08 Hz附近。

为更深入地研究带齿圆柱瞬时升力和阻力的脉动变化，取图4中瞬时时间为278.05 s时刻，在带齿圆柱高度为3.884 5D处，取其外表面观察其压力系数Cp的分布规律。该位置处的截面如图6(a)所示，圆柱截面左侧为来流方向。经流场分析后，得到其相对应的压力系数Cp曲线，如图6(b)所示。

图6 带齿圆柱受力

由图6可知：在垂直于来流方向上，左侧半圆板驻点处压力系数最大，随后开始降低；当流动位于齿条处，压力系数陡然增加，但齿条处的压力系数仍然小于圆柱驻点处的压力系数；当齿条处压力系数陡然增加后，沿着尺高方向，压力系数又迅速下降，流动到达齿顶处，压力系数已下降为负值，表明在该位置产生边界层分离并形成速度回流区；在齿条后方的压力系数均为负值，与前端正压力系数一起形成作用于带齿圆柱上较大的阻力系数。

3.2 时均流场分析

为方便观察带齿圆柱附近尾流速度分布情况[7]，提取不同特征位置处的速度进行分析。带齿圆柱周围特征位置分布如图7所示，所取的特征位置平面位于z=3.884 5D的xOy平面，y/D=0表示该截面处带齿圆柱中心沿x轴正向的直线。具体的特征位置分别位于x/D=0.58、1.06、1.54、2.02处，即该截面上与y轴平行的4条直线段[3]，所有直线段的长度为6D，且分别关于x轴对称。

图7 带齿圆柱周围特征位置

图8为计算后得到的带齿圆柱中心线y/D=0上的x方向时均流速分布，加粗曲线为平均速度分布，其余线条为不同瞬时速度分布。

图8 带齿圆柱中心线上来流方向速度分布

由图8可知：在靠近带齿圆柱的附近，产生负速度，这主要是因为在此距离内带齿圆柱旋涡脱落引起回流区，与圆柱在此回流区速度相比，带齿圆柱在此回流区形成的形状更趋向于V形；当x/D=2时负向速度消失，这是因为该距离超过其回流区长度，其速度也越来越接近入口速度；在x/D=3后，速度趋于平缓。

图9分别给出x/D=0.58、1.06、1.54、2.02处x方向多个瞬时速度和剖面时均速度，在这里共取30个瞬时速度，时间间隔为0.05 s，瞬时速度曲线用细线表示，时均速度曲线用粗线条表示。

图9 不同x/D下x方向瞬时速度和时均速度

由图9可知：在靠近带齿圆柱壁面附近处，x方向的速度曲线呈现V形，且底部较尖锐，负向速度达最大，这是因为在靠近带齿圆柱的地方，流动十分复杂，且此处的回流区速度也较大，因此产生较为尖锐的V形；在距离稍远处，x流向的速度逐渐趋于平缓，速度曲线逐渐向U形过渡，并最终呈现U形，负向速度也在变小，说明在带齿圆柱稍远处，回流区速度减小；距离带齿圆柱的距离再远一些，速度变化逐渐平缓，并逐渐趋于初始速度；在x/D=2.02处，时均速度分布表现出不对称性，这是因为在此处带齿圆柱对流场的控制减弱。

图10 给出在不同x/D下y向瞬时速度和时均速度分布，取30个瞬时速度，时间间隔为0.05 s。取x/D=1.06、1.54、2.02处分析y方向瞬时速度和时均速度。

图10 不同x/D下y方向瞬时速度和时均速度

由图10可知，y向瞬时速度在平均速度上下脉动：在x/D=1.06和x/D=1.54位置处，y向时均速度关于原点对称，与圆柱绕流流场在此处的平均速度相比，带齿圆柱绕流的平均速度波动更大一些[7]；在x/D=2.02处，带齿圆柱对其尾流速度的控制有减弱的趋势，时均图像也趋于平缓，且不再具有关于原点对称的特点。

3.3 瞬时流场分析

在t=278.05 s时刻，流动趋于稳定，可较好地观察流场特点，取此时瞬时时间，观察瞬时流场特点。图11为带齿圆柱后方的旋涡脱落图，所提取的2个图为该瞬时时刻下2个不同视角下的旋涡脱落图。

图11 带齿圆柱瞬时速度涡量图

由于齿的存在，周围流场变得更复杂，在齿周围也有旋涡脱落，这与圆柱的旋涡脱落有很大区别。同时，由于齿的存在，旋涡脱落地更不规则，在带齿圆柱后方依然存在速度回流区。为观察带齿圆柱流场的三维特性，取该瞬时时刻下4个不同展向位置的平面，4个位置平面如图12所示，平面之间的间隔距离为0.8D，分别为z/D=3.884 5、3.084 5、2.284 5、1.484 5平面(每个平面均通过齿顶端)，观察其流场分布特点。图13即为带齿圆柱在同一时刻下不同位置处的流线图。

图12 带齿圆柱不同位置

图13 不同z/D瞬时速度流线图

由图13可知：流体在流经带齿圆柱时，在齿条位置处发生边界层分离，齿条后方随即产生速度回流区，形成旋涡，并开始脱落，且在脱落位置处其速度也在变小；在同一时刻，不同展向位置处，带齿圆柱流场也有很大的不同，这也体现了带齿圆柱流场的三维特性。

4 结 论

采用雷诺平均k-ωSST湍流模型对Re为3 900时的三维带齿圆柱绕流流场进行分析，得到升阻力系数时历曲线，总结旋涡脱落、瞬时和时均速度分布的特征，得到以下结论：

(1)在同一时刻，不同展向高度下，速度流线图均不相同，表现出带齿圆柱绕流的三维效应。

(2)在带齿圆柱的顺流向和横流向，均产生负向速度，即产生回流区。其瞬时速度始终围绕时均速度上下波动，并在壁面处由近及远分别呈现V、U形等特点。

(3)带齿圆柱的升阻力系数在一定范围内波动，很好地体现了湍流状态下流动的不稳定性，同时观察到旋涡脱落现象，涡泄频率约0.08 Hz。

工程界在评估海洋工程结构的疲劳损伤时，往往只考虑往复运动的波浪水质点的影响。对带齿圆柱的绕流进行分析，升力、阻力的波动都较为明显，尤其是横流向升力的波动，不仅频率较快，而且波动幅值也较大，必将对结构带来明显的疲劳损伤，因此建议在评估此类结构物的疲劳损伤时，将涡激因素考虑进去。