扭曲圆柱绕流特性的大涡模拟

2021-05-08郭春雨郭航胡健

哈尔滨工程大学学报 2021年3期

郭春雨，郭航，胡健

(哈尔滨工程大学船舶工程学院，黑龙江哈尔滨 150001)

圆柱绕流一直是流体力学中的经典问题，其涉及到的边界层转捩、流动分离、再附着及旋涡脱落等物理现象都是流体力学中的研究难点。圆柱绕流流动既有不固定的分离点,又有分离后的尾流和脱体涡。随着雷诺数的增加,尾流性质,脱体涡的形态有很大的变化,具有丰富的流动现象。桥梁结构物,如拱桥的吊杆、斜拉桥的拉索,具有圆柱形状,物体表面边界层流动在逆压梯度下分离,在一定的雷诺数范围产生规则的旋涡结构,这些交替的旋涡脱落对圆柱受到的升阻力,涡激振动以及产生的噪声问题都有很大影响，因此对圆柱绕流涡脱落以及振动等问题的研究具有十分重要的意义[1]。

Schewe[2]在密封风洞中对从亚临界至超临界雷诺数范围进行了阻力、升力及斯特劳哈尔数的测量。在临界雷诺数范围内观察到了2个不连续的转变，这可以解释为2个临界雷诺数的分叉点。在这2种情况下，这些转变都伴随着临界值波动、对称性破裂(出现稳定升力)和迟滞现象。Lam等[3]利用各种实验测量技术对雷诺数为3 000～9 000扭曲圆柱的近尾迹进行了研究，结果表明扭曲圆柱的平均涡形成长度超过了一半波长且大于光滑圆柱，并且涡形成长度也直接导致阻力以及相关波动升力的减小。Lam等[4]对雷诺数为3 000的不同形状的扭曲圆柱进行了数值模拟研究，结果表明由于展长方向波浪形状的分离线导致扭曲圆柱近尾流涡结构沿展向有周期性的变化。Jung等[5]对雷诺数为3 000的3种(光滑、波浪形、扭曲形)圆柱体的绕流特性进行了数值模拟，可以得知扭曲圆柱阻力和升力与光滑圆柱相比分布减小了13%和96%，沿展长方向产生了周期性的横向涡结构，分离剪切层更加稳定且更靠近下游。Zhou等[6]在水槽以及拖曳水池中对雷诺数为7.4×103～8×104的光滑表面、沟槽状表面以及洼坑状表面3种圆柱进行了研究。结果表明沟槽状表面以及洼坑状表面圆柱的阻力及均方根升力系数比光滑圆柱要更小，而PIV结果表明圆柱表面的粗糙度导致了圆柱尾流区域以及涡脱落强度相比于光滑圆柱也要更小。Jie等[7]主要研究了触角形圆柱与光滑、椭圆以及波浪圆柱体之间的力与涡结构等方面的不同。结果表明触角形圆柱在升力方面比普通圆柱减小了79.2%，并在斯特劳哈尔数为0.2时只观察到一个相对较低幅度的波峰。

本文计算采用的雷诺数范围为亚临界雷诺数，这是由于圆柱表面流体分离后尾流是完全湍流的[8]。与上述文献中的方法不同，本研究采用的是相同直径不同扭曲角度的扭曲方式，着重于改变圆柱表面形状以达到减小圆柱受力波动以及涡激振动的目的。本研究通过三维大涡模拟方法讨论了亚临界雷诺数范围内扭曲圆柱的扭曲角度及旋转角度对其水动力性能的影响。

1 数值控制方程及模拟方法

1.1 控制方程

本文使用三维大涡模拟湍流模型来对圆柱绕流进行数值模拟[9]。大涡模拟的基本思想是直接计算大尺度涡，而通过亚格子模型来对小涡进行模拟。与时间平均不同，大涡模拟方法通过空间滤波操作来分成大涡以及小涡。因此滤波后不可压缩流体的纳维-斯托克斯以及连续方程即为大涡模拟控制方程：

(1)

(2)

通常情况下亚格子模型是基于涡粘模型的，它基于人工涡流粘度的方法，其中湍流的影响被集中到湍流粘度。该方法将亚格子尺度下的动能耗散视为类似于分子扩散。因此τij的表达形式为：

(3)

(4)

本文采用的亚格子涡粘模型为WALE(wall-adapting local-eddy viscosity)模型[10]，它是一个更现代的亚格子尺度模型，在其公式中使用了速度梯度张量的一种新形式。WALE亚格子模型提供的混合长度形式的亚格子涡粘系数为：

νt=ρΔ2Sw

(5)

式中：ρ为流体密度；Δ为长度尺度或网格过滤宽度[11]；Sw为扭曲参数。长度尺度是通过网格体积来定义的：

(6)

式中：Cw为模型系数，取值为0.544；κ为冯卡门常数，取值为0.41。扭曲参数Sw定义为：

(7)

(8)

1.2 数值方法

本文模拟采用的是三维、不可压缩、分离流、恒定密度以及非定常模拟，通过将有限体积法应用于非结构化网格上来对纳维-斯托克斯方程进行求解。在大涡模拟中，应用SIMPLE算法对压力速度耦合方程进行求解，空间离散采用的是有界中心差分离散格式[12]，时间项采用的是二阶隐式时间离散格式，梯度计算采用的是Hybrid Gauss-LSQ方法。

2 计算模型及网格和边界条件

2.1 计算模型

本模拟选择了4种不同的圆柱进行数值模拟，其中包括光滑圆柱以及3种不同扭曲角度的扭曲圆柱。3种圆柱的直径和展长相同，分别为D和6D，其中扭曲圆柱是由3段波长λ=2D的相同圆柱组成。扭曲圆柱的形状是通过扭曲角度α来定义的，而α为圆柱沿展长方向截面绕圆周上一点旋转的角度，并规定沿逆时针方向旋转为正，如图1(a)所示，且α与截面展向坐标z具有以下关系：

(9)

式中αmax为最大扭曲角度，本模拟选取的扭曲圆柱A1、A2、A3分别对应为10°、25°、40°。

由于扭曲圆柱不具有光滑圆柱那样的几何对称性，因此需要对其进行旋转不同的角度进行计算。旋转角度θ定义为圆柱绕其中心轴旋转的角度，本文选取了θ=0°，45°，90°，135°，180°，225°，270°，315°共8个角度。为了更好地对圆柱绕流特性进行分析，在展向选择了2个截面，节面和鞍面如图1(b)所示。

图1 扭曲圆柱Fig.1 The twisted cylinder

2.2 网格和边界条件

计算域网格划分采用Star CCM+生成的切割体网格，通过使用几何体曲面剪切六面体模板网格来生成体积网格。图2(a)展示的是扭曲圆柱计算域及圆柱边界层附近的体网格，为了满足大涡模拟计算要求，近壁面y+取值小于1，为了更好地捕捉圆柱尾流区的涡结构，对下游方向的网格进行了加密，且边界层处的网格尺寸以一个合适的比率向加密区逐渐增大。计算采用的无量纲时间步长定义为：T=U∞Δt/D，对于目前的所有模拟，计算了至少500个无量纲时间步长，对应大约100个旋涡脱落周期以获得更可靠的统计信息，这也要大于Lam[4]计算的时间步长以及脱落周期。

计算域和边界条件如图2(b)所示。计算域的尺寸在固定笛卡尔坐标系下在x、y、z轴方向上分别为24D×16D×6D，其中x轴方向上游边界距离原点为8D，下游边界距离原点为16D，原点距离y轴方向上下2个边界距离均为8D，展向长度则与圆柱展长相同。坐标系原点位于圆柱底部圆心上，x轴与入口来流方向相同(顺流方向)，z轴与圆柱中心轴方向平行，y轴则与x轴和z轴垂直(横流方向)。本模拟基于圆柱直径D和来流速度U∞的雷诺数Re=U∞D/ν=28 712。计算域的边界条件设置如下：入口边界设置为速度入口，出口边界设置为压力出口，由于考虑是无限长圆柱结构，因此其他4个边界均设置为对称平面，圆柱则设置为无滑移壁面。

图2 数值模拟设定Fig.2 Schematic of numerical simulation setting

3 结果与分析

3.1 网格及时间步长无关性验证

(10)

(11)

式中：ρ为流体密度；U∞为来流速度；D为圆柱截面直径；L为圆柱展向长度；FD和FL分别为总阻力和升力。斯特劳哈尔数为无量纲化的涡脱落频率(fs)，其表达式为：

(12)

式中涡脱落频率fs是通过对升力系数的时历曲线作快速傅里叶变换得到的。

首先对光滑圆柱进行网格无关性验证，表1展示了3套不同数量的网格，这里Si(i=1,2,3)分别代表粗、中和细网格对应的计算结果。以一种类似于文献[13]中处理非结构化网格的方式，网格加细比是需要在进行网格无关性验证之前确定的一个重要参数，网格加细比rG定义为：

(13)

对光滑圆柱进行时间步长无关性验证，表2展示了3个不同无量纲化时间步长，这里Si(i=4,5,6)分别代表小、中和大时间步长对应的计算结果。本文3个不同时间步长网格数量均为435万的粗网格。结果表明，只有S6与实验结果相差较大，而S4和S5则与实验数据吻合较好。因此，在保证计算精度和效率的前提下，本文最终选择了网格数量为435万的粗网格，而无量纲时间步长则为T=0.02。

表1 光滑圆柱网格无关性验证Table 1 Grid independence verification on smooth cylinder

表2 光滑圆柱时间步长无关性验证

表3 扭曲圆柱A1网格无关性验证

表4 扭曲圆柱A1时间步长无关性验证

3.2 力系数和斯特劳哈尔数

为了节省大涡模拟所需的计算资源以及时间，对扭曲圆柱不同角度之间的相似性进行验证。以α=10°的扭曲圆柱A1为例，选取了8个旋转角度(θ=0°，45°，90°，135°，180°，225°，270°，315°)进行计算，计算结果如表5所示，从计算结果可以看出，扭曲圆柱间隔角度180°的力系数以及斯特劳哈尔数近似相同，可以初步判断扭曲圆柱不同角度之间具有部分的相似性。

表5 扭曲圆柱A1不同旋转角度的力系数和St

图3 3种扭曲圆柱与光滑圆柱的对比Fig.3 Comparison of three twisted cylinders and the smooth cylinder

3.3 尾涡形成长度

圆柱尾涡形成长度是一个很重要的量，因为它影响圆柱尾部压力以及力系数。对于尾涡形成长度的定义不是普适的，不同学者有不同的方法[15-19]。通常用尾流中心线(y=0平面)上无量纲化的平均流向速度为零的时间平均闭合点位置(PU0)以及无量纲化的流向速度波动均方根的最大值点位置(PUrms)来定义涡形成长度。因此尾涡形成长度包括平均流向速度为零(U/U∞=0)对应的尾流闭合长度Lfc和流向速度波动均方根(u′/U∞)最大值对应的最大湍流强度长度Lfu。

图4为3种扭曲圆柱沿展长方向一个波长长度的尾流闭合长度Lfc和最大湍流强度长度Lfu与光滑圆柱的对比。结果表明，圆柱A1的Lfc和Lfu与光滑圆柱相比几乎没有增长。当θ=0°时，扭曲圆柱A2的Lfc和Lfu与光滑圆柱相比最大分别增大了14%和24%，而在其他旋转角度Lfc增大幅值达到23%～96%，Lfu增大幅值达到31%～105%。当θ=0°时，扭曲圆柱A3的Lfc和Lfu与光滑圆柱相比最大分别增大了22%和27%，相比光滑圆柱有较大增幅，在其他旋转角度Lfc最大增幅为98%，Lfu增大幅值达到39%～114%。

综合图4分析可知，Lfc和Lfu在θ=0°时波动幅值较小，而在其他旋转角度时沿展长方向基本呈现增大的趋势，并且随着圆柱扭曲角度α的增大，Lfc和Lfu的波动幅值也逐渐增大，特别的扭曲圆柱A3甚至出现了Lfc的最小值小于光滑圆柱的情况，原因可能为随着圆柱表面曲率逐渐增大而导致尾流闭合长度的波动范围增大。当θ=0°时，Lfc和Lfu随着α的增大也逐渐增大；当θ为其他值时，Lfc和Lfu随着α的增幅变化不是很大，并且均大于θ=0°时的Lfc和Lfu的增幅。通过对比可以得出结论扭曲圆柱升阻力随着涡形成长度Lfc和Lfu的增大而减小，当θ=0°时,涡尾流长度相对于光滑圆柱增幅比较小，因此对圆柱尾流涡结构的控制以及抑制幅度也比较小，使得升阻力下降也比较小；相反其他角度尾流长度得到很大的延伸表明能够对尾涡结构很好地控制，从而能够使升阻力波动最小化。

3.4 压力系数分布

图4 扭曲圆柱与光滑圆柱的尾流闭合长度Lfc和最大湍流强度长度Lfu对比Fig.4 Comparison of the wake closure length Lfc and maximum turbulence intensity length Lfu of twisted cylinders and the smooth cylinder

图5 扭曲圆柱与光滑圆柱的平均圆周压力系数对比Fig.5 Comparison of the average circumferential pressure coefficient of twisted cylinders and the smooth cylinder

3.5 湍流动能

图6 在x-y平面的无量纲化湍流动能TKE分布对比Fig.6 Comparison of the normalized turbulent kinetic energy (TKE) distributions in the x-y plane

从图6(a)可以看出，光滑圆柱最大TKE区间位于圆柱近尾流区，且该区间内TKE值很大。对比图6(b)～(e)可以发现，当θ=0°时，扭曲圆柱A2在节面和鞍面与光滑圆柱相比TKE的分布基本相同。而当θ=90°时，圆柱A2的最大TKE区间相比θ=0°时有较大的后移，同时区间内的TKE值有大幅度的减小，并且鞍面的TKE值要比节面更小。根据图6(f)～(i)可以看出，当θ=0°时，圆柱A3的最大TKE区间相比A2也有后移，区间内的TKE值也有大幅度的减小。当θ=90°时，圆柱A3的最大TKE区间位置与圆柱A2基本相同，区间内的TKE值略有减小，并且鞍面的TKE值仍然比节面要小。

综合图6分析可知，随着α的增大，圆柱的最大值TKE区间位置逐渐后移，远离圆柱表面，沿流向的TKE分布变化逐渐变小，且最大及整体TKE值相比光滑圆柱也在逐渐减小，尤其在α=40°时有很大程度的减幅。而随着α的增大，鞍面的TKE值相比于节面也越来越小。当θ=90°时整体的TKE值均小于θ=0°，这也与3.2节中升阻力系数大小以及3.3节中尾涡形成长度曲线相一致，因此可以得出结论，圆柱尾流TKE的显著减小能够导致圆柱的阻力以及波动升力的减小。

3.6 三维涡结构

图7展示了不同圆柱的三维瞬时涡结构等值面，对比了自由剪切层发展和涡形成长度的显著差异。

图7 三维瞬时涡结构等值面对比Fig.7 Comparison of the three-dimensional instantaneous vortex structure isosurfaces

图7(a)为光滑圆柱的三维涡结构云图，可以观察到在圆柱近尾流区域存在大尺度的卡门涡街，并且振荡幅度较大。通过观察图7(b)～(g)能够发现，圆柱A1尾流区域涡振荡幅值与光滑圆柱相比略有减小。当θ=0°时，圆柱A2尾流区域涡分布情况与光滑圆柱基本相同，而随着α的增大，圆柱尾涡的脱落位置延迟得更加靠近下游，同时还能观察到卡门涡街在远尾流区域重新出现，但是其波动幅值较光滑圆柱减小了。当θ=90°时，其他2种圆柱剪切层沿下游方向延伸得更长了，并且涡的脱落位置也更加远离圆柱表面。同时由于初始涡结构的形成沿下游方向更加延迟了，因此卡门涡街强度沿下游方向也更加弱化了。

综合对图7进行分析可知，与光滑圆柱相比，随着α的增大，3种圆柱尾涡脱落位置以及剪切层长度沿顺流方向的延伸逐渐增加，卡门涡街在更远离圆柱表面的下游位置重新出现，但是其波动幅值也逐渐减小。并且可以看到，θ=90°时涡的抑制效果均比θ=0°时要更好，延伸长度也更大，这也与3.2节中力系数以及3.5节中TKE的分布情况相吻合。由于圆柱剪切层的延伸伴随着较弱的尾涡波动，从而直接导致了阻力以及升力振荡幅值的减小。