对于非定常分离流动，雷诺平均模拟(RANS)[2]无法准确捕捉湍流涡的瞬态脉动，而大涡模拟(LES)[3]则由于对近壁面网格分辨率的严格要求导致其在工程应用中的计算成本太高。为此，Menter等[4-5]提出了尺度自适应模拟(SAS)方法，其思路是在传统RANS模型的基础上引入冯卡门尺度作为湍流模型的第二特征长度尺度。该方法虽然基于RANS方程，但在大分离区域却能表现出类似LES的性能，原因在于冯卡门尺度的构造中包含了速度场的高阶导数，能够自适应地随着当地湍流脉动大小变化，从而在小尺度的高频脉动区域降低湍流模型的涡黏性及耗散性[6]。同时，由于冯卡门尺度的构造中不包含网格尺度，SAS模型能够实现从RANS到LES区域的光滑过渡，避免了大多数RANS/LES混合方法所面临的网格敏感性相关问题[7-8]。然而，冯卡门尺度中高阶导数的离散仍然受到网格分辨率的影响。初步研究发现，不同网格分辨率下冯卡门尺度存在差异，因此有必要开展SAS的网格尺度关联分析，为SAS在复杂工程问题中的应用提供网格划分指导[9]。
对于雷诺数3 900的圆柱绕流，Norberg[10]实验测量了圆柱表面的压力分布，Ong和Wallace[11]利用热线风速计(HWA)精确测量了远壁面回流区的时均速度和雷诺应力分布，Lourenco和Shih[12]则使用粒子图像测速(PIV)方法获得了近壁面回流区的湍流统计数据，为数值模拟提供了验证基础[13-18]。SAS方法也应用于圆柱绕流的计算研究。最早的两方程SST-SAS模型由Menter和Egorov[14]提出，随后Shim等[19]采用该模型对雷诺数3 900的圆柱绕流进行了数值模拟并与Parnaudeu等[20]的LES和PIV数据进行比较，研究发现该模型预测的近壁面回流区二阶湍流统计量与实验数据有差异。在进一步的开发和验证过程中，Egorov和Menter[5]提出了新的两方程SST-SAS模型，目前鲜有文献针对新模型展开网格尺度关联分析。
本文的研究目的是对新的两方程SST-SAS模型进行网格尺度关联性分析和评估，研究对象是雷诺数3 900的圆柱绕流，对比目前广泛应用的分离涡模拟(DES)[21]和实验数据[12, 20]，分析不同网格分辨率和展向计算域对尾迹区一阶、二阶湍流统计量及冯卡门尺度的影响。
1 湍流模型 1.1 SST-SAS模型 该模型的构造方法是在SST模型的ω方程中添加源项QSAS：
 | (1) |
式中:Ui为速度矢量；?和β为常数系数；ρ为密度; k为湍动能; ω为比耗散率; Pk为k方程生成项; β*为k方程耗散项系数; Ui为i方向速度; σω为ω方程黏性项系数；t为时间; S为应变率; μ为层流黏性系数; σk为k方程黏性项系数; μt为湍流黏性系数; F1为混合函数；QSAS为
 | (2) |
 | (3) |
式中: ζ2=3.51；σΦ=2/3；C=2; κ为冯卡门常数; Φ为模型常数; L为湍流模型的第一特征长度尺度; LvK为冯卡门尺度，是湍流模型的第二特征长度尺度，它能够自适应地随当地求解的湍流涡结构而变化。
 | (4) |
式中: U为速度矢量; U、V、W为三方向速度分量。
为了提供合适的高波数耗散，需要限制LvK的最小值以防止SAS的涡黏性系数小于LES的亚格子黏性系数。
 | (5) |
式中:CS为模型常数；Δ为网格尺度。
以上为新版本的SST-SAS模型，相比于早期版本的区别之一是加入了高波数衰减限制器，另一个区别则是采用湍流尺度之比的平方(L/Lvk)2代替了原有的线性比形式。
1.2 SST-DES模型 SST-DES在SST模型的基础上修改了k方程源项耗散项的系数。
 | (6) |
式中：
 | (7) |
其中：CDES为模型常数。
1.3 控制方程和离散方法 控制方程为笛卡儿坐标系下积分形式的三维可压缩Navier-Stoker方程，基于单元中心有限体积法的自研求解器。无黏通量采用Roe格式进行空间离散，为了获得高阶的数值解，对计算界面两侧的变量进行五阶WENO重构。黏性通量采用四阶中心差分格式，时间离散则选择二阶后向差分格式的双时间步法。
2 计算网格 表 1列出了所有的计算状态。为了研究网格分辨率的影响，设置了A、B、C三套粗、中、细网格；为了研究展向计算域的影响，在中等网格的基础上设置了B2、B3、B4 3个不同展向长度(Lz/D，D为圆柱直径)并保持展向网格间距(Δz/D)不变；为了比较不同湍流模型的性能，采用DES(A1、B1、C1)进行对比研究，数值格式与SAS完全一致。
表 1 计算状态 Table 1 Computational setup
| 编号 | 湍流模型 | 网格数量 | Lz/D | Δz/D | Δt×U∞/D |
| A1 | DES | 137 137 31 | π | 0.105 | 0.01 |
| A2 | SAS | 137 137 31 | π | 0.105 | 0.01 |
| B1 | DES | 193 193 48 | π | 0.067 | 0.005 |
| B2 | SAS | 193 193 24 | 0.5π | 0.067 | 0.005 |
| B3 | SAS | 193 193 48 | π | 0.067 | 0.005 |
| B4 | SAS | 193 193 96 | 2π | 0.067 | 0.005 |
| C1 | DES | 249 249 61 | π | 0.052 | 0.005 |
| C2 | SAS | 249 249 61 | π | 0.052 | 0.005 |
表选项
XY平面内的中等网格划分如图 1所示，外边界设定为压力远场边界，圆柱表面采用无滑移壁面，壁面第一层网格高度为5×10-4D，假设流动在展向具有周期性。粗网格采用0.01的无量纲时间步长(D/U∞，U∞为来流速度)，中等和细网格采用更小的0.005时间步长。迭代20个涡脱落周期后达到稳定，湍流统计量在60个涡脱落周期内进行时间平均和展向平均。
 |
| 图 1 XY平面的中等网格 Fig. 1 Medium grid in XY-plane |
| 图选项 |
3 结果分析 3.1 时均速度分布 图 2显示了不同展向长度和网格分辨率下时均流向速度沿尾迹区中心线的分布规律。从图 2(a)可以看出，在相同展向网格分辨率下不同的展向长度对SAS的计算结果影响不大，也证明了大多数文献所采用的Lz/D=π展向长度是合适的。从图 2(b)可以看出，在相同的网格分辨率下SAS预测的回流区长度总体来说小于DES。SAS方法表现出较好的网格收敛性：粗网格过低地预测了回流区长度，中等网格得到的回流区长度与Lourenco & Shih的实验结果相接近，在此基础上继续加密网格对计算结果的改善不明显。对于DES方法，回流区长度随网格加密而不断增加：粗网格下的尾迹区中心线速度分布与Lourenco & Shih的实验结果相接近，中等网格下的尾迹区中心线速度分布逐渐趋近于Parnaudeau等的实验结果，而细网格预测的回流区长度偏大。前期研究表明[22]，Lourenco & Shih开展的早期实验可能受到环境扰动的影响，导致剪切层过早转捩。由此推测，SAS预测的剪切层失稳早于DES，网格分辨率对SAS的剪切层失稳特征影响较小，而对DES则影响较大。
 |
| 图 2 时均流向速度沿尾迹区中心线的分布规律 Fig. 2 Distribution law of streamwise mean velocity along centerline of wake region |
| 图选项 |
图 3给出了近壁面回流区内沿轴线不同站位的时均流向和法向速度分布。图 3~图 4中的曲线和符号含义与图 2(b)完全相同。图 3(a)中不同站位处SAS预测的时均流向速度均呈V型分布，回流速度随着网格加密逐渐增大，趋近于Lourenco & Shih的实验结果。粗网格下，DES预测的时均流向速度也呈V型分布，趋近于Lourenco & Shih的实验结果；随着网格加密，由x/D=1.06处的U型分布逐渐过渡到x/D=2.02处的V型分布，DES得到的剪切层失稳延迟、回流速度增大；中等网格下DES的计算结果与Parnaudeau等的实验数据非常接近。图 3(b)中时均法向速度的分布规律与上述流向速度的分析基本一致。由于时均速度分布与速度脉动直接相关，下面将进一步讨论圆柱绕流尾迹区雷诺应力的分布规律。
 |
| 图 3 近壁面回流区内沿轴线不同站位的时均速度分布 Fig. 3 Mean velocity profiles at different locations along axis in near wake region |
| 图选项 |
 |
| 图 4 近壁面回流区内沿轴线不同站位的时均雷诺应力分布 Fig. 4 Mean Reynolds stress profiles at different locations in near wake region |
| 图选项 |
3.2 时均雷诺应力分布 图 4为近壁面回流区内沿轴线不同站位的时均流向、剪切及法向雷诺应力分布。近壁面的x/D＝1.06站位处，相同网格分辨率下SAS比DES预测的流向和法向雷诺应力峰值更高，尤其是在粗网格下，由此导致求解湍动能偏大；随着网格加密，2种湍流模型预测的雷诺应力峰值降低，SAS的计算结果仍高于实验数据而DES更接近Parnaudeau等的实验结果。x/D=1.54站位处，网格加密后SAS预测的流向雷诺应力与实验数据吻合较好，但法向雷诺应力峰值仍偏高，造成剪切应力偏离实验结果；中等网格下DES得到的雷诺应力分布与Parnaudeau等的实验数据非常吻合，而细网格下DES得到的雷诺应力峰值略低于实验结果。在下游的x/D=2.02站位处，对于SAS模型，不同网格分辨率之间的差异几乎消失，流向雷诺应力峰值低于实验结果，法向雷诺应力与实验数据相接近；对于DES方法，网格加密后的流向雷诺应力与实验结果比较吻合，而流向雷诺应力峰值低于实验数据。
3.3 瞬态涡结构和冯卡门尺度 图 5对比了中等网格分辨率下DES和SAS的尾迹区流场，通过Q准则显示瞬态涡结构，定义为
 | (8) |
 |
| 图 5 Q准则显示的瞬态涡结构 Fig. 5 Instantaneous vortex structure plotted by iso-surface of Q-criterion |
| 图选项 |
式中：S和Ω分别为涡量和应变率张量。
图 5中SAS和DES都能清晰预测圆柱绕流的边界层分离、自由剪切层的形成以及尾迹区涡结构的发展过程。其中，DES方法通过引入网格尺度可以捕捉到与当地网格同量级的大量小尺度脉动，而SAS模型则过滤了许多小尺度脉动而捕捉对流动起主导作用的大尺度涡结构。
SAS模型的上述特性因为在ω方程的附加源项QSAS中引入了冯卡门尺度LvK，因此图 6显示了LvK和QSAS在XY平面的瞬态分布规律。图 6中LvK可以根据当地湍流脉动尺度自适应动态调整，在涡结构内部逐渐减小，QSAS由于包含LvK倒数的平方而明显增大，从而导致SAS模型预测的涡粘性系数降低，展现出类似LES的性能。
 |
| 图 6 XY平面的QSAS和LvK的瞬态分布 Fig. 6 Instantaneous distributions of QSAS and LvK in XY-plane |
| 图选项 |
4 结论 本文针对雷诺数3 900的圆柱绕流对新版本的两方程SST-SAS模型展开了网格尺度关联分析，重点分析了网格分辨率和展向计算域对SAS计算结果的影响，为SAS在复杂工程问题中的应用提供网格划分指导。
1) ?研究发现，在相同的网格分辨率下，相比于DES方法，SAS预测的剪切层失稳更早，导致回流区长度小、流向和法向雷诺应力峰值高、求解湍动能偏大。网格分辨率对SAS的剪切层失稳特征影响较小，而对DES则影响较大。
2) ?随着网格加密，SAS模型计算的回流区长度和回流速度增大、雷诺应力峰值降低，逐渐趋近于Lourenco & Shih的实验结果。对于三维瞬态流动，SAS模型过滤了许多小尺度脉动而捕捉对流动起主导作用的大尺度涡结构，能够清晰预测圆柱绕流的尾迹区涡结构发展过程。
3) ?此外，在相同的网格分辨率下增大展向区域尺寸对SAS的计算结果影响很小，为展向具有周期性的流场研究提供了方便。
雷诺数3 900的圆柱绕流位于亚临界流动区域，转捩发生在分离的自由剪切层之中，本文采用的全湍流模拟准确性会受到制约。下一步工作将开展带边界层转捩的尺度自适应模拟，并研究时间步长对计算结果的影响。
参考文献
| [1] | 杜磊, 宁方飞. 高亚临界雷诺数圆柱绕流的尺度自适应模拟[J]. 力学学报, 2014, 46(4): 487-496. DU L, NING F F. Scale adaptive simulation of flows around a circular cylinder at high sub-critical Reynolds number[J]. Chinese Journal of Theoretical and Applied Mechanics, 2014, 46(4): 487-496. (in Chinese) |
| [2] | TRAVIN A, SHUR M, SPALART P, et al.On URANS solutions with LES-like behavior[C]//Congress on Computational Methods in Applied Sciences and Engineering, 2004. |
| [3] | FUREBY C. Towards the use of large eddy simulation in engineering[J]. Progress in Aerospace Sciences, 2008, 44(6): 381-396. DOI:10.1016/j.paerosci.2008.07.003 |
| [4] | MENTER F R, EGOROV Y. The scale-adaptive simulation method for unsteady turbulent flow predictions.Part 1:Theory and model description[J]. Flow, Turbulence and Combustion, 2010, 85(1): 113-138. DOI:10.1007/s10494-010-9264-5 |
| [5] | EGOROV Y, MENTER F R, LECHNER R, et al. The scale-adaptive simulation method for unsteady turbulent flow predictions.Part 2:Application to complex flows[J]. Flow, Turbulence and Combustion, 2010, 85(1): 139-165. DOI:10.1007/s10494-010-9265-4 |
| [6] | 李钊, 陈海昕, 张宇飞. 基于k-kL两方程湍流模式的尺度自适应模拟[J]. 工程力学, 2016, 33(12): 26-35. LI Z, CHEN H X, ZHANG Y F. Scale adaptive simulation based on a k-kL two-equation turbulence model[J]. Engineering Mechanics, 2016, 33(12): 26-35. (in Chinese) |
| [7] | 许常悦, 倪竹青, 孙智, 等. 尺度自适应模拟和大涡模拟的关联性分析[J]. 气体物理, 2018, 14(2): 49-58. XU C Y, NI Z Q, SUN Z, et al. Analysis of the relationships between scale-adaptive and large-eddy simulation[J]. Physics of Gases, 2018, 14(2): 49-58. (in Chinese) |
| [8] | 郑玮琳, 阎超. XY-SAS模型中不同网格尺度限制器的影响分析[J]. 北京航空航天大学学报, 2014, 40(12): 1725-1729. ZHENG W L, YAN C. Influence analysis on grid scale limiter of XY-SAS model[J]. Journal of Beijing University of Aeronautics and Astronautics, 2014, 40(12): 1725-1729. (in Chinese) |
| [9] | 杨振东, 谷正气. 基于尺度自适应模拟的汽车天窗风振噪声特性分析[J]. 机械工程学报, 2016, 12(52): 107-117. YANG Z D, GU Z Q. Analysis of car sunroof buffeting noise based on scale-adaptive simulation[J]. Journal of Mechanical Engineering, 2016, 12(52): 107-117. (in Chinese) |
| [10] | NORBERG C.Pressure forces on a circular cylinder in cross flow[C]//International Union of Theoretical and Applied Mechanics Symposium Bluff-Body Wakes, Dynamics and Instabilities.Berlin: Springer, 1993: 275-278. |
| [11] | ONG L, WALLACE J. The velocity field of the turbulent very near wake of a circular cylinder[J]. Experiments in Fluids, 1996, 20(6): 441-453. DOI:10.1007/BF00189383 |
| [12] | LOURENCO L M, SHIH C.Characteristics of the plane turbulent near wake of a circular cylinder, a particle image velocimetry study[Z].Private Communication, 1993. |
| [13] | KRAVCHENKO A G, MOIN P. Numerical studies of flow over a circular cylinder at ReD=3 900[J]. Physics of Fluids, 2000, 12(2): 403-417. DOI:10.1063/1.870318 |
| [14] | FRANKE J, FRANK W. Large eddy simulation of the flow past a circular cylinder at ReD=3 900[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2002, 90(10): 1191-1206. DOI:10.1016/S0167-6105(02)00232-5 |
| [15] | MANI A, MOIN P, WANG M. Computational study of optical distortions by separated shear layers and turbulent wakes[J]. Journal of Fluid Mechanics, 2009, 625(7): 273-298. |
| [16] | OUVRARD H, KOOBUS B, DERVIEUX A, et al. Classical and variational multiscale LES of the flow around a circular cylinder on unstructured grids[J]. Computers & Fluids, 2010, 39: 1083-1094. |
| [17] | WORNOM S, OUVRARD H, SALVETTI M V, et al. Variational multiscale large eddy simulations of the flow past a circular cylinder:Reynolds number effects[J]. Computers & Fluids, 2011, 47(1): 44-50. |
| [18] | DMITRY A L, IVAR S E, KJELL E R. Large-eddy simulation of the flow over a circular cylinder at Reynolds number 3900 using the OPENFOAM toolbox[J]. Flow, Turbulence and Combustion, 2012, 89(4): 491-518. DOI:10.1007/s10494-012-9405-0 |
| [19] | SHIM Y M, SHARMA R N, RICHARDS P J.Numerical study of the flow over a circular cylinder in the near wake at Reynolds number 3900[C]//39th AIAA Fluid Dynamics Conference.Reston: AIAA, 2013: 2009-4160. |
| [20] | PARNAUDEU P, CARLIER J, HEITZ D, et al. Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900[J]. Physics of Fluids, 2008, 20(8): 085101. DOI:10.1063/1.2957018 |
| [21] | STRELETS M.Detached eddy simulation of massively separated flows[C]//39th AIAA Aerospace Science Meeting and Exhibit.Reston: AIAA, 2001: 1-18. |
| [22] | XU C Y, SUN Z, ZHANG Y T, et al. Improvement of the scale-adaptive simulation technique based on a compensated strategy[J]. European Journal of Mechanics/B Fluids, 2020, 81: 1-14. DOI:10.1016/j.euromechflu.2020.01.002 |
