五圆柱体结构群绕流特性及互扰效应研究

2020-04-10涂佳黄张志豪肖毅凡梁经群

海洋工程 2020年1期

涂佳黄，张志豪，曹波, 3，肖毅凡，梁经群

(1. 湘潭大学土木工程与力学学院，湖南湘潭 411105; 2. 岩土力学与工程安全湖南省重点实验室，湖南湘潭 411105; 3. 中车株洲电力机车有限公司，湖南株洲 412000)

自然界与工程界中，柱体结构大多以结构群的形式出现，如高层建筑群、冷却塔群、海洋平台立柱和海洋立管等。随着柱体结构数目与排列方式的改变，多柱体结构群布置形式更加多样化，从而导致其绕流问题也变得更加复杂。

现阶段，一些学者对圆柱体结构群绕流问题开展了一系列的物理试验和数值计算研究，主要集中于三/四圆柱体群工况[1-8]。Lam和Zou[1-2]题组运用数值模拟与试验方法对三/四圆柱体群绕流问题进行了研究，分析了布置方式、间距比与雷诺数等参数对圆柱体群的流体力系数和流场流动模式的影响。同时，将数值结果与物理试验结果进行了对比验证，发现两者吻合较好，验证了数值方法的正确性与可靠性。另外，顾志福和孙天风[3]与Bao等[4]分别采用风洞试验与有限元方法对层流范围内等边排列三圆柱体群绕流问题进行了研究，并分析了来流攻角与间距比参数对圆柱体的流体力系数和流场分布特性的影响。最近，Gao等[5]与Tong等[6]对低雷诺数范围内方形排列四圆柱体绕流问题进行了数值模拟，分析了雷诺数对圆柱体群泄涡模式的影响并将其进行了归纳分类。Zou等[7]运用有限体积法研究低雷诺数下间距比对菱形排布四圆柱体群三维绕流特性影响，并与水槽试验结果对比分析验证其数值方法的正确性，研究发现狭窄间隙流模式到涡流冲击模式的转变发生在间距比为3.0工况。Wang等[8]对Re=8 000时正方形排列四圆柱体群绕流问题进行了物理试验研究，分析了来流攻角和间距比对四圆柱体绕流特性及其互扰效应的影响。

近几年，国内外学者们对阵列式排布柱体群绕流问题开展了一些研究工作[9-13]。Nicolle和Eames[9]采用高精度计算方法对亚临界雷诺数下Re=2 100不同数量单圆柱体组成的固定直径比的阵列圆柱体群绕流问题进行了二维数值计算，分析了空隙流对阵列圆柱体群绕流特性的影响，并观察到三种典型的流场模态。李强[10]对九圆柱体组合结构群绕流特性问题进行了水槽试验和数值计算，发现圆柱体群直径与单个圆柱体直径比和雷诺数对圆柱体群尾流形态影响显著。Ma等[11]对正方形排列九圆柱体绕流问题进行了试验研究，分析间距比和雷诺数对圆柱体的流体力系数及其尾流场特性的影响。与雷诺数相比，尾流场特性和流体力对间距比参数的变化更为敏感。

综上，现阶段对多柱体结构群绕流特性及互扰机理的认识还未成熟。文中对“X”形排列形式的五圆柱体在亚临界雷诺数下三维绕流问题进行了数值计算，主要分析来流攻角和间距比对五圆柱体结构群三维绕流特性的影响，并揭示“X”形排列五圆柱体结构群的互扰效应及其内在机理。

1 数值计算理论

1.1 流体控制方程

基于大涡模拟湍流模型，通过滤波函数将大尺度涡和小尺度涡分离，大尺度涡用N-S方程直接求解，小尺度涡通过亚格子尺度模型建立与大尺度涡的关系并进行模拟。滤波函数为：

(1)

(2)

式中：Δ为计算网格体积，G(x,x′)滤波函数。经滤波分析后，不可压缩流体的Navier-Stokes方程为；

(3)

(4)

采用有限体积法对控制方程进行离散，并采用SIMPLE算法对压力与速度耦合，空间离散压力项和动量项分别为二阶格式和边界中心差分格式，瞬态方程为边界二阶隐式格式。

1.2 无量纲参数

文中相关参数计算式如下：阻力系数CD=2FD/(ρU2HD)，升力系数CL=2FL/(ρU2HD)，压力系数CP=2(P-P∞)/(ρU2)，斯特劳哈尔数Sr=fsD/U。其中，FD与FL分别为阻力与升力；H和D分别为柱体轴向高度和直径，P为静止压力，P∞为无穷远处压力且P∞=0，fs为漩涡脱落的频率。

2 计算模型和网格划分

计算域的边界条件设置：速度入口边界，即ux=U=1.0 m/s，uy=uz=0；压力出口边界，即p=0；上下壁面和前后壁面边界为对称边界条件；圆柱壁面边界为无滑移壁面条件，即ux=uy=0。

文中的计算流域采用结构化六面体网格来划分，网格总数约6.2×106，柱体径向长度上的节点数为42，采用圆柱为中心范围内的十字型加密区间采用Otype型，加密区域为5D×5D，加密层数为8层。对靠近圆柱体壁面的网格进行局部加密，这样可以更好地刻画边界层处的流场。同时圆柱体后方的尾流区域也进行局部加密处理，如图2所示。另外，从图2可知网格加密区未产生网格畸变。

图1 计算模型与来流攻角示意Fig. 1 Schematic diagram of the computational model and incident angle

图2 网格划分示意Fig. 2 Schematic diagram of the mesh in different cases

3 计算参数选取及验证

3.1 计算参数分析

本节以L/D=4.0工况下正方形排列四圆柱体结构群(α=0°)作为研究对象进行计算参数研究，分析时间步长、圆周节点数和壁面网格尺度对四圆柱体的阻力系数平均值和斯特劳哈尔数的影响，并引入文献[12]的结果进行对比。文中时间步长符合克朗数(Cr，且取Cr=1)计算公式要求：Δt=(Cr) Δx/U，通过计算可知时间步长取0.01可满足要求。通过对比发现，当Δt较大时会导致计算结果发生偏差，而较小的Δt会消耗较大的计算资源，其中平均阻力系数最大误差约为3%，Sr最大误差约为7%，详细见表1。综合计算资源的消耗、频域采样率和计算结果有效性，本算例的网格与时间步长均选用Case4方案。同时边界层第一层高度均满足y+≤1.0的条件。本文结果与文献[12]结果存在一定细小差异，这主要是由于文献[12]采用了二维模型对问题进行求解，其次可能是两者选取了不同的湍流模型。

表1 时间步长、圆周节点数和壁面网格尺度对正方形布置四圆柱体平均阻力系数、斯特劳哈尔数的影响(L/D=4.0)

Tab. 1 Effect of time-step, nodes of circumference and wall grid scale on mean drag coefficient and Strouhal number of square-arranged four circular cylinders (L/D=4.0)

网格模型时间步长/s环向节点数首层边界厚度/mm圆柱1圆柱4CDSrCDSrCase10.0012400.318 21.0900.1940.6680.194Case20.0052400.318 21.1070.1930.6700.196Case30.0152400.318 21.0930.2020.6710.202Case40.0102400.318 21.1300.2080.6800.208Case50.0102400.100 01.1160.2050.7010.203Case60.0102401.000 01.1030.2110.6780.211Case70.0101600.318 20.8780.3570.4960.357Case80.0103200.318 21.1130.2050.6980.205文献[12]1.0040.1810.6500.181

一些学者们采用阻力系数平均值积分曲线(average drag coefficient intergral curve)判断湍流流动是否达到稳定状态[13]。基于该方法，当柱体结构的阻力系数平均值时程曲线由陡峭逐渐趋于平稳，则可认为湍流流动达到了稳定状态，即可以开始提取足够多的相关物理量进行数据处理与分析。为使涡脱落充分发展，文中模拟总时间为300 s，且数据采样的时间为t>250 s即在525个涡脱落之后，共采样了113个涡脱落周期；其中数据采集的时间间隔为0.01 s，涡脱落周期为0.47 s。

3.2 算例验证

将单圆柱体和正方形排列四圆柱体结构群(α=0°)两个经典算例的数值计算结果与已有文献结果进行对比分析，结果见图3。

图3 单圆柱尾流区时均速度和正方形排列四圆柱体的阻力系数平均值与文献结果对比Fig. 3 Comparisons between mean streamwise velocity, mean drag coefficients and the existing data of the cases of a single cylider and four square-arranged ones

从图3(a)可以看出，单圆柱体的尾流区不同特征位置处的时均速度分布与相关文献[15-16]吻合较好。从图3(b)可以看出，正方形排列四圆柱体的阻力系数平均值随间距比的变化趋势与文献[8, 12]记载的内容基本一致。从而验证了这里选取的计算模型、计算方法和相关参数设置的准确性和可靠性。为方便描述圆柱体之间的相对位置关系，特此将五个圆柱体分别命名为C1～C5。

4 计算结果与分析

在流攻角为0°、22.5°和45°工况下，主要对L/D=2.0、3.5、5.0和7.0四种间距比工况下，五圆柱体结构群绕流问题进行数值研究，分析了五圆柱体的三维涡结构、圆周平均压力及流体力系数随来流攻角和间距比的变化规律，并解释其互扰效应和流动控制机理。

4.1 流场特性分析

结构表面剪切层的发展与相互作用及尾涡结构的分布会影响流场分布特性，导致柱体表面所受的压力分布发生变化，进一步改变圆柱体表面受到的流体力。图4给出了L/D=2.0、3.5、5.0和7.0四种间距比工况下及不同来流攻角下的五圆柱体结构群瞬时涡量图。

图4 不同工况下五圆柱结构群瞬时涡量Fig. 4 Three-dimensional vorticity of five circular cylinders in different cases

当间距比较小时，由于间隙区域过于狭窄，圆柱体群的剪切层相互干扰耦合，五圆柱体结构群之间的互扰作用以临近效应为主，并且可观察到波状剪切层(corrugated shear layer，简称CSL)和尾涡第一次卷起位置FRU (the first roll up)[14]。如图4所示。在中间圆柱(C5)的阻挡作用与间隙流的高速流动共同影响下，上游两圆柱体(C1和C2)的上下表面剪切层分别斜向下和斜向上卷起CSL并形成回流区。同时，间隙流冲击C5后形成稳定的回流区，沿圆周向周围流动并冲击下游两圆柱体(C3和C4)，使得下游两圆柱体的剪切层在尾流区较远处卷起产生较大且单一的宽尾流涡街，如图4(a)和图4(b)所示。另外，C3、C4和C5处于背流面被屏蔽状态且其上下表面剪切层均未发生摆动现象。这种流动模式与文献[12]所记载的类似，称之为稳定屏蔽模式。

随间距比的增大，上游两圆柱体上下表面剪切层斜向卷起偏离角度逐渐减小，其尾流涡街发展越来越充分，撞击下游圆柱体并改变下游圆柱体的涡脱落模式，并更容易观察到CSL和FRU。下游两圆柱体受到间隙流的冲击作用及其他圆柱体尾流结构的干扰作用逐渐减弱，如图4(d)和图4(g)所示。当α=22.5°时，均匀来流斜向流经狭窄间隙区域使得C3完全处于背流面被屏蔽，从而会影响其顺流向阻力系数的大小。需要注意的是，C1和C2上下两侧流体流动特性几乎一致，而其他圆柱体流体绕流特性均存在较大差异性。该流动现象会导致五圆柱体所受的流体力发生变化。当α=45.0°时，随间距比逐步增大，两侧圆柱体(C2和C4)尾流涡街的发展受其他圆柱体的影响逐渐减小且其涡脱落模式趋于稳定。当L/D=3.5时，C1和C5的剪切层向下游延展的位移逐步增大，五圆柱体结构群之间的互扰效应以临近效应和剪切层干扰为主，如图4(f)所示。当L/D=5.0时，中游圆柱体C5尾流区会发生漩涡脱落现象。然而，C1剪切层仍会再附着在C5表面并影响其漩涡的形成和脱落，五圆柱体结构群之间的互扰效应由剪切层干扰和尾涡冲击共同作用，如图4(i)所示。

当间距比大于临界值时，上游两圆柱体尾涡结构的发展几乎不再受中间圆柱体及下游两圆柱体的干扰，其绕流特性接近于单圆柱体工况；中间与下游两圆柱体的尾流分布特性会逐渐趋于稳定。另一方面，除受到上游圆柱体尾流撞击作用外，下游两圆柱体还会受到间隙流的冲击作用，如图4(j)所示。当L/D=7.0，α=22.5°时，五圆柱体结构群之间的互扰作用以偏斜尾流冲击干扰为主。上游圆柱体的涡脱落并非直接撞击下游圆柱体，而会改变下游圆柱体两侧或一侧的流体流动特性，导致C3和C5上下表面所受的横流向流体力不一致，如图4(k)所示。另外，当α=45.0°时，C1和C5尾流区均产生周期性卡门涡街，五圆柱体结构群之间的互扰效应转变至尾流效应，如图4(l)所示。

4.2 柱体表面压力分布特性分析

图5 当来流攻角0.0°时，z/H=0.5截面处五圆柱体群时均压力系数分布特性随间距比的变化Fig. 5 Mean pressure coefficient distributions along the circumference of five cylinders versus L/D at α=0.0° and z/H=0.5

值得注意的是，L/D=2.0工况下，C3与C4的压力分布存在显著的差异性。这主要与柱体周围流场分布特性有关，如图6(a)所示。一方面是上游柱体尾流与间隙流对迎流面的冲击作用不同，另一方面是柱体结构附近尾流区漩涡形态及其影响的变化。

图6 L/D=2.0时，α=0.0°与22.5°工况下z/H=0.5截面处五圆柱体群时均流线Fig. 6 The time-averaged streamline with z/H=0.5 axial height at α=0.0° and α=22.5° with L/D=2.0

当α=22.5°时，圆柱体群之间的位置呈错列形式，内在的互扰效应更为复杂，导致柱体结构表面的压力分布特性会发生显著变化，如图7所示。间距比的变化对C1的驻点位置及其平均压力系数的影响较小，对其最小平均压力系数位置及相应数值影响较大，如图7(a)所示。特别在小间距比工况下，C2与C3的最小平均压力系数分别达到了-2.17和-1.98。同时，C3下表面没有出现正压区。这是因为柱体群间的相互作用，导致产生了间隙流高速流动现象，使得尾流区涡结构呈偏斜状，如图6(b)所示。

图7 当来流攻角22.5°时，z/H=0.5截面处五圆柱体群时均压力系数分布特性随间距比的变化Fig. 7 Mean pressure coefficient distributions along the circumference of five cylinders versus L/D at α=22.5° and z/H=0.5

图8 当来流攻角45.0°时，z/H=0.5截面处五圆柱体群时均压力系数分布特性随间距比的变化Fig. 8 Mean pressure coefficient distributions along the circumference of five cylinders versus L/D at α=45.0° and z/H=0.5

4.3 流体力系数特性分析

本节分析了α=0°、22.5°和45°工况下，五圆柱体的阻力系数和升力系数平均值及其均方根值随来流攻角和间距比的变化规律，并与单圆柱工况的结果进行了对比分析，见图9～图11所示。

当来流攻角α=0°时，该五圆柱体结构群布置形式使得上下两排并列布置双圆柱体关于y=0的中心面对称，导致上游两圆柱体(C1和C2)的流体力系数统计值几乎一致，同时，下游两圆柱体(C3和C4)也会呈现相同的趋势，见图9(a)。值得注意的是，当L/D≥5.0时，五圆柱体的尾流均已发展充分，如图4(g)和图4(j)所示，诱发其阻力系数平均值保持稳定。

图9 来流攻角0°工况，五圆柱体流体力系数随间距比的变化Fig. 9 Variation of fluid force coefficients of five circular cylinders versus L/D at α=0°

由图9(b)可知，上游两圆柱体的升力系数平均值大小几乎相同但符号相反，下游两圆柱体的升力系数平均值亦是如此，这种情况反映了上排两圆柱体与下排两圆柱体之间的相互排斥作用。特别地，当L/D=2.0时，上排两圆柱体与下排两圆柱体的升力系数平均值绝对值较大，说明此时圆柱体之间有较强的排斥作用。

不同间距比工况下，下游两圆柱体的升力系数均方根始终大于上游两圆柱体，如图9(c)所示。随间距比的增大，除L/D=2.0工况外，上游两圆柱体的升力系数均方根呈现单调递增，并逐渐接近于单圆柱体工况的结果，与文献[9]中计算结果基本一致。从图9(c)可知，当L/D≥3.5时，下游和中间圆柱体的升力系数的脉动性较强。这是因为五圆柱体的尾涡结构相互耦合与干扰从而使得C3、C4和C5三圆柱体的升力脉动性增强。

当α=22.5°时，随间距比逐渐变大，C3、C5与C2三者之间的干扰作用逐渐减弱，进而导致其阻力逐渐减小。同时，C3的阻力系数平均值单调增大至1.21。另一方面，当间距比较小时，由于C3的迎流面与背流面均受到一对漩涡冲击作用，如图6(b)所示，导致其顺流向的阻力接近零。

由图10(b)可知，因为小间距比工况下，高速间隙流的存在，使得C2和C4上下表面的流体流速有较大差别，进一步导致其上下表面存在较大压力差。随间距比逐步增大，C2和C4上下表面的流体流速差异逐步减小，其压力差也逐步减小，导致升力系数平均值变为零。特别，当L/D=7.0时，C3和C5升力系数平均值仍然不为零，这是由上游圆柱体(C1和C2)涡脱落往下游传播使得C3和C5两侧流体流动特性发生改变而引起的，如图4(k)所示。