圆柱绕流一直是流体力学中的经典问题，其涉及到的边界层转捩、流动分离、再附着及旋涡脱落等物理现象都是流体力学中的研究难点。圆柱绕流流动既有不固定的分离点, 又有分离后的尾流和脱体涡。随着雷诺数的增加, 尾流性质, 脱体涡的形态有很大的变化, 具有丰富的流动现象。桥梁结构物, 如拱桥的吊杆、斜拉桥的拉索, 具有圆柱形状, 物体表面边界层流动在逆压梯度下分离, 在一定的雷诺数范围产生规则的旋涡结构, 这些交替的旋涡脱落对圆柱受到的升阻力, 涡激振动以及产生的噪声问题都有很大影响，因此对圆柱绕流涡脱落以及振动等问题的研究具有十分重要的意义^[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×10³~8×10⁴的光滑表面、沟槽状表面以及洼坑状表面3种圆柱进行了研究。结果表明沟槽状表面以及洼坑状表面圆柱的阻力及均方根升力系数比光滑圆柱要更小，而PIV结果表明圆柱表面的粗糙度导致了圆柱尾流区域以及涡脱落强度相比于光滑圆柱也要更小。Jie等^[7]主要研究了触角形圆柱与光滑、椭圆以及波浪圆柱体之间的力与涡结构等方面的不同。结果表明触角形圆柱在升力方面比普通圆柱减小了79.2%，并在斯特劳哈尔数为0.2时只观察到一个相对较低幅度的波峰。

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

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

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

$ {\frac{{\partial {{\bar u}_i}}}{{\partial t}} + \frac{{\partial {{\bar u}_i}{{\bar u}_j}}}{{\partial {x_j}}} = - \frac{1}{\rho }\frac{{\partial \bar p}}{{\partial {x_i}}} + \nu \frac{{{\partial ^2}{{\bar u}_i}}}{{\partial {x_j}\partial {x_j}}} - \frac{{\partial {\tau _{ij}}}}{{\partial {x_j}}}} $

(1)

$ {\frac{{\partial {{\bar u}_i}}}{{\partial {x_i}}} = 0} $

(2)

式中：u_i是沿笛卡尔坐标系x_i方向的滤波速度分量; t为时间; p为压力; ρ为流体密度; ν为流体运动粘度。亚格子应力张量τ_ij定义为 $ {{\boldsymbol{\tau}} _{ij}} = \overline {{u_i}{u_j}} - {\bar u_i}{\bar u_j}$ ，和雷诺应力相仿，亚格子应力是过滤掉小尺度脉动和可解尺度湍流间的动量输运。要实现大涡数值模拟，必须构造亚格子应力的封闭模式。亚格子应力湍流模型是实现大涡数值模拟的关键。

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

$ {\mathit{\boldsymbol{\tau }}_{ij}} - \frac{1}{3}{\mathit{\boldsymbol{\tau }}_{kk}}{\delta _{ij}} = - 2{\nu_t}{\mathit{\boldsymbol{\bar S}}_{ij}} $

(3)

式中：ν_t为湍流涡粘系数；S_ij为可解尺度的扭曲率张量，定义为：

$ {{\mathit{\boldsymbol{\bar S}}}_{ij}} = \frac{1}{2}\left( {\frac{{\partial {{\bar u}_i}}}{{\partial {x_j}}} + \frac{{\partial {{\bar u}_j}}}{{\partial {x_i}}}} \right) $

(4)

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

$ {\nu _t} = \rho {\mathit{\Delta }^2}{\mathit{\boldsymbol{S}}_w} $

(5)

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

$ \varDelta=\left\{\begin{array}{ll} C_{w} V^{1 / 3}, & \text { 不应用长度尺度限制 } \\ \min \left(\kappa d, C_{w} V^{1 / 3}\right), & \text { 应用长度尺度限制 } \end{array}\right. $

(6)

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

$ \boldsymbol{S}_{w}=\frac{\left(\boldsymbol{S}_{i j}^{d} \boldsymbol{S}_{i j}^{d}\right)^{3 / 2}}{\left(\overline{\boldsymbol{S}}_{i j} \overline{\boldsymbol{S}}_{i j}\right)^{5 / 2}+\left(\boldsymbol{S}_{i j}^{d} \boldsymbol{S}_{i j}^{d}\right)^{5 / 4}} $

(7)

其中张量S_ij^d定义为：

$ \boldsymbol{S}_{i j}^{d}=\frac{1}{2}\left(\bar{g}_{i j}^{2}+\bar{g}_{j i}^{2}\right)-\frac{1}{3} \delta_{i j} \bar{g}_{k k}^{2}, \bar{g}_{i j}=\frac{\partial \bar{u}_{i}}{\partial x_{j}} $

(8)

1.2 数值方法

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

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

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

$ \alpha = {\alpha _{{\rm{max}}}}\sin \left( {\frac{{{\rm{ \mathsf{ π} }}z}}{D}} \right) $

(9)

式中α_max为最大扭曲角度，本模拟选取的扭曲圆柱A₁、A₂、A₃分别对应为10°、25°、40°。

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

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个边界均设置为对称平面，圆柱则设置为无滑移壁面。

3 结果与分析 3.1 网格及时间步长无关性验证

网格数量及时间步长的设置对大涡模拟方法的计算精度有很大的影响，因此需要对光滑圆柱和扭曲圆柱分别进行网格以及时间步长无关性验证。为了定量验证数值计算方法的准确性，将光滑圆柱时间平均阻力系数(C_D)，升力系数的波动均方根值(C_{L, rms})以及斯特劳哈尔数(St)的计算结果与实验值进行对比。其中阻力系数(C_D)和升力系数(C_L)的定义为：

$ {{C_D} = \frac{{2{F_D}}}{{\rho U_\infty ^2DL}}} $

(10)

$ {{C_L} = \frac{{2{F_L}}}{{\rho U_\infty ^2DL}}} $

(11)

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

$ St = \frac{{{f_s}D}}{{{U_\infty }}} $

(12)

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

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

$ {r_G} = {\left( {\frac{{{N_{{\rm{fine}}}}}}{{{N_{{\rm{coarse}}}}}}} \right)^{1/d}} $

(13)

式中：N为总的网格数; d为计算问题的维度，本文为三维问题，所以d取3。本文3套网格的加细比约为1.2，无量纲化时间步长均取为T=0.02。表 1给出了雷诺数Re=28 712时光滑圆柱网格收敛性验证的计算结果以及与实验值^[14]的对比。从表中可以看出，通过LES方法计算得到的C_D、C_{L, rms}与St对3套网格的密度敏感度不高，同时与实验值吻合得也比较好。

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

对于扭曲圆柱同样进行了网格以及时间步长无关性验证, 选择α=10°的扭曲圆柱A₁(θ=0°)进行计算，结果分别如表 3和表 4所示。为了更好地与光滑圆柱进行对比，扭曲圆柱采用的网格数量以及时间步长与光滑圆柱均相同。从表中可以看出网格数量以及时间步长对C_D、C_{L, rms}与St计算结果的影响均比较小，因此可以得出扭曲圆柱的网格以及时间步长无关性验证结果是很好的。最终选择网格数量为435万的粗网格，无量纲时间步长取T=0.02。

3.2 力系数和斯特劳哈尔数

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

由于扭曲圆柱的部分相似性，因此只选取了4个旋转角度(θ=0°，45°，90°，135°)进行计算，图 3展示了Re=28 712时光滑圆柱和3种扭曲圆柱的力系数以及斯特劳哈尔数。从图 3(a)可以看出，3种扭曲圆柱的斯特劳哈尔数均在0.2附近，与光滑圆柱相比结果相差不大，Sumer等^[14]中光滑圆柱在亚临界雷诺数范围内的斯特劳哈尔数均在0.2附近，因此可以得出结论雷诺数或圆柱表面几何形状的变化对其尾涡脱落频率的影响很小。从图 3(b)可以看出，圆柱A₁与光滑圆柱相比阻力系数相差不大。而除了θ=0°，其他3个旋转角度对应的2种圆柱的平均阻力系数C_D相比光滑圆柱都减少了14%~18%左右；当θ=0°时，扭曲圆柱A₂与光滑圆柱相差不大，而A₃减小达到了11%。图 3(c)与图 3(b)相似，圆柱A₁与光滑圆柱相比升力系数相比略有减小。除了θ=0°，其他3个旋转角度对应的2种圆柱的升力波动均方根值C_{L, rms}相比光滑圆柱都减少了60%~85%左右。当θ=0°时，扭曲圆柱A₂减小了约9%，而A₃减小达到了68%。根据图 3(b)和图 3(c)可知，仅扭曲圆柱A₃相比光滑圆柱阻力以及升力系数减小幅度比较大，其他2种扭曲圆柱在θ=0°时阻力以及升力系数减小幅度很小甚至增大。

3.3 尾涡形成长度

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

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

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

3.4 压力系数分布

对于不同圆柱的阻力减额可以通过分析圆周压力分布来进行预报^[20]。图 5为不同扭曲圆柱在节面和鞍面处的平均圆周压力系数C_p的分布图，从图中可以看到，不同截面处最小压力系数位置以及大小均不同。同时将光滑圆柱也包含进去以进行对比，光滑圆柱最小压力系数点位置为70°左右。对于扭曲圆柱A₁，其最小压力系数位置以及大小与光滑圆柱相差不大。而对于圆柱A₂，在θ=0°时2个截面的最小压力系数位置和大小与光滑圆柱相差不大，随着θ逐渐增大，节面位置随之增大，而鞍面位置则逐渐减小，最小压力系数大小也逐渐增大，且与圆柱A₁相比均有所增加。圆柱A₃的4个旋转角度在节面和鞍面处的最小压力系数位置及大小与光滑圆柱相比均有较大的增幅，其他规律则与圆柱A₂基本相同。

综合图 5可以看出，除圆柱A₁与光滑圆柱基本相同外，当θ=0°时其他2种扭曲圆柱在节面和鞍面处的最小压力系数位置和大小相比光滑圆柱均随着α的增大而逐渐加，而其他3个旋转角度变化规律则基本相同。而在亚临界雷诺数范围内，平均阻力系数与-C_p是成正比的^[5]，表明压力系数越大圆柱的平均阻力越小，这也与3.2、3.3节中圆柱的力系数以及尾涡形成长度的变化规律相吻合。

3.5 湍流动能

物体周围的湍流波动与施加在物体上的波动力有关^[21]。图 6展示了x-y平面上不同圆柱尾涡的无量纲化湍流动能云图(TKE=(u′²+v′²+w′²)/2U_∞²)。

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

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

3.6 三维涡结构

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

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

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

4 结论

1) 扭曲圆柱A₃相比光滑圆柱阻力以及升力系数减小幅度比较大，斯特劳哈尔数基本没有变化。

2) 尾流长度得到很大的延伸，同时最小压力系数位置和大小也增大了。

3) 圆柱尾流湍动能显著减小，三维涡结构的波动也得到了抑制。

因此扭曲圆柱A₃对尾涡结构有很好的控制以及使自由剪切层更加稳定的发展，并且也为圆柱提供了更大的阻力减额以及更好的抑制升力波动，能够显著减小圆柱的涡激振动现象。