引言

方柱绕流作为一种非常典型的钝体绕流流动, 在航空航天、能源、建筑等领域中广泛存在, 其结构形式简单, 却蕴含了丰富且复杂的流体力学现象, 长久以来一直备受研究者关注.随着实验及数值模拟技术的不断进步, 国内外学者围绕该问题开展了大量工作, 取得了许多进展.

Lyn等^[1-2]用实验方法对Reynolds数为21 400的方柱绕流进行了较为细致的测量, 得到了柱体表面、两侧分离区、尾迹区的时均流场、Reynolds应力分布等, 该实验结果在后续的其他研究中作为对照被广泛引用.基于上述数据, Voke^[3], Rodi等^[4]陆续组织了两次研讨会, 采用大涡模拟方法对Reynolds数为21 400或22 000的方柱绕流进行计算, 分析对比不同研究小组的结果, 发现各小组的数据偏差较大, 没有一个让人完全满意, 计算域尺寸、网格分辨率及边界条件均对计算结果有很大影响, 而来流湍流度、数值格式的耗散等也可能是LES结果与实验数据相偏离的原因.尽管和RANS方法的对比表明LES结果更为可信^[5-7], 但对这一较高Reynolds数钝体绕流问题的准确预测显然还面临着诸多挑战. Murakami等^[8]和Sohankar等^[9]分别评估了几种亚格子模型对Reynolds数22 000的方柱绕流LES结果的影响, 前者发现Lagrange型动力模型的结果与实验数据吻合最好, 后者则认为一方程动力模型最优; Grigoriadis等^[10]采用浸没边界方法对同一问题进行大涡模拟, 发现足够的网格分辨率是提高预测精度的关键因素. Minguez等^[11], Brun等^[12]同时运用实验和LES方法对Reynolds数21 400的方柱绕流进行了更为深入的研究, 探讨了Karman涡街与自由剪切层K-H不稳定性的关系.考虑到过去一些算例存在着网格过于粗糙、计算周期过少及计算域不够大等问题, Trias等^[13]采用直接数值模拟方法, 扩大计算域并提高网格分辨率, 增加统计时间, 对Reynolds数22 000的方柱绕流进行了更精细的模拟, 结果与已有数据吻合得很好.

国内研究者对该问题也开展了不少工作, 童兵等^[14]对Reynolds数22 000的方柱进行了大涡模拟, 所得气动力和时均速度场与实验结果吻合较好; 邓小兵等^[15]采用拟压缩方法结合双时间步求解不可压LES方程模拟同一流动, 通过与实验结果的对比验证了该计算方案的可行性.谢志刚等^[16]采用有限体/有限元混合格式、非结构网格和LES方法对Reynolds数22 000的方柱绕流进行了计算, 所得结果都与实验数据符合得较好, 甚至优于前人的密网格结果.张海涛等^[17]使用了开源CFD软件OpenFOAM的两种动力模型对Reynolds数22 000的方柱进行计算, 对比发现一方程动力模型的结果与实验数据吻合得更好.唐鹏等^[18]运用两种RANS/LES混合方法在粗细两套网格上对同一方柱进行了计算, 发现DES方法的整体表现最优.

DNS方法虽然精确, 但需耗费大量的计算资源, 代价巨大, LES方法对该问题的预测结果还不能令人满意.本文使用自主发展的LES程序, 以远少于DNS的网格对来流Mach数M=0.3, Reynolds数Re_D=22 000的孤立方柱绕流进行模拟, 通过与已有实验和DNS数据的比较, 进一步评估LES方法对此类流动的预测能力.

1 数值方法

计算工作使用叶建^[19-20]发展的大涡模拟程序MPLES, 该程序能够处理任意界面匹配的多块结构化网格, 通过计算域分解和消息传递接口(message passing interface, MPI)实现并行, 具有很高的计算效率. MPLES求解Favre滤波的无量纲化可压缩N-S方程组, 其具体形式如下

$ \partial \bar{\rho }/\partial t+\partial \left( \bar{\rho }{{{\tilde{u}}}_{j}} \right)/\partial {{x}_{j}}=0 $

(1)

$ \partial \left( \bar{\rho }{{{\tilde{u}}}_{j}} \right)/\partial t+\partial (\bar{\rho }{{{\tilde{u}}}_{i}}{{{\tilde{u}}}_{j}}+\bar{p}{{\delta }_{ij}})/\partial {{x}_{j}}=\partial ({{{\tilde{\sigma }}}_{ij}}-{{\tau }_{ij}})/\partial {{x}_{j}} $

(2)

$ \begin{align} &\partial \left( \bar{\rho }\tilde{E} \right)/\partial t+\partial \left( \bar{\rho }\tilde{E}+\bar{p} \right){{{\tilde{u}}}_{j}}/\partial {{x}_{j}}= \\ &-\partial ({{{\tilde{q}}}_{j}}+{{Q}_{j}}/\left( \gamma - 1 \right)M_{0}^{2})/\partial {{x}_{j}}+\partial {{{\tilde{u}}}_{i}}({{{\tilde{\sigma }}}_{ij}}-{{\tau }_{ij}})/\partial {{x}_{j}} \\ \end{align} $

(3)

要封闭上述方程组, 必须对亚格子应力和亚格子热通量进行模化, MPLES采用如下形式的动力涡黏模型:

$ {{\tau }_{ij}}-{{\delta }_{ij}}{{\tau }_{kk}}/3=-2C\bar{\rho }{{\Delta }^{2}}|\tilde{S}|({{{\tilde{S}}}_{ij}}-{{\delta }_{ij}}{{{\tilde{S}}}_{kk}}/3) $

(4)

$ {{\tau }_{kk}}=2{{C}_{\text{I}}}\bar{\rho }{{\Delta }^{2}}|\tilde{S}{{|}^{2}} $

(5)

$ {{Q}_{j}}=-C{{\Delta }^{2}}\bar{\rho }|\tilde{S}|(\partial \tilde{T}/\partial {{x}_{j}})/P{{r}_{\text{T}}} $

(6)

其中的系数C_I, C和Pr_T通过公式动态计算

$ {{C}_{\text{I}}}=\left\langle {{L}_{kk}} \right\rangle /\left\langle 2{{\Delta }^{2}}\left( \hat{\bar{\rho }}{{\alpha }^{2}}|\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\tilde{S}}{{|}^{2}}-\widehat{\bar{\rho }|\tilde{S}{{|}^{2}}} \right) \right\rangle $

(7)

$ C=\langle ({{L}_{ij}}-{{\delta }_{ij}}{{L}_{kk}}/3){{M}_{ij}}\rangle /\langle {{M}_{kl}}{{M}_{kl}}\rangle $

(8)

$ P{{r}_{\text{T}}}=C\langle {{T}_{k}}{{T}_{k}}\rangle /\langle {{K}_{j}}{{T}_{j}}\rangle $

(9)

上述各项的具体含义参见文献[19].将封闭的方程组(1)~(3) 写成积分形式, 采用网格中心的有限体积法离散.对流项使用4阶反对称型的中心格式^[21], 黏性项使用2阶中心格式, 时间推进则是3阶3步的紧致Runge-Kutta方法.为消除因对流项中心型格式而在流场中出现的非物理振荡, 还使用11点的显式选择性滤波器对流场变量进行滤波处理^[22], 相关细节见文献[20].

2 算例设计

本文所用计算程序主要针对可压缩流动, 为了和实验及数值模拟结果对比, 计算的Reynolds数取22 000, Mach数取0.3.图 1给出了二维计算域及计算网格的示意图, 以方柱中心为坐标原点, x, y分别代表流向和横向, 二维计算域沿展向z拉伸得到三维计算域.以方柱的边长D为参考长度, 结合文献[3-4, 13], 计算域大小设定如下:入口位于x=-10处, 有效计算域出口位于x=15处, 其后有长度为5的缓冲区, 横向计算域高度为20, 展向宽度为π.

计算域沿流向网格单元数为420(包含出口缓冲区的20个单元), 横向为330, 展向为80, 总的网格单元数为1.044×10⁷.柱体附近的网格分布如图 1(b)所示, 壁面每条边长分布80个网格单元, 第1层网格的法向高度为0.003, 通过恰当布置网格点, 可以满足LES方法对网格分辨率的要求.

算例的边界条件设置如下:方柱表面为绝热无滑移壁面, 所有外边界均为远场, 采用特征边界条件处理, 展向则是周期性边界.计算依托国家超算广州中心的天河二号系统完成, 根据实际需要, 本算例使用了188个核心并行, 计算的时间步长取为6.0×10^-4, 初场设定为均匀场, 待流场演化足够的时间(2×10⁵步, 无量纲时间120) 后, 初场的影响早已消失, 流动进入统计定常状态, 此后按照每200步输出一个结果的频率输出瞬态流场.本文记录了约1.8×10⁵步的结果, 从中取出15个旋涡脱落周期的数据(约900个瞬态流场)用于统计, 将这些瞬态流场沿时间和展向进行平均, 即可获得时均流场.如无说明, 本文展示的均是无量纲化后的结果.

3 计算结果分析 3.1 统计结果

表 1给出了Strouhal数St及阻力系数C_d的对比, 可以看到, 本文计算结果St取0.139, 较多数的实验及DNS结果稍微偏大一些.相比St, 阻力系数的取值为2.17, 与实验及DNS的结果吻合较好.可以看出, 即便对这类统计参数, 要得到完全一致的结果也是非常困难的, 这从一个侧面反映了实现本算例精确预测的难度.

3.2 时均流场

图 2给出了时均流场中柱体附近的压力分布及流线图, 分别为文献[13]的DNS结果及本文LES结果, 尽管没有给出压力的量级, 但显然, 二者在分布上是非常相似的; 从时均流线看, LES与DNS预测的分离区形态及尺寸基本一致.

图 3展示了方柱4个表面上压力系数C_p分布的对比, 可以看到, LES与DNS曲线的分布规律完全一致, 但二者的数值存在一些差异, 这种差异应该来源于本文算例取M=0.3, 而DNS是针对不可压流动进行的.

图 4展示了x=-0.125和0.125位置处流向速度u, Reynolds应力分量$\overline{{u}'{u}'}$和$\overline{{u}'{v}'}$的对比, 除DNS结果外, 图中还给出了Lyn等^[1-2]、Minguez等^[11]的实验数据.可以看到, 对于流向速度u和$\overline{{u}'{u}'}$, 本文LES和DNS与实验数据均符合得较好; 对于$\overline{{u}'{v}'}$, LES曲线与DNS结果的差异较大, 但与Minguez等^[11]的激光Doppler测速(laser Doppler velocimetry, LDV)结果符合得较好.

图 5展示了柱体两侧壁面不同位置处的平均速度剖面, 其中实线为LES结果, 短点线为DNS结果.由流向速度u和横向速度v的结果都可以发现, LES与DNS的曲线吻合得很好.图 6展示了方柱尾迹区不同位置处(x=0~8) 的平均速度u及Reynolds应力$\overline{{u}'{u}'}$剖面.流向速度剖面图中, 靠近方柱时LES和DNS结果吻合较好, 朝向下游, 中心线附近的结果出现了一定的偏离.对于$\overline{{u}'{u}'}$, 除了x=1附近有一些偏差, 整个尾迹区的结果都吻合得较好.

图 7给出了尾迹中心线上压力系数C_p和Reynolds应力$\overline{{u}'{u}'}$的分布, 由图可知, LES的C_p计算结果略高于DNS, 而对于$\overline{{u}'{u}'}$的结果LES则低于DNS和实验数据.这可能是由于尾迹区LES的网格不够密, 导致流动结构耗散加快, 因此尾迹区速度恢复得更快, Reynolds应力也偏小.

3.3 瞬态流场

图 8给出了计算所得某一瞬态流场展向截面的涡量云图, 从图中可以看到, 涡量场存在从方柱左侧表面两个角点分离形成的自由剪切层, 剪切层在扰动作用下很快失稳, 产生展向涡结构.图 9给出了同一瞬态流场中的大尺度结构, 旋涡识别采用Q准则, 展向涡量着色, 容易看到, 角点两侧的展向涡继续快速失稳, 从二维演化出复杂的三维结构.图 10为方柱尾迹区的大尺度结构, 通过调小Q值, 从图中可以清楚地观察到Karman涡街的存在, 相邻展向大涡之间有着拉伸的准流向涡结构.

4 结论

采用高精度大涡模拟方法, 对来流Mach数0.3, Reynolds数22 000的孤立方柱绕流进行了计算, 从时均流场、2阶统计量、瞬态流场等角度分析了计算结果, 主要结论如下:

(1) 对于较高Reynolds数下的孤立方柱流动而言, 大涡模拟方法能够以相对很少的计算资源获得与实验和DNS吻合很好的计算结果, 验证了方法的可靠性;

(2) LES的瞬态流场展示了计算条件下方柱绕流分离转捩及尾迹区旋涡交替脱落形成Karman涡街的全过程, 为更细致的流动机理研究提供参考;

(3) LES的预测结果仍然存在尾迹中心流向速度恢复较快的问题, 其原因有待进一步的探讨.

致谢: 本文所有计算工作均在国家超级计算广州中心的天河二号系统上完成, 感谢广州中心提供的超算资源及相关技术支持.