一种基于虚拟格栅的LES湍流生成方法及圆柱绕流研究

An LES Turbulence Generation Method Based on Virtual Grid and Its Application of Flow around a Circular Cylinder

摘要

为提高大涡数值模拟（large eddy simulation, LES）中生成格栅湍流风场效率，对一种LES湍流生成方法进行改进，在入口边界处设置虚拟格栅网获取网状射流，并将其成功应用于Re=3 900的圆柱绕流模拟.结果表明：基于虚拟格栅生成的湍流风场具有各向同性特征，各湍流强度沿流向呈指数衰减趋势，湍流积分尺度沿流向呈线性增长趋势，且脉动风速功率谱与von Kármán谱吻合较好.此外，相较于均匀来流，湍流来流下圆柱的平均阻力系数减小，而阻力系数和升力系数的标准差均大幅增加，即湍流会加剧圆柱升阻力的脉动效应.来流湍流对圆柱表面的流动分离有抑制作用，可使分离点位置后移，导致回流长度和斯脱罗哈数减小，同时使尾流涡旋脱落不稳定、涡脱能量更强.

Abstract

To further improve the efficiency of grid-generated turbulent wind field for large eddy simulation （LES）， an improved turbulence generation method based on LES is proposed in this paper， in which grid-jet can be obtained on the inlet boundary. Then， LES of the flow around a circular cylinder at Re=3 900 is successfully carried out. The results indicate that the turbulent wind field generated based on virtual grid exhibits isotropic characteristics. The turbulence intensity decreases exponentially along the flow direction， and the turbulence integral scale increases linearly along the flow direction. The fluctuating wind power spectrum is consistent with the von Kármán spectrum. In addition， compared to a smooth inflow， the mean drag coefficient of the cylinder decreases in the turbulent inflow， whereas the standard deviation of the drag coefficients and lift coefficients increase significantly， indicating stronger pulsations of the lift and drag forces in turbulent flow. Turbulence also delays flow separation on the cylinder surface and inhibits separation， resulting in a shorter recirculation length and a lower Strouhal number. Vortex shedding becomes unstable and the shedding energy is stronger.

Graphical abstract

关键词

虚拟格栅 / 湍流 / 圆柱绕流 / 大涡模拟 / 流场

湍流是自然界中流体常见的流动状态，其流场中压力和速度等物理量都伴随着强烈的变化，土木工程结构大多位于湍动剧烈的大气对流层内. 在进行结构抗风研究时，通常将自然风分解为不随时间变化的平均分量和随时间变化的脉动分量两部分. 早期研究大多关注平均风特性对工程结构的影响，却忽略了来流风的脉动分量的影响. 但实际工程中，脉动风会影响结构的疲劳寿命^［1］和气动稳定^［2］等性能，甚至引发风致振动等严重后果. Vickery^［3］在湍流下方柱气动力研究中发现，湍流对气动力影响显著，当湍流积分尺度与模型长度相近时，方柱升力的相关性增大. Kawatani等^［4］针对多种基础断面在湍流下的涡激振动研究表明，湍流对涡振振幅的影响表现为有时放大，有时抑制. 因此，在结构抗风设计中脉动风的影响不可忽略.

湍流下结构抗风研究中，均匀各向同性是湍流研究中的重要理想假设，如何高效且准确地生成均匀各向同性湍流，仍是当前研究亟待解决的问题. 风洞试验中生成湍流的方法主要分为被动模拟^［5-8］和主动模拟^［9-11］两种. 被动模拟通过摆放格栅、尖劈、粗糙元等阻碍装置生成实验需要的湍流风场，具有低成本、易操作等特点，是风洞试验中生成湍流风场的常用方法之一. 受被动模拟方法自身限制，其流场中的湍流积分尺度通常难以满足实验缩尺比的需要，在模拟不同流场时需要重新调整摆放. 主动模拟利用震动栅板、旋转翼等装置，通过向流场输入适当频率的能量来调节流场指标，这虽然弥补了被动模拟方法中积分尺度不满足缩尺比要求的缺点，但也使得经济成本和技术要求较高. 随着计算流体动力学（computational fluid dynamics，CFD）技术和计算机硬件的飞速发展，基于CFD的数值模拟方法已在桥梁及高层建筑抗风中得到广泛的应用，并迅速成为空气动力研究的有力工具. 与传统风洞试验相比，数值模拟具有成本低、速度快和能提供可视的信息等诸多优点. LES是目前数值模拟中一种广泛使用的方法，周桐等^［12］认为目前在运用LES产生湍流的方法中，主要有预前模拟法^［13-16］和人工合成法^［17-21］两大类. 预前模拟法主要分为被动模拟法和“回收-变换”法两种；通过在计算域内建立格栅、尖劈、粗糙元等实体模型，一定程度上再现了风洞试验中被动模拟方式产生的风场，具有技术成本低、流场结构相对真实的特点. 但在模拟不同的目标风场时，需要重新建立实体模型并重新划分网格，严重降低了工作效率. 同时，引入实体模型意味着需要进行额外的建模和网格划分工作，导致时间和计算成本大幅提高. 人工合成法通过在计算域入口处设置数学函数来获得目标湍流流场，具有灵活、高效率等优点，但在计算域中湍流的自保持性较差.

为弥补预前模拟法生成格栅湍流时计算成本高的缺陷，本文在格栅湍流生成方式上进行了创新，摒弃了实体格栅建模思路，提出虚拟格栅生成格栅湍流的方法. 该方法通过调整数值输入即可获得不同的湍流风场，省去传统方法需在计算域入口建立实体模型的繁琐操作，且无需额外网格划分，大幅提升生成效率并节省计算资源，同时兼顾了预前模拟法技术成本低、流场结构相对真实的特点. 基于此湍流生成方法，进行了湍流风场下Re=3 900的圆柱绕流研究，得到了湍流下圆柱绕流特性，并对结果进行分析.

1 湍流生成方法

通过设置虚拟格栅函数表达式的方式，调整计算域中入口面上的速度分布，从而得到目标虚拟格栅速度入口. 虚拟格栅的函数表达式为

U=\left\{\begin{array}{l}{U}_{\mathrm{0}},i\in \left[\frac{a_i}{\mathrm{2}}+\left(n_i-\mathrm{1}\right)M_i,n_i{M}_i-\frac{a_i}{\mathrm{2}}\right]\\ \mathrm{0},i\notin \left[\frac{a_i}{\mathrm{2}}+\left(n_i-\mathrm{1}\right)M_i,n_i{M}_i-\frac{a_i}{\mathrm{2}}\right]\end{array}\right.（1）

式中：i=x，y；U为入口面上的速度；U_{\mathrm{0}}为虚拟格栅孔口处速度；a_i为i方向上的格栅条宽度；M_i为i方向上单个格栅宽度；n_i为沿i方向上布置的格栅数量，如图1所示. 与实体格栅通过阻碍和扰动流场生成湍流的原理不同，虚拟格栅通过在计算域入口模拟格栅射流的方式形成湍流. 本文虚拟格栅沿y方向均匀布置4个，z方向均匀布置3个. 格栅宽度均设为1 m，格栅条宽度均为0.4 m，虚拟格栅孔口处的速度为 U₀=1 m/s.

2 数值方法

2.1 湍流模型和求解器设置

采用LES模型来描述三维不可压流体流动，利用滤波函数分离大尺度的涡和小尺度的涡，大尺度的涡通过N-S方程直接求解，小尺度的涡用亚格子模型（sub-grid scale， SGS）来封闭. 滤波后的瞬时不可压N-S方程如下：

\frac{\partial {\tilde{u}}_i}{\partial x_i}=\mathrm{0}（2）

\frac{\partial {\tilde{u}}_i}{\partial t}+\frac{\partial {\tilde{u}}_i{\tilde{u}}_j}{\partial x_j}=-\frac{\mathrm{1}}{\rho }\frac{\partial \tilde{p}}{\partial x_i}+v\frac{{\partial }^{\mathrm{2}}{\tilde{u}}_i}{\partial x_j\partial x_j}-\frac{\partial {\tau }_{ij}}{\partial x_j}（3）

式中：i,j=\mathrm{1,2},\mathrm{3}；x_i分别表示为笛卡尔坐标系中三个轴坐标；{\tilde{u}}_i表示为笛卡尔坐标系中滤波后三个方向上的速度分量；t为时间；v和\rho分别为流体的运动黏性系数和密度；\tilde{p}为滤波后的压力； {\tau }_{ij}^{}为亚格子尺度雷诺应力（SGS Renolds stress），基于Boussinesq假设，可表示为

{\tau }_{ij}-\frac{\mathrm{1}}{\mathrm{3}}{\tau }_{kk}{\delta }_{ij}=-\mathrm{2}{v}_{\mathrm{s}\mathrm{g}\mathrm{s}}{\tilde{S}}_{ij}（4）

式中：{\delta }_{ij}为克罗内克符号；v_{\mathrm{s}\mathrm{g}\mathrm{s}}为SGS涡黏系数；{\tilde{S}}_{ij}为滤波后的应变率张量，表示为

{\tilde{S}}_{ij}=\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\partial {\tilde{u}}_i}{\partial x_j}+\frac{\partial {\tilde{u}}_j}{\partial x_i}\right)（5）

采用局部涡黏度壁面自适应模型（wall-adapting local eddy viscosity model，WALE），WALE模型的涡黏系数为

v_{\mathrm{s}\mathrm{g}\mathrm{s}}=L_{\mathrm{s}}^{\mathrm{2}}\frac{{\left(S_{ij}^{\mathrm{d}}{S}_{ij}^{\mathrm{d}}\right)}^{\mathrm{3}/\mathrm{2}}}{{\left({\tilde{S}}_{ij}{\tilde{S}}_{ij}\right)}^{\mathrm{5}/\mathrm{2}}+{\left(S_{ij}^{\mathrm{d}}{S}_{ij}^{\mathrm{d}}\right)}^{\mathrm{5}/\mathrm{4}}}（6）

式中：亚格子尺度的混合长度L_{\mathrm{s}}和应变率张量S_{ij}^{\mathrm{d}}的表达式分别为

L_{\mathrm{s}}=\mathrm{m}\mathrm{i}\mathrm{n}\left(kd,C_{\mathrm{w}}{\Delta }^{\mathrm{1}/\mathrm{3}}\right)（7）

S_{ij}^{\mathrm{d}}=\frac{\mathrm{1}}{\mathrm{2}}\left({\tilde{g}}_{ij}^{\mathrm{2}}+{\tilde{g}}_{ji}^{\mathrm{2}}\right)-\frac{\mathrm{1}}{\mathrm{3}}{\delta }_{ij}{\tilde{g}}_{kk}^{\mathrm{2}}（8）

式中：C_{\mathrm{w}}=\mathrm{0.325}为WALE模型系数；d为到最近壁面的距离；k为冯卡门常数；\Delta =\Delta x\Delta y\Delta z为滤波尺度.其中：

{\tilde{g}}_{ij}=\frac{\partial {\tilde{u}}_i}{\partial x_j}（9）

数值模拟基于开源CFD软件OpenFOAM，采用自适应时间步长，选用pimpleFoam求解器（PIMPLE算法，一种混合迭代方法，通过PISO和SIMPLE算法的结合，实现压力和速度场的高效耦合计算）进行求解，残差收敛标准为1.0×{\mathrm{10}}^{-\mathrm{6}}，数值离散各项参数格式设置如表1所示.

2.2 计算域和边界条件

图2为计算域示意图，为了使计算域中的流体充分发展，沿x方向长度取20 m，计算域y方向长度为4 m，z方向长度为3 m. 计算域入口采用速度入口，入口处的速度函数分布如第1节所示，出口面为压力出口边界，四周壁面均设置为滑移边界条件^［22-24］，这种边界条件可使流体在边界上的切向速度不变，法向速度为零，此时的边界为绝对光滑的理想壁面. 采用这种设置是为了方便计算，简化模型，同时最大限度地减少湍流耗散.

2.3 网格无关性验证

网格大小显著影响模拟结果的准确性，网格尺寸过大导致计算结果不准确，过小则浪费计算资源，因此网格无关性验证极其重要. 本文中计算域内无实体模型，为减小网格对湍流流场发展的影响，直接采用六面体网格单元对计算域进行均匀划分，如图3所示. 测试网格有数量为192万、375万、889万三种网格，对应网格单元边长分别为0.05 m、0.04 m和0.03 m. 在计算域流动方向的几何中心线上，每隔0.5 m设置一个监测点，采集1 000 s的速度数据. 为确保湍流已充分发展，选择后800 s的数据作为统计样本.

图4和图5为不同网格分辨率下流向湍流强度和湍流积分尺度的发展曲线（图中x/M为格栅宽度对横坐标的归一化，M是格栅宽度）.由图4、图5可以看出，网格分辨率影响流场参数大小，但参数曲线的衰减幂律不受网格分辨率的影响. 随着网格分辨率增加，湍流强度发展曲线逐渐向下平移，375万和889万网格的湍流强度在x/M=7.5后的区域偏差逐渐减小. 192万网格的湍流积分尺度发展曲线波动较大，随着网格分辨率增加，湍流积分尺度发展曲线波动逐渐减弱，375万和889万网格下的湍流积分尺度波动较小，曲线更平顺.

图6为9 m处不同网格分辨率下流向脉动风速功率谱.由图可以看出，375万网格的脉动风速功率谱形状与889万网格大致相同，而192万网格与889万网格模拟结果差异较大. 随着网格分辨率增加，低频部分功率略微减小，高频部分明显增大，这是因为网格分辨率增大使得可求解的网格尺度减小，高频部分的解析会更加精确，网格尺寸的减小也降低了数值耗散带来的能量损耗. 因此，选用375万的网格数量可满足研究要求.

3 流场检验

为检验该方法产生的湍流能否应用于数值模拟研究，需要对流场中湍流强度、湍流积分尺度和脉动风速功率谱等参数进行验证. 图7展示了流向中心线上平均风速的衰减情况.由图7可知，三个方向上的平均速度在x/M=7位置后基本不变. 湍流的特征是不规则性，流场中各点处速度均具有脉动特性，这些脉动源于流场中涡旋的无规则运动.图8为计算域几何中心点处风速时程曲线. 由图8可以看出，三个方向的速度分量均有明显的脉动特性.

图9展示了流场的三维瞬时Q准则等值面图. 由图9可知，流场中涡旋具有不规则性、随机性和多尺度性，随着入口距离增加，涡旋的尺度不断增大，并且沿流向的尺度大于其他两个方向的尺度，这与湍流积分尺度的发展趋势一致. 为进一步了解流场内部的湍流发展情况，截取z=1.5 m处x-y平面的涡量图，如图10所示. 随着入口距离增大，涡量逐渐减弱，能量主要集中在靠近入口侧，这与湍流强度的发展趋势一致. 图11展示了计算域内不同位置处y-z平面的涡量图.由图11可知，当流体流入后，流入流体与域中流体发生掺混，由于两者存在速度差，在黏性剪切力作用下，域中静止的流体随流入流体开始运动，涡旋也随之产生. 起初，该区域的流动表现为层流状态. 随入口距离增加，涡旋的尺寸增大直到完全阻断射流，并充斥整个截面，这个过程流动从层流逐渐变为湍流，伴随着强烈的能量传递，且能量较为集中，此时流动呈非均匀各向同性，如图11（a）~（c）所示. 随后涡旋开始坍塌，能量迅速耗散，这是导致湍流强度迅速衰减的重要原因，如图11（c）~（e）所示. 随涡旋的发展，能量传递开始减弱，流动逐步稳定且呈均匀各向同性，湍流相关的气动力实验通常在这个区域进行，如图11（f）~（h）所示.

3.1 湍流强度及积分尺度

湍流强度是衡量流体湍动程度的量度，反映了湍流流场中某一方向上风的平均脉动分量占平均分量的比例. 湍流强度定义为脉动风速的标准差与平均风速的比值，即

I_i={\sigma }_i/U（10）

式中：i=\mathrm{u},\mathrm{v},\mathrm{w}；{\sigma }_i为脉动风速i的标准差；U为测点平均风速. 湍流强度越高表明流体中的扰动和波动越强，对于均匀各向同性湍流，存在I_{\mathrm{u}}=I_{\mathrm{v}}=I_{\mathrm{w}}的关系. 图12给出了三个速度分量的湍流强度发展曲线，从图中可以看出，湍流强度随着入口距离的增加呈指数衰减，x/M≥9后衰减趋势渐缓，且趋于均匀各向同性. Mohamed等^［25］提出湍流强度衰减曲线的幂律关系如下：

I^{\mathrm{2}}=A{\left[\left(x-x_{\mathrm{0}}\right)/M\right]}^{-n}（11）

式中：n为衰减指数；x为距离格栅下游的距离；x_{\mathrm{0}}为原点；A为衰减系数. 对顺风向湍流强度进行拟合，原点取零，得到衰减指数n=1.12，如图13所示. Yu等^［15］获得的衰减指数为1.09<n<1.26.

湍流积分尺度是衡量湍流场中某一方向上涡旋平均尺度的指标，定义为湍流相关系数的无穷积分，可以用时间或者空间来量度，即

L_i^r=\frac{\mathrm{1}}{{\sigma }_i^{\mathrm{2}}}{\int }_{\mathrm{0}}^{\mathrm{\infty }}{C}_{i_{\mathrm{1}}{i}_{\mathrm{2}}}\left(r\right)\mathrm{d}r（12）

式中：r=x,y,z；{\sigma }_i^{\mathrm{2}}为脉动风速i的方差；C_{i_{\mathrm{1}}{i}_{\mathrm{2}}}\left(r\right)为相距r的两点脉动风速分量i_{\mathrm{1}}=i\left(x_{\mathrm{1}},y_{\mathrm{1}},z_{\mathrm{1}},t\right)和i_{\mathrm{2}}=i\left(x_{\mathrm{1}}+x,y_{\mathrm{1}},z_{\mathrm{1}},t\right)之间的互协方差函数.C_{i_{\mathrm{1}}{i}_{\mathrm{2}}}\left(\mathrm{0}\right)={\sigma }_i^{\mathrm{2}}. 湍流积分尺度越大表明起主导作用的涡旋越大，对于均匀各向同性湍流，则有L_{\mathrm{u}}^x=\mathrm{2}{L}_{\mathrm{v}}^x=\mathrm{2}{L}_{\mathrm{w}}^x. 根据Taylor冻结假设^［26］，仅需要一个点的速度信息，即可得到该点的湍流积分尺度，即

L_i^r=\frac{U}{{\sigma }_i^{\mathrm{2}}}{\int }_{\mathrm{0}}^{\mathrm{\infty }}{R}_i\left(\tau \right)\mathrm{d}\tau（13）

式中：R_i\left(\tau \right)是脉动风速分量i\left(x_{\mathrm{1}},t+\tau \right)的自相关函数，R_i\left(\mathrm{0}\right)={\sigma }_i^{\mathrm{2}}. Flay等^［27］建议积分上限取R_i\left(\tau \right)=\mathrm{0.05}{\sigma }_i^{\mathrm{2}}. 图14给出了湍流积分尺度沿流向发展曲线，可以看出，随着距离入口距离增大，湍流积分尺度近似线性增大，且具有较好的均匀各向同性. 湍流积分尺度的幂律关系为

L=B{\left(\frac{x-x_{\mathrm{0}}}{M}\right)}^k（14）

式中：B为常数；k为积分尺度的增长指数. 本文得到的流向湍流积分尺度增长指数k≈0.7，如图15所示. 这与邱旭等^［16］得到的湍流积分尺度随入口距离增加线性增大结论相吻合.

以上结果表明，利用本文方法生成的湍流，流场中的湍流强度和积分尺度均满足相应的幂律关系，且具有较为良好的均匀各向同性，能够满足湍流相关研究的要求.

3.2 湍流功率谱

脉动风速功率谱反映了脉动风能量在频域上的分布情况，本节采用周期图谱法计算脉动风速的功率谱，为确保湍流风场的质量，同与大气边界层湍流风速谱吻合的von Kármán谱进行拟合比较，von Kármán谱方程为：

\frac{f{S}_{\mathrm{u}}(f)}{{\sigma }_{\mathrm{u}}^{\mathrm{2}}}=\frac{\mathrm{4}\left(\frac{f{L}_{\mathrm{u}}}{U}\right)}{{\left[\mathrm{1}+\mathrm{70.8}{\left(\frac{f{L}_{\mathrm{u}}}{U}\right)}^{\mathrm{2}}\right]}^{\mathrm{5}/\mathrm{6}}}（15）

\begin{array}{l}\frac{f{S}_i(f)}{{\sigma }_i^{\mathrm{2}}}=\frac{\mathrm{4}\left(\frac{f{L}_i}{U}\right)\left[\mathrm{1}+\mathrm{755}{\left(\frac{f{L}_i}{U}\right)}^{\mathrm{2}}\right]}{{\left[\mathrm{1}+\mathrm{238}{\left(\frac{f{L}_i}{U}\right)}^{\mathrm{2}}\right]}^{\mathrm{11}/\mathrm{6}}}\\ (i=\mathrm{v},\mathrm{w})\end{array}（16）

式中：f为频率；U为平均风速；S_i(f)表示脉动风速i的功率谱；L_i表示脉动风速i的湍流积分尺度；{\sigma }_i^{\mathrm{2}}是脉动风速i的方差. 功率谱曲线的面积大小反映了该点处脉动风速的总功率，频率的高低与涡旋尺寸的大小呈负相关. 图16展示了流向脉动风速功率谱，由图可以看出，随着入口距离增加，脉动风速功率谱高频部分的功率逐渐下降，低频部分的功率略有升高，曲线围成的总面积不断减小，这表明随着入口距离增加，能量从小涡逐渐传递到大涡，同时脉动风的能量不断降低（能量耗散），且能量耗散逐渐减慢；这与3.1节中湍流强度的衰减和湍流积分尺度的增加相对应. 图17展示了流向9 m位置的截面中心处流向、展向和垂直方向的脉动风速功率谱. 由图17可知，功率谱在低频区与von Kármán谱拟合较好，高频部分出现了偏离. 这是因为受网格尺寸的影响，小于滤波尺度的高频涡旋无法被网格捕捉，只能用亚格子模型解析，这使得高频部分无法准确模拟，因此亚格子模型的求解精度和网格尺寸决定了高频部分的求解精度. 由于在实际工程中，通常主要考虑大尺度涡旋（低频部分）对结构气动力特性的影响^［28］，本文与von Kármán谱拟合较好的频段为0.01~0.9 Hz，后续研究中模型的理论涡脱频率f_vs=0.44 Hz，因此该频段能够满足后续研究需求.

4 圆柱绕流研究

圆柱绕流是流体力学的一个经典问题，工程结构中如桥塔、拉索、输电线等都具有一定的圆形截面特征，亚临界雷诺数区下的圆柱绕流，因具有丰富的流场特征而被广泛研究. 周强等^［29］、战庆亮等^［30］及Ma等^［31］使用LES方法对经典亚临界雷诺数Re=3 900的圆柱展开了均匀来流下绕流特性的研究. 为研究湍流下圆柱的绕流特性，利用虚拟格栅产生的湍流风场，进行格栅湍流下Re=3 900的圆柱绕流研究. 将圆柱放置在距离入口9 m处，该位置处的流向平均风速U_x =V=0.352 1 m/s，稳定区的流向平均速度约为入口风速的3.50%，计算域如图18所示.

4.1 数值模型设置

首层最小网格高度取为d=0.03D，以确保圆柱近壁面网格y⁺ ≤1. Jing等^［32］对Re=3 900的圆柱绕流网格划分方式进行了详细研究，认为圆周节点数为160时可满足精度要求，因此本研究的圆周节点数设为160. 全域中最大网格尺寸不超过0.04 m，圆柱附近网格采用O型切分，流体运动黏度v=1.460 8×10^-5 m²/s，为确保雷诺数Re=3 900，圆柱直径D选为0.161 8 m，雷诺数计算公式为Re=VD/v，计算域的阻塞率为4%. 为减小引入实体模型对流场的影响，流场区与绕流加密区采用循环边界条件（cyclicAMI），两交界面单元间的物理量采用插值过渡. 网格划分如图19所示，圆柱表面采用无滑移边界，其余壁面的边界条件均设置为滑移边界.

4.2 数值验证

为验证该数值方法的准确性，基于4.1节数值模型设置，进行均匀来流下圆柱绕流计算，将计算结果与已有研究结果^［29-33］进行对比. 表2展示了均匀来流和湍流中圆柱的归一化回流长度〈L_r〉/D、斯脱罗哈数Sr、分离角〈θ_sep〉（分离点与前驻点所对应的圆心角）、平均阻力系数〈C_d〉、归一化最小平均速度-u_min/V、升力系数均方根（C₁）_rms和阻力系数的平均脉动值C_{\mathrm{d}}^{\text{'}}的计算结果. 由表2可知，均匀来流下的计算结果与已有文献［29-33］的差异很小，验证了本文的数值方法是准确和有效的.

4.3 结果和讨论

图20展示了不同来流下涡脱稳定后的升阻力系数时程图，升力系数和阻力系数分别定义为C_{\mathrm{l}}=\mathrm{2}{F}_y/\left(\rho V^{\mathrm{2}}D\right)和C_{\mathrm{d}}=\mathrm{2}{F}_x/\left(\rho V^{\mathrm{2}}D\right)，式中：F_y和F_x分别表示升力（垂直流向）和阻力（顺流向）. 由图20可知，均匀流下的阻力系数基本在定值上下脉动，阻力系数的均值为0.95；升力系数脉动幅值变化较为平缓，升力系数的均值为零. 而湍流下的升阻力系数变化均体现出很强的不规则性，并且脉动变得更加强烈，湍流下阻力系数的均值为0.94，相比均匀流略微减小. Cheung等^［34］在湍流强度对圆柱气动力特性影响的研究中发现，在亚临界雷诺数区，圆柱阻力系数随湍流强度增大而减小，阻力系数的标准差和升力系数的均方根分别反映了各自的平均脉动程度，计算发现湍流下阻力系数的标准差大约是均匀流时的5.6倍，湍流下升力系数的均方根约为均匀流的2.4倍.

图21为通过快速傅里叶变换得到的升力系数频谱图. 圆柱的升力系数波动是由圆柱尾部涡旋脱落形成，通常用斯脱罗哈数（Strouhal number）Sr来描述均匀流下涡旋脱落现象，圆形断面的斯脱罗哈数约为0.20，斯脱罗哈数公式为Sr=f_vsD/V. 式中：f_vs为涡旋脱落频率，为升力系数频谱图中最大振幅处的频率；D为物体在来流平面上的特征长度. 对横坐标进行无量纲化处理，可以看到，均匀流下的升力频谱曲线具有明显的主峰，频谱曲线在达到涡脱频率时幅值急剧增加，其余频率处幅值趋于零，幅值分布集中，该位置处Sr=0.22. 而湍流下频谱曲线并未出现明显主峰，幅值在频率0.2附近离散分布，湍流下的Sr=0.19，湍流涡脱振幅约为均匀流的2.3倍. 这表明湍流下涡漩脱落得更慢，湍流对圆柱表面的涡旋脱落产生了一定的抑制，且湍流使圆柱的涡漩脱落频率不再稳定，尾流涡旋呈现多频脱落，湍流下尾流的涡脱能量更强. Bruun等^［35］发现的来流湍流会抑制尾流涡旋脱落. 这可能是由于不同的来流频率影响了圆柱的涡脱频率^［36］.

对于亚临界雷诺数下的圆柱，其阻力主要源于边界层分离而产生的压差阻力. 压力系数定义为C_{\mathrm{p}}=\mathrm{2}\left(〈p〉-p_{\mathrm{\infty }}\right)/\left(\rho V_{\infty }^{\mathrm{2}}\right)，式中：p为压力，\rho为流体密度，V_{\mathrm{\infty }}为远场速度，p_{\mathrm{\infty }}为远场压力. 图22展示了圆柱表面平均压力分布，θ表示圆柱表面某点相对于前驻点的圆心角.由图22可知，本文均匀流的结果与已有文献数据具有良好的一致性，具体为，平均压力系数随角度增加先减小后增加，在θ=90°附近增长逐渐减弱，随后趋于稳定. 这是因为当流体流经圆柱时，受圆柱几何阻碍作用，驻点处流速降至零，该位置压力达最大值. 对于不可压缩流体，在压力作用下，流动于驻点处发生分离并沿圆柱表面开始发展；此阶段压力梯度与流动方向一致，流动速度逐渐提升，压力部分转化为流体动能，部分用于克服黏滞力而损耗，直至某一位置速度达峰值，压力同步降至最小值. 此后流体进入减速阶段，该过程中惯性力需克服黏滞力和逆压差，随着动能持续耗散，压力逐渐增加. 当流动停滞时，惯性力无法平衡黏滞力和逆压差，在逆压差作用下，流动方向发生逆转并伴随回流现象，进而引发边界层分离. 回流区域内，压力受回流涡旋影响而最终趋于稳定.对比湍流和均匀流中圆柱表面的平均压力系数分布可知，湍流工况下，圆柱驻点处的压力系数低于均匀流，而在θ为45°~80°区间内，其压力系数略大于均匀流. 平均背压系数定义为θ=180°处的平均压力系数，由图22可知，湍流下的平均背压系数绝对值略大于均匀流. 湍流对圆柱表面平均压力系数的影响主要集中在驻点、平均压力系数最低点和背压点三个区域，但总体而言，在试验雷诺数范围内，湍流对圆柱表面的平均压力系数分布的影响较小.

涡量是描述涡旋强度的重要参数，即速度的旋度，表达式为\Omega =\nabla \times \vec{V}=\mathrm{2}\omega，其中ω为旋转角速度. 图23给出了不同来流下圆柱表面涡量分布. 可以看出，湍流下的分离角相对于均匀流出现明显的后移. 这是因为边界层分离需要达到一定的逆压梯度，而湍流的流动总是伴随不同层之间的掺混，各层间的动量交换导致不同层间的流体相互拖动，这使得边界层可以承受更大的逆压梯度而不分离，进而造成了湍流流动分离延后，分离点延后是造成平均阻力系数减小的重要原因.

回流区的长度定义为圆柱θ=180°壁面位置与尾流中心线上平均流向速度为零处的直线距离. 图24为尾流中心线上平均流向速度分布. 由图可知，均匀流的结果与已有文献［29-30，32］数据具有良好一致性. 在圆柱后方存在速度为负值的区域，即圆柱尾流回流区；此区域内速度绝对值先增大后减小，其最大值所在位置定义为回流中心. 此外，湍流下回流中心的速度绝对值和回流长度更小，且回流中心更接近圆柱壁面. 在回流区外、速度稳定之前的区段，湍流下平均流向速度大于均匀流. 进一步分析表明，湍流下背压点附近速度梯度更大，这一现象源于回流中心更靠近壁面——尽管均匀流下回流中心的压强更低，但其对壁面的吸附作用弱于湍流. 此外，回流区外更高的速度意味着湍流下涡旋脱落的能量更高，由于尾流涡旋脱落是造成升力系数脉动的原因，这一特征在一定程度上阐明了湍流下升力系数频谱振幅增大的原因.

为更直观了解湍流对圆柱流场的影响，图25给出了不同来流下流场的时间平均流线图. 由图可知，流体流经圆柱时，在驻点位置处流速降至零，随角度的增加，速度先增大后减小，分离点后的尾流区域速度较为稳定. 由于圆截面的对称性，图中物理量也沿流向对称分布. 图中尾流区内两个对称的涡旋是由于边界层速度逆转形成的回流涡，对比均匀流和格栅湍流，可以发现均匀流的回流涡更大，并且沿流向的长度差异更显著，同时边界层外速度分布更加均匀. 产生这种现象的原因是，湍流的不稳定造成流体更易相互混合，抑制了边界层发展，且离圆柱越远抑制作用越强，因此格栅湍流下回流区更小. 由于均匀流下流体间不相互干扰，因此边界层得到充分发展，尾流涡旋更宽更长. 此时涡旋中心具有更高的速度，压力更小. 而湍流下回流涡旋中心更靠近壁面，因此涡旋的吸附力更强，这可能是导致湍流下平均背压系数更小的原因.

图26展示了在x-z平面上圆柱附近的时间平均流线图. 可以发现，均匀流下圆柱表面存在小尺度的涡旋，这可能是由于无滑移壁面边界限制了圆柱底部流体在垂直方向上的运动，从而产生了涡旋，因此圆柱高度的选取十分重要. 由于本文圆柱在垂直方向高度为18.5D，在一定程度上减弱了端部壁面对计算结果造成的影响，从与已有文献数据的对比结果来看，这部分影响可以忽略. 图26中来流为湍流时，对涡旋的抑制效果显著，圆柱壁面的涡旋基本消失，流线更加平顺.

5 结 论

本文提出了一种在计算域入口设置虚拟格栅函数生成格栅湍流的数值方法，该方法具有较高的模拟精度，且避免了数值模拟方法生成格栅湍流时计算成本高的问题. 利用该方法，获得了流场内湍流参数的空间变化规律，并进行了Re=3 900的圆柱绕流研究，分析了湍流对圆柱气动力与流场特性的影响机理. 主要结论如下：

1）在虚拟格栅生成的湍流风场中，湍流强度和湍流积分尺度满足均匀各向同性. 湍流强度的衰减指数n=1.12，积分尺度增长率为k≈0.7，该方法得到的湍流流场可以较好地满足传统实体格栅湍流风场中的幂律关系.

2）虚拟格栅生成的湍流风场中脉动风速功率谱均符合von Kármán谱.

3）湍流会导致圆柱的平均阻力系数和回流长度减小，但升阻力系数的脉动振幅增大，分离点延后. 同时湍流会抑制尾涡的发展和脱离，导致涡旋呈多频脱落，且涡脱能量更强.

以上结果表明，采用虚拟格栅生成湍流的方法是有效、可行的，该方法可以用于湍流下圆柱绕流研究. 但这种方法仍存在一定缺陷，如无法直接定义下游生成的湍流特性，因此需要进行大量流场调试，从而获得目标流场. 此外，研究中尚未给出下游湍流特性与输入参数（如速度、格栅宽度、格栅数量等）间的关系，没有定量分析湍流特性对圆柱气动力及流场结构的影响，因此，后续有必要进一步研究.

参考文献

[1] 孙大奇，朱颖，刘晓光，等 ．基于非概率集合理论的屋面板风致疲劳寿命估计［J］．工程力学，2021，38（4）：211-220．

[2] SUN D Q， ZHU Y， LIU X G，et al ．Fatigue life analysis of panels under wind load based on non-probabilistic set theory［J］．Engineering Mechanics，2021，38（4）：211-220．（in Chinese）

[3] LI M， LI M S， SUN Y G ．Effects of turbulence integral scale on the buffeting response of a long-span suspension bridge［J］．Journal of Sound and Vibration，2021，490：115721．

[4] VICKERY B J ．Fluctuating lift and drag on a long cylinder of square cross-section in a smooth and in a turbulent stream［J］．Journal of Fluid Mechanics，1966，25（3）：481-494．

[5] KAWATANI M， KIM H， UEJIMA H，et al ．Effects of turbulent flows on vortex-induced oscillation of bridge girders with basic sections［J］. Journal of Wind Engineering and Industrial Aerodynamics，1993，49（1/2/3）：477-486．

[6] 白桦，何晗欣，刘健新，等 ．格栅紊流风特性参数模拟规律研究［J］．振动与冲击， 2016，35（22）：209-214．

[7] BAI H， HE H X， LIU J X，et al ．Wind characteristic parameters in grille turbulent flow［J］．Journal of Vibration and Shock，2016，35（22）： 209-214．（in Chinese）

[8] SIMMONS L F G， SALTER C. Experimental investigation and analysis of the velocity variations in turbulent flow［J］. Proceedings of the Royal Society of London．Series A，Containing Papers of a Mathematical and Physical Character，1934，145（854）： 212-234．

[9] CHANG Y， ZHAO L， CHEN L，et al ．Effect of free stream turbulence in critical Reynolds number regime （1.6∶105-6.1∶ 105） on flow around circular cylinder［J］．Physics of Fluids，2022，34（11）：115126．

[10] MARYAMI R， ALI S A S， AZARPEYVAND M，et al ．Turbulent flow interaction with a circular cylinder［C］//25th AIAA/CEAS Aeroacoustics Conference. 20-23 May 2019，Delft，The Netherlands．Reston，Virginia：AIAA，2019：2503．

[11] KANG H S， CHESTER S， MENEVEAU C ．Decaying turbulence in an active-grid-generated flow and comparisons with large-eddy simulation［J］．Journal of Fluid Mechanics，2003，480：129-160．

[12] 李春光，陈政清，张记 ．主动格栅紊流场对典型主梁颤振导数影响的研究［J］．实验力学，2016，31（1）：75-86．

[13] LI C G， CHEN Z Q， ZHANG J ．On the influence of active grid turbulence field on flutter derivatives of typical main girder［J］．Journal of Experimental Mechanics，2016，31（1）：75-86．（in Chinese）

[14] 胡尚瑜，李秋胜，张明 ．低矮建筑标模风荷载的主动湍流模拟试验研究［J］．实验流体力学，2020，34（4）：22-29．

[15] HU S Y， LI Q S， ZHANG M. Active turbulence simulation study of wind loads on standard low-rise building［J］. Journal of Experiments in Fluid Mechanics，2020，34（4）：22-29．（in Chinese）

[16] 周桐，杨庆山，闫渤文，等 ．大气边界层大涡模拟入口湍流生成方法综述［J］．工程力学，2020， 37（5）：15-25．

[17] ZHOU T， YANG Q S， YAN B W，et al ．Review of inflow turbulence generation methods with large eddy simulation for atmospheric boundary layer［J］．Engineering Mechanics， 2020，37（5）：15-25．（in Chinese）

[18] PHUC P V， NOZU T， KIKUCHI H，et al ．Wind pressure distributions on buildings using the coherent structure Smagorinsky model for LES［J］．Computation，2018，6（2）：32．

[19] 周桐，闫渤文，杨庆山，等 ．大气边界层大涡模拟入口湍流生成方法研究［J］．工程力学，2020，37（7）：68-76．

[20] ZHOU T， YAN B W， YANG Q S，et al ．Study of inflow turbulence generation methods with large eddy simulation for atmospheric boundary layer［J］．Engineering Mechanics，2020，37（7）：68-76．（in Chinese）

[21] YU J H， LI M S， STATHOPOULOS T ．Strategies for modeling homogeneous isotropic turbulence and investigation of spatially correlated aerodynamic forces on a stationary model［J］．Journal of Fluids and Structures，2019，90：43-56．

[22] 邱旭，王昊，李峰，等 ．格栅湍流场湍流特性的数值模拟［J］．沈阳建筑大学学报（自然科学版），2016，32（5）：863-870．

[23] QIU X， WANG H， LI F，et al ．Numerical simulation of wind characters of grid-generated wind field［J］．Journal of Shenyang Jianzhu University （Natural Science），2016，32（5）：863-870．（in Chinese）

[24] 胡晓兵，杨易 ．基于NSRFG方法的标准地貌风场大涡模拟研究［J］．工程力学，2020，37（9）：112-122．

[25] HU X B， YANG Y ．Research on NSRFG-based LES simulation for standard wind terrains［J］．Engineering Mechanics，2020，37（9）：112-122．（in Chinese）

[26] 王天，孙东，郭启龙，等 ．可压缩边界层的入口合成湍流生成方法［J］．力学学报， 2024， 56（1）： 45-57．

[27] WANG T， SUN D， GUO Q L，et al ．Inlet synthetic turbulence generation method for compressible boundary layer［J］．Chinese Journal of Theoretical and Applied Mechanics，2024，56（1）：45-57．（in Chinese）

[28] ROCCHIO B， MARIOTTI A， SALVETTI M V ．Flow around a 5∶1 rectangular cylinder：effects of upstream-edge rounding［J］．Journal of Wind Engineering and Industrial Aerodynamics，2020，204：104237．

[29] ABOSHOSHA H， ELSHAER A， BITSUAMLAK G T，et al ．Consistent inflow turbulence generator for LES evaluation of wind-induced responses for tall buildings［J］. Journal of Wind Engineering and Industrial Aerodynamics，2015，142：198-216．

[30] LUO Y， LIU H J， XUE H L，et al. Large-eddy simulation evaluation of wind loads on a high-rise building based on the multiscale synthetic eddy method［J］. Advances in Structural Engineering，2019，22（4）：997-1006．

[31] REN H H， LAIMA S， CHEN W L，et al ．Numerical simulation and prediction of spatial wind field under complex terrain［J］．Journal of Wind Engineering and Industrial Aerodynamics，2018，180：49-65．

[32] JING H M， LI W K， SU Y，et al ．Numerical study of wind characteristics at a long-span bridge site in mountain valley［J］．Physics of Fluids，2024，36（3）： 035131．

[33] 李启，杨庆山，朱伟亮 ．湍流入口条件下建筑非定常风场的大涡模拟［J］．工程力学，2012，29（12）：274-280．

[34] LI Q， YANG Q S， ZHU W L ．Large eddy simulation of unsteady wind field around building using turbulent inflow［J］. Engineering Mechanics，2012，29（12）：274-280．（in Chinese）

[35] MOHAMED M S， LARUE J C ．The decay power law in grid-generated turbulence［J］．Journal of Fluid Mechanics，1990，219：195-214．

[36] TAYLOR G I ．Statistical theory of turbulence-II［J］. Proceedings of the Royal Society of London. Series A-Mathematical and Physical Sciences，1935，151（873）：444-454．

[37] FLAY R G J， STEVENSON D C ．Integral length scales in strong winds below 20 m［J］. Journal of Wind Engineering and Industrial Aerodynamics，1988，28（1/2/3）：21-30．

[38] 张伟峰，张志田，张显雄，等 ．节段模型气动导纳的数值模拟与试验［J］．中国公路学报，2018，31（6）：207-216．

[39] ZHANG W F， ZHANG Z T， ZHANG X X，et al ．Experimental and numerical investigations on aerodynamic admittances of section models［J］．China Journal of Highway and Transport，2018，31（6）：207-216．（in Chinese）

[40] 周强， 曹曙阳， 周志勇 ．亚临界雷诺数下圆柱体尾流结构的数值模拟［J］．同济大学学报（自然科学版），2013，41（1）：33-38．

[41] ZHOU Q， CAO S Y， ZHOU Z Y ．Numerical studies of wake characteristics on a circular cylinder at sub-critical Reynolds number［J］．Journal of Tongji University （Natural Science），2013，41（1）： 33-38．（in Chinese）

[42] 战庆亮，周志勇，葛耀君 ． Re=3 900圆柱绕流的三维大涡模拟［J］．哈尔滨工业大学学报，2015，47（12）：75-79．

[43] ZHAN Q L， ZHOU Z Y， GE Y J ．3-Dimensional large eddy simulation of circular cylinder at Re=3 900［J］．Journal of Harbin Institute of Technology，2015，47（12）：75-79．（in Chinese）

[44] MA X， KARAMANOS G S， KARNIADAKIS G E ．Dynamics and low-dimensionality of a turbulent near wake［J］．Journal of Fluid Mechanics，2000，410：29-65．

[45] JING H， ZHANG J， LIU Q， et al. LES study of flow field and aerodynamic forces on a circular cylinder at Re=3 900 with focus on grid resolution［J］. Wind and Structures， 2023， 36（3）：175-200.

[46] NORBERG C ．An experimental investigation of the flow around a circular cylinder：influence of aspect ratio［J］．Journal of Fluid Mechanics，1994，258：287-316．

[47] CHEUNG J C K， MELBOURNE W H ．Turbulence effects on some aerodynamic parameters of a circular cylinder at supercritical numbers［J］．Journal of Wind Engineering and Industrial Aerodynamics，1983，14（1/2/3）：399-410．

[48] BRUUN H H， DAVIES P O A L ．An experimental investigation of the unsteady pressure forces on a circular cylinder in a turbulent cross flow［J］．Journal of Sound and Vibration， 1975，40（4）： 535-559．

[49] 于春放，靖洪淼，王仰雪，等 ．顺向正弦来流条件下圆柱气动力和绕流流场数值模拟研究［J］．振动与冲击，2024，43（6）：216-224．

[50] YU C F， JING H M， WANG Y X，et al ．Numerical study on the aerodynamic forces and flow around a circular cylinder under streamwise sinusoidal inflow condition［J］．Journal of Vibration and Shock，2024，43（6）：216-224．（in Chinese）

基金资助

国家自然科学基金资助项目（52208494）; 河北省自然科学基金资助项目（E2024210037,E2021210063）;河北省自然科学基金创新研究群体项目（E2022210078）; 河北省高端人才项目（[2019]63号）