1.   引　言

钝体绕流问题无论在理论上还是在工程实践中, 都有着重要的研究价值. 所谓亚临界区雷诺数钝体绕流, 即边界层为层流状态, 而尾流则为湍流状态的一类流动现象. 在理论上亚临界区雷诺数钝体绕流场研究中, 最具挑战性的难题当属其复杂三维时空瞬态发展的相干结构特性问题, 包括上游低强度湍流、自由剪切层中相干结构的产生与发展、高强度湍流溃灭, 以及随后产生的涡泄现象等^[1].

对圆柱绕流问题, 其亚临界区雷诺数的范围一般地定义为$400 \leqslant Re \leqslant 2.0 \times {10^5}$. 在这个雷诺数区圆柱绕流场中除了会出现大尺度涡泄这类相干结构外, 还会出现其他一些更为复杂的三维时空瞬态发展的相干结构现象, 包括剪切层小尺度Kelvin-Helmholtz不稳定性(K-H不稳定性)结构, 在圆柱尾部会出现两类不同长度的回流区结构现象, 以及在圆柱尾部近壁面区的平均流向速度剖面会出现V形和U形两个流动分支等^[2].

研究表明, 在圆柱绕流雷诺数位于亚临界区的情况, 当雷诺数Re从400增大到约1200时, 从圆柱表面分离的剪切层开始出现不稳定性现象^[3], 称为K-H不稳定性. 研究进一步表明, 圆柱绕流场中周期性涡泄(Karman涡街)现象发生于雷诺数Re~190时^[4], 但当雷诺数Re达到5000附近时, 圆柱绕流涡泄出现突然转变现象^[3], 其主要特征表现为泄涡结构开始变得不再具有周期性, 这种现象可以维持到雷诺数$Re = 2.0 \times {10^5}$.

所谓圆柱绕流场中的K-H不稳定性, 是指一条速度不连续的切变线上产生涡度集中而导致的流动不稳定性现象. Karman涡街属于一类频率相对较低(频率记为${f_{{\text{vs}}}}$)的大尺度相干结构, 而K-H不稳定性则属于一类频率相对较高(频率记为${f_{{\text{kh}}}}$)的小尺度相干结构, 其主要特征表现为宽频的信号特性, 且其峰值频率受雷诺数Re的影响而显著变化. 在亚临界区雷诺数不大于5000的情况, 这两种不稳定性模式通常可以共存, 且两者的频率近似满足^[3]${{{f_{{\text{kh}}}}} \mathord{\left/ {\vphantom {{{f_{{\text{kh}}}}} {{f_{{\text{vs}}}}}}} \right. } {{f_{{\text{vs}}}}}} = 0.023 R{e^{0.67}}$.

在圆柱绕流泄涡为周期性大尺度Karman涡街的情况, 利用RANS模型通常能够获取其较为准确的涡泄频率${f_{{\text{vs}}}}$的值, 以及Karman涡街的流态结构. 由于RANS模型只能提供圆柱绕流场的时均量信息, 而不能获得其三维空间瞬态发展的信息, 因此对圆柱绕流大规模流动分离及剪切层小尺度K-H不稳定性等这类复杂流动问题, RANS模型并不适用^[5].

直接数值模拟(DNS)^[6-12]和大涡模拟 (LES)^{[6–8,13-23]}以及部分平均N-S(PANS)^[2,24–26]等尺度解析模拟(SRS)方法则可以弥补RANS的这种缺陷, 并已成为研究这类复杂常钝体绕流问题的主要手段. 总体上, 通过基于DNS, LES及PANS等大量CFD数值模拟研究工作, 并结合相关的模型实验^[18,27]研究工作, 对亚临界区雷诺数下圆柱绕流场涉及的层流分离、层流-湍流转捩、周期性涡脱落及剪切层不稳定性等复杂流动现象的形成机理及其特征等问题, 有了较为深入的认识.

在雷诺数Re = 3900下的圆柱绕流是典型的亚临界区雷诺数流动, 在其流动中除了存在大尺度涡泄及剪切层小尺度K-H不稳定性这两类相干结构外, 无论在模型实验中还是在CFD数值模拟中, 都发现在其流动中还会出现两类特殊的流动结构现象. 其中, 第一个特殊流动现象为在其绕流场中会出现两种不同长度的回流区结构, 而第二个特殊流动现象是在圆柱尾部近壁面区的平均流向速度剖面会出现V和 U形两个流动分支结构.

Lourenco和Shih ^[27]的实验发现, 该雷诺数下圆柱绕流的回流区长度为${{{L_{\text{r}}}} \mathord{\left/ {\vphantom {{{L_{\text{r}}}} D}} \right. } D}$ = 1.18, 并且在${{{x_1}} /D} = 1.06$处的平均流向速度剖面呈V形. 然而, 在同样的雷诺下, Parnaudeau等^[18]的实验发现, 圆柱绕流的回流区长度为${{{L_{\text{r}}}} \mathord{\left/ {\vphantom {{{L_{\text{r}}}} D}} \right. } D}$ = 1.51, 且在${{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} D}} \right. } D}=1.06$处的平均流向速度剖面呈U形. 同时, 两个实验所测得的不同站位处平均流向和横向速度及各向同性和异性雷诺应力剖面特性等均不相同. 为此, 众多学者采用CFD数值模拟方法对此问题进行了长时间的研究.

Tremblay^[8]采用DNS和LES, Breuer^[15]及Wong和Png^[19]采用LES, Pereira等^[2]及Luo等^[24]采用PANS, 均复现了Lourenco和Shih^[27]的实验结果. Parnaudeau等^[18]及Franke等^[17]采用LES, Song等^[11]及Ooi等^[12]采用DNS, 则均复现了Parnaudeau等^[18]的实验结果. Tremblay^[8]认为Lourenco和Shih^[27]的实验结果有误, Afgan等^[20]认为Lourenco和Shih^[27]的统计周期不够进而导致结果未收敛, Kravchenko等^[16]认为Lourenco和Shih^[27]的实验由于受到了外部干扰而导致剪切层过早转捩.

从目前的文献资料看, 许多学者认为Parnaudeau等^[18]的实验结果更具有可信度, 因此大部分学者都采用Parnaudeau等^[18]的实验结果作为CFD数值模拟的基准参考数据, 然而这种说法很难令人信服. 考虑到从实验、DNS、LES及PANS等诸多方面均可获得V形和U形结果, 出现这两种结果的背后必然有其更深层次的理论和物理方面的机制, 如圆柱体长径比、展向网格分辨率及湍流模型对绕流场湍动能的解析能力等.

Ma等^[7]采用DNS及LES, 对圆柱体展向长度的影响进行了研究, 发现当展向长度为2πD时产生V形结果, 而当展向长度为πD时产生U形结果, 并将产生两种结果的原因归于展向长度设置问题. Dong等^[9]采用DNS得出了类似结论. 然而, Kravchenko等^[16]基于LES的研究发现, 当展向网格分辨率为πD/8时, 产生V形结果, 而当展向网格分辨率为πD/48时, 产生U形结果, 但在保持足够密的展现网格分辨率时, 如果将展向长度从πD增大至2πD时并不会产生明显区别.

Xia等^[6]采用展向很疏的网格分辨率开展DNS研究(在展向仅布置四层网格), 却惊奇地复现了V形结果. Kim^[13]基于LES研究获得了与Kravchenko等^[16]一致的结论, 即展向网格分辨率是导致U形和V形两个分支结构的主要原因. 然而, Breuer^[15]采用LES的研究发现, 在设置展向πD/32及πD/64两种网格分辨率的情况下均出现V形结果.

对PANS方法, 其对钝体绕流场湍动能的解析能力通过一个滤波参数${f_{\text{k}}}\;( = {{{k_{\text{u}}}} \mathord{\left/ {\vphantom {{{k_{\text{u}}}} k}} \right. } k})$进行控制. 其中, $k$为总的湍动能, ${k_{\text{u}}}$为不可解湍动能. Pereira等^[2]及Kořínek等^[26]采用PANS对V形和U形问题进行了研究. 研究发现, 当${f_{\text{k}}}$大于某个值${f_{{\text{k0}}}}$时, 会产生V形结果, 而当${f_{\text{k}}}$小于某个值${f_{{\text{k0}}}}$时, 会产生U形结果, 他们将此现象归因于湍流模型对湍动能的解析能, 即具有高湍动能解析能力的湍流模型产生U形结果, 而具有低湍动能解析能力的湍流模型产生V形结果.

为准确获取亚临界区雷诺数圆柱绕流中剪切层小尺度K-H不稳定性结构高频成分的频率${f_{{\text{kh}}}}$, PANS模型中的滤波参数一般地需要满足${f_{\text{k}}} \leqslant 0.25$, 即PANS对湍动能的解析能力至少需要达到75%^[2]. 这意味着, 对PANS而言, 为准确获取亚临界区雷诺数下圆柱绕流场中剪切层小尺度K-H不稳定性结构高频成分的频率${f_{{\text{kh}}}}$, 其所需网格量应当与LES的网格量相当.

DNS需要解析边界层内所有尺度的湍流, 对网格分辨率的要求特别高, 需要的网格量特别巨大, 不适合高雷诺数钝体绕流CFD计算. LES直接解析大尺度的湍流, 而小尺度湍流则用亚格子应力模化, 虽然与DNS相比网格量要少很多, 但对高雷诺数钝体绕流的工程化应用仍很遥远^[28]. 由此可见, 对PANS而言, 其工程化应用前景亦是如此.

在诸多湍流模型中, 兼顾计算精度与资源消耗的混合RANS/LES(HRL)模型已成为当今CFD领域研究与应用的热点, 包括脱体涡模拟(DES)、延迟脱体涡模拟(DDES)及增强版脱体涡模拟(IDDES)等^[29]. 下文将这三类模型统称为DES类模型.

D'Alessandro等^[30]基于DES, 对不同网格分辨率及湍流模型的能力进行了研究与评估, 认为V形及U形两种结果与网格分辨率及相应湍流模型性能关系密切. 研究表明, 标准SA-DES模型在不同加密程度的网格下仅能预测V形结果, 而SA-IDDES模型在很密网格下可预测U形结果, 在较疏网格下则可预测V形结果.

综上所述, 目前的研究仍存在诸多未解问题, 主要为对V形及U形两个平均流向速度剖面分支结构产生的机理尚不能形成统一的认识. 特别地, 为何Ma等^[7]的结果与其他相关文献的结果矛盾? 第二, 为何Breuer^[15]所用展向网格分辨率与Kravchenko等^[16]相同, 但获得的结果并不一致? 第三, SA-DES与SA-IDDES模型所采用的RANS和LES模型一样, 其区别主要在于RANS与LES之间的转换方式不同, 为何网格分辨率对其结果却有如此大的影响^[30]?

另一方面, 根据目前所能查阅到的文献资料, 利用混合RANS/LES模型来研究亚临界区雷诺数圆柱绕流的模型主要为DES及DDES^[24,30]. 虽然这两类模型能够较为准确地获得与实验结果一致的一阶统计量(压力系数、流向及横向速度剖面分布等)的计算结果, 但对二阶统计量(各向同性及异性雷诺应力剖面分布等)的计算结果与实验结果仍有较大差异. 同时, 由于这类模型对边界层中小尺度流动结构的解析能力有限, 因此难以准确获取圆柱绕流剪切层小尺度K-H不稳定性结构的精细信息, 特别是其频谱结构及频率${f_{{\text{kh}}}}$的准确值.

有鉴于此, 本文发展了一种壁面模化混合RANS/LES模型(WM-HRL), 该方法也属于一类HRL方法. WM-HRL模型与传统DES类模型的主要不同之处在于, 可实现自剪切层小尺度K-H不稳定性结构发生区域即做具有至少80%湍动能解析能力的完全LES计算, 不仅可有效地减少计算网格的数量, 而且还可以有效解析剪切层小尺度K-H不稳定性结构特征, 并准确地捕捉其频谱结构及特征频率等信息.

在此基础上, 本文以亚临界区雷诺数Re = 3900下的圆柱绕流问题为对象, 对该WM-HRL模型的能力进行系列数值模拟和评估研究, 包括对圆柱绕流场中剪切层小尺度K-H不稳定性和两类不同长度回流区精细结构解析与捕捉能力的评估, 以及对圆柱尾部${{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} D}} \right. } D}$ = 1.06处平均流向速度剖面出现V形和U形两个流动分支结构形成机理的进一步认识等.

2.   理论模型

考虑不可压缩流体的钝体绕流问题. 设流 体密度为${\rho _0}$, 运动黏度系数为$\nu$. 建立直角坐标系为$o{\text{-}}{x_1}{x_2}{x_3}$, 其中$o{x_3}$轴垂直向上为正, ${\boldsymbol{u}} = ({u_1}, {u_2}, {u_3})$为流体运动的速度矢量. 对各类混合RANS/LES模型, 虽然其在钝体近壁面区取采用RANS进行计算, 而在远离钝体近壁面的区域采用LES进行计算, 但两者均采用如下RANS统一框架下的控制方程对流场进行计算:

$$\begin{split} &\frac{\partial {\overline{u}}_{k}}{\partial {x}_{k}}=0,\\ \frac{\partial {\overline{u}}_{i}}{\partial t}+{\overline{u}}_{k}\frac{\partial {\overline{u}}_{i}}{\partial {x}_{k}} =\;&-\frac{1}{{\rho }_{0}}\frac{\partial \overline{p}}{\partial {x}_{i}}+\frac{\partial }{\partial {x}_{k}}\left(\nu \frac{\partial {\overline{u}}_{i}}{\partial {x}_{k}}\right)+\frac{\partial {\tau }_{ki}}{\partial {x}_{k}}, \\ \end{split}$$

(1)

其中, 顶标“—”表示雷诺时均. ${\tau _{ki}}$为Reynolds应力, 采用如下Boussinesq近似进行计算:

$${\tau _{ki}} = 2{\nu _{\text{t}}}{\bar S_{ki}} - \frac{2}{3}k{\delta _{ki}},$$

(2)

其中, $k$为湍动能, $\omega$为比耗率, $\varepsilon = {\beta ^*}k\omega$为耗散率, ${\bar S_{ki}}$为形变率张量, ${v_{\text{t}}}$为涡黏系数.

为封闭RANS方程(1), 需要引入相应的湍流模型, 本文采用SST k-ω模型. 基于该湍流模型的混合RANS/LES模型, 可通过修改其k方程中的色散项而建立, 具体如下^[31]:

$$\frac{{\partial k}}{{\partial t}} + {\bar u_i}\frac{{\partial k}}{{\partial {x_i}}} = {P_{\text{k}}} - \frac{{{k^{{{3} \mathord{\left/ {\vphantom {{3} {2}}} \right. } {2}}}}}}{{{{\tilde l}_{{\text{hyb}}}}}} + \frac{\partial }{{\partial {x_j}}}\left[ {\left( {\nu +{\sigma _k}{\nu _{\text{t}}}} \right)\frac{{\partial k}}{{\partial {x_j}}}} \right],$$

(3)

$$\begin{split}&\frac{{\partial \omega }}{{\partial t}} + {\bar u_i}\frac{{\partial \omega }}{{\partial {x_i}}} = \frac{{\tilde \gamma }}{{{\nu _{\text{t}}}}}{P_{\text{k}}} - \beta {\omega ^{2}}+\frac{\partial }{{\partial {x_j}}}\left[ {\left( {\nu +{\sigma _\omega }{\nu _{\text{t}}}} \right)\frac{{\partial \omega }}{{\partial {x_j}}}} \right]\\ & ~~~~~~~~~~~~~~~ ~~~~~~~~~~ +2\left( {{1} - {F_{1}}} \right)\frac{{{\sigma _{\omega {2}}}}}{\omega }\frac{{\partial k}}{{\partial {x_j}}}\frac{{\partial \omega }}{{\partial {x_j}}},\\[-1pt] \end{split}$$

(4)

其中, ${\tilde l_{{\text{hyb}}}}$为混合长度, ${P_{\text{k}}}$为产生项, 相关参数见文献[32]. ${v_{\text{t}}}$为湍流涡黏系数, 定义为

$${\nu _{\text{t}}} = \frac{{{\alpha _1}k}}{{\max \left( {{\alpha _1}\omega ,\left| {\bar S} \right|{F_2}} \right)}},$$

(5)

其中, ${\alpha _1}$为模型系数, 取值为0.31; F₂为混合函数, $\left| {\bar S} \right|$为形变率张量的模.

为保证迭代求解的稳定性, Menter等^[32]加入了湍流黏度的限制性条件, 对湍流动能生成项${P_{\text{k}}}$进行了如下的修正:

$${P_{\text{k}}} = \min \left\{ {{v_{\text{t}}}\left| {\bar S} \right|{}^2,10{\beta ^*}k\omega } \right\},$$

(6)

其中, ${\beta ^{\text{*}}}$为SST的模型常数, 取值为0.09.

对DES模型^[29], ${\tilde l_{{\text{hyb}}}}$定义为

$${\tilde l_{{\text{hyb}}}} = \min \left( {{l_{{\text{RANS}}}},{l_{{\text{LES}}}}} \right),$$

(7)

其中, ${l_{{\text{RANS}}}}$为RANS尺度, ${l_{{\text{LES}}}}$为LES尺度, 可分别表示为

$${l}_{\text{RANS}}={k}^{1/2}/\left({\beta }^{\text{*}}\omega \right),\;{l}_{\text{LES}}={C}_{\text{DES}}\varDelta ,$$

(8)

其中, C_DES为模型常数, $\varDelta$为LES滤波宽度.

在(8)式中, ${l_{{\text{RANS}}}}$和${l_{{\text{LES}}}}$与当地网格中湍流特征尺度相关. 当${l_{{\text{RANS}}}} < {l_{{\text{LES}}}}$时, ${\tilde l_{{\text{hyb}}}} = {l_{{\text{RANS}}}}$, 此时DES为RANS模式; 当${l_{{\text{RANS}}}} > {l_{{\text{LES}}}}$时, ${\tilde l_{{\text{hyb}}}} = {l_{{\text{LES}}}}$, 此时DES为LES模式. 由此可见, DES模型没有明确的RANS/LES混合转换界面. 该模型从RANS到LES的转换完全由RANS和LES尺度的相对大小决定, 或者说主要由网格尺度的空间分布决定. 因此, 在使用该模型时, 通常要求流向和展向网格尺度不能同时小于边界层的厚度.

然而, 当边界层内流向网格和展向网格加密至某个中间情况时, 即当网格尺度小于RANS尺度而又远大于LES计算壁湍流所需的网格尺度时, 此时DES中的LES尺度会提前进入边界层, 导致边界层内模化的雷诺应力不足(MSD), 缺失的湍流脉动又不足以被解析出来, 当流向遇到一定的逆压梯度时, 则产生非物理的流动分离现象, 即网格诱导非物理分离(GIS), 甚至发生对数律层不匹配(LLM)现象^[33].

有鉴于此, Spalart等^[33]提出可在混合长度尺度${\tilde l_{{\text{hyb}}}}$中引入边界层的识别函数来延迟RANS到LES的转换, 从而防止LES提前进入边界层, 进而得到DDES方法, 具体如下:

$${\tilde{l}}_{\text{hyb}}={l}_{\text{RANS}}-{f}_{\text{d}}{\mathrm{max}}\left(0,{l}_{\text{RANS}}-{l}_{\text{LES}}\right),$$

(9)

其中, ${f_{\text{d}}}$为一个经验性混合函数, 具体形式如下:

$$\begin{split}&{f}_{\text{d}}=1-{\mathrm{tanh}}\left[{\left(8{r}_{\text{d}}\right)}^{3}\right],\\ &{r}_{\text{d}}= \dfrac{\nu +{v}_{\text{t}}}{{\kappa }^{2}{d}_{\text{w}}^{2}\sqrt{\left(\partial {\overline{u}}_{k}/\partial {x}_{i}\right)\left(\partial {\overline{u}}_{k}/\partial {x}_{i}\right)}},\end{split}$$

(10)

其中, ${d_{\text{w}}}$为网格计算点到壁面的距离, $\kappa$为冯卡门常数, 取值为0.41.

函数${f_{\text{d}}}$的分布特点为, 在近壁层的某个范围内(与网格及当地流场有关)等于0, 即在此区域内混合长度${\tilde l_{{\text{hyb}}}} = {l_{{\text{RANS}}}}$, 此时DDES模型还原为RANS模式. 在此区域外的计算区域, 混合长度决定于两个尺度${l_{{\text{RANS}}}}$和${l_{{\text{LES}}}}$的相对大小, 此时DDES与DES的表现是一样的.

DDES模型虽然解决了原始DES模型存在的MSD及GIS等缺陷, 但其仍继承了DES存在的其他问题, 如LLM缺陷. 为此, Shur等^[34]提出了增强版的DDES模型(IDDES), 但在原始IDDES定义的混合长度中, 含有一个所谓的提升函数f_e, 这是一个纯人工构造的函数. Gritskevich等^[31]指出, 该函数的作用仅为增加涡黏系数, 但这种人为增加涡黏系数的作用是否合理有待商榷, 因此建议采用如下的简化版本:

$${\tilde{l}}_{\text{hyb}}={\tilde{f}}_{\text{d}}\cdot {l}_{\text{RANS}}+\left(1-{\tilde{f}}_{\text{d}}\right){l}_{\text{LES}},$$

(11)

其中

$${\tilde{f}}_{\text{d}}={\mathrm{max}}\left\{{f}_{\text{B}},{f}_{\text{dt}}\right\},$$

(12)

$$\begin{split} &{f}_{\text{dt}}=1-{f}_{\text{d}},{f}_{\text{B}}={\mathrm{min}}\left(2{\mathrm{exp}}\left(-9{\alpha }^{2}\right),1.0\right),\\ &\alpha =0.25-{d}_{\text{w}}/{h}_{{\mathrm{max}}},\end{split}$$

(13)

${h_{\max }}$为计算单元的最大网格步长.

对如上构造的IDDES模型, 当${\tilde f_{\text{d}}} = 1$时, ${\tilde l_{{\text{hyb}}}} = {l_{{\text{RANS}}}}$, 此时IDDES为完全RANS模式. 当${\tilde f_{\text{d}}} = 0$时, $\tilde l = {l_{{\text{LES}}}}$, 此时IDDES为完全LES模式. 当$0 < {\tilde f_{\text{d}}} < 1$时, IDDES执行RANS/LES混合解析模式. IDDES模型通过修改混合长度的定义和引入具有壁面模化能力的保护层函数${f_{\text{B}}}$, 相当大程度上缓解了MSD, GIS及LLM的问题.

对DES, DDES和IDDES模型, 当其进入LES后, 在局部平衡流条件下, SST模型k方程中的生成项与色散项相等^[35], 即

$${v_{\text{t}}}{\left| {\bar S} \right|^2} = {{{k^{{3 / 2}}}}/ {({C_{{\text{DES}}}}\varDelta }}).$$

(14)

此外, 由文献[35]可得

$${v_{\text{t}}}{\left| {\bar S} \right|^2} = 0.3k\left| {\bar S} \right|.$$

(15)

由(14)式和(15)式可得

$${v_{\text{t}}}{\left| {\bar S} \right|^2} = 0.3k\left| {\bar S} \right| = {{{k^{{3 / 2}}}} / {({C_{{\text{DES}}}}\varDelta }}),$$

(16)

由(16)式可得

$$k = \frac{{{v_{\text{t}}}\left| {\bar S} \right|}}{{0.3}} = \frac{{{v_{\text{t}}}\left| {\bar S} \right|}}{{{{\left( {{\beta ^{\text{*}}}} \right)}^{{1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}}}}.$$

(17)

由(16)式和(17)式, 可得:

$${v_{\text{t}}} = {\left( {{{\left( {{\beta ^{\text{*}}}} \right)}^{{3/ 4}}}{C_{{\text{DES}}}}\varDelta } \right)^2}\left| {\bar S} \right| = {\left( {{C_{\text{S}}}\varDelta } \right)^2}\left| {\bar S} \right|,$$

(18)

其中

$${C_{\text{S}}}={\left( {{\beta ^{\text{*}}}} \right)^{{3 \mathord{\left/ {\vphantom {3 4}} \right. } 4}}}{C_{{\text{DES}}}}.$$

(19)

(19)式正好与经典Smagorinsky亚格子涡黏模型一致. 因此, DES, DDES和IDDES模型都是利用SST模型k方程中产生项与色散项平衡这种极限情况下来间接等效LES模式而实现的.

在(8)式中, 滤波尺寸$\varDelta$控制LES能否在Kolmogorov能量谱尺度下解析尽可能多含能尺度的湍流场. 在经典DES, DDES及IDDES模型(DES类模型)中, $\varDelta$一般取为最大网格尺寸${\varDelta _{\max }}$. 根据这3类DES类模型中“局部平衡流”的假定, 即在边界层外的区域, 要实现DES类模型从RANS模式转换为LES模式, 其中SST湍流模型k方程中生成项与色散项需要达到平衡, 此时DES类模型相当于经典的Smagorinsky型亚格子涡黏模型.

Breuer等^[36]认为, 在DES类模型中采用最大网格尺寸并不合适, 进而建议采用如下的体积立方根尺度:

$${\varDelta _{{\text{vol}}}} = {V^{{1 \mathord{\left/ {\vphantom {1 3}} \right. } 3}}},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} V = {\varDelta _1}{\varDelta _2}{\varDelta _3},$$

(20)

其中, ${\varDelta _1},\; {\varDelta _2}$和${\varDelta _3}$分别为3个坐标方向的网格尺寸.

Reddy等^[37]针对DDES模型建议了一个新的网格混合形式如下:

$$\varDelta = {f_{\text{d}}}{\varDelta _{{\text{vol}}}} + \left( {1 - {f_{\text{d}}}} \right){\varDelta _{\max }}.$$

(21)

Shur等^[34]则建议采用如下的定义:

$$\varDelta = \min \left\{ {\max \left[ {{C_{\text{w}}}{d_{\text{w}}},{C_{\text{w}}}{h_{\max }},{h_{{\text{wn}}}}} \right],{h_{\max }}} \right\},$$

(22)

其中, ${C_{\text{w}}} = 0.15$, ${h_{{\text{wn}}}}$为壁面法向网格步长.

对两方程的SST湍流模型, (8)式中C_DES的取值可采用如下加权形式^[31]:

$${C_{{\text{DES}}}} = {F_1}{C_{{\text{DES,in}}}} + \left( {1 - {F_1}} \right){C_{{\text{DES,out}}}}.$$

(23)

在(23)式中, ${C_{{\text{DES, in}}}}$为IDDES内层RANS分支的系数, ${C_{{\text{DES, out}}}}$为IDDES类模型外层LES分支的系数, 一般地可按下式取值:

$${C}_{\text{DES,in}}=0.78,\;\;{C}_{\text{DES,out}}={\mathrm{0.61.}}$$

(24)

对IDDES模型, 在边界层内一般为完全RANS模式, 而其完全LES模式一般发生在边界层外^[38]. 由此可见, IDDES除了可以避免MSD及GIS的问题外, 也可以避免LLM的问题. 本文研究的亚临界雷诺数圆柱绕流问题, 其剪切层小尺度K-H不稳定性发生在对数律层区内, 由于RANS模式难以准确地捕捉到这类小尺度K-H不稳定性结构的三维时空瞬态发展流动的精细结构, 因此IDDES同样难以高置信度地解析这类非定常、非平衡流动现象的精细结构特性.

克服IDDES模型缺陷的一种有效途径是使其具有壁面模化大涡模拟的能力. 有鉴于此, 本文构造一种新的混合函数${f_{{\text{nd}}}}$, 使新所构造的混合RANS/LES模型, 除了具有延迟脱体涡模拟(DDES)的能力外, 同时具有壁面模化大涡模拟(WMLES)的能力, 将其称为WM-HRL模型.

为此, 首先引入湍动能解析度指标概念, 即${r_{\text{k}}} = {{{k_{\text{u}}}} \mathord{\left/ {\vphantom {{{k_{\text{u}}}} k}} \right. } k}$, 其中${k_{\text{u}}}$为未解湍动能. 当${r_{\text{k}}} = 1$时, 可得${k_{\text{u}}} = k$, 即湍动能被完全模化, 此时WM-HRL为完全RANS模式. 当${r_{\text{k}}} = 0$时, 即${k_{\text{u}}} = 0$, 即湍动能被完全解析, 此时WN-HRL为完全DNS模式. 对LES来说, 为准确地捕捉剪切层小尺度K-H不稳定性结构这类特殊流动现象, 其对湍动能的解析能力至少需要达到75%^[2], 即${r_{\text{k}}} \leqslant 0.25$.

定义${k_{\text{c}}}$为截断波数, 它表示LES的滤波宽度, 可由当地网格尺寸${\varDelta ^{\text{*}}}$确定如下:

$${k_{\text{c}}} = {{\text{π }}/{{\varDelta ^*}}}.$$

(25)

根据Kolmogorov的${{ - 5} \mathord{\left/ {\vphantom {{ - 5} 3}} \right. } 3}$谱定律, 当${k_{\text{c}}}$位于惯性亚区时, ${k_{\text{u}}}$可由下式进行计算:

$$\begin{split} {k_{\text{u}}} =\;& \int_{{k_{\text{c}}}}^{ + \infty } {E(k){\text{d}}k} = \int_{{k_{\text{c}}}}^{ + \infty } {{C_{\text{k}}}{\varepsilon ^{{2 / 3}}}{k^{{{ - 5} / 3}}}} {\text{d}}k \\ =\;& \frac{3}{2}{C_{\text{k}}}{\varepsilon ^{{2 / 3}}}k_{\text{c}}^{{{{{-2}}} / {3}}}, \end{split}$$

(26)

其中, ${C_{\text{k}}}~( \approx 1.5)$为Kolmogorov常数. 由此可得:

$${r_{\text{k}}} = \frac{{{k_{\text{u}}}}}{k}=\frac{{\dfrac{3}{2}{C_{\text{k}}}{\varepsilon ^{{2 / 3}}}k_{\text{c}}^{{{{{ - 2}}} /{3}}}}}{k} = {\left( {\frac{{{L_{{\text{sgs}}}}}}{{{L_{{\text{tur}}}}}}} \right)^{{2 / 3}}},$$

(27)

其中, ${L_{{\text{sgs}}}}$为亚格子滤波尺度, ${L_{{\text{tur}}}}$为湍流含能尺度, 它们可分别表示如下:

$${L}_{\text{sgs}}=\left[{\left(3{C}_{\text{k}}/2\right)}^{3/2}/\text{π}\right]{\varDelta }^{*}\approx {\varDelta }^{*},{L}_{\text{tur}}={k}^{3/2}/\varepsilon .$$

(28)

由(27)式和(28)式可得:

$${r_{\text{k}}} = {\left( {{{{\varDelta ^*}} /{{L_{{\text{tur}}}}}}} \right)^{{2 /3}}}.$$

(29)

对当地网格尺度${\varDelta ^{\text{*}}}$, 一种合理的选择^[39]为${\varDelta ^{\text{*}}}={h_{\max }}$. 此时, ${r_{\text{k}}}$可以改写为

$${r_{\text{k}}} = {\left( {{{{h_{\max }}} \mathord{\left/ {\vphantom {{{h_{\max }}} {{L_{{\text{tur}}}}}}} \right. } {{L_{{\text{tur}}}}}}} \right)^{{2 \mathord{\left/ {\vphantom {2 3}} \right. } 3}}}.$$

(30)

根据Nyquist-Shannon样本定理, 为了能够准确解析相关尺度的湍流结构, 湍流含能尺度${L_{{\text{tur}}}}$需要满足如下条件^[39]:

$${L_{{\text{tur}}}} > {N_{\text{r}}}{h_{\max }},$$

(31)

其中, ${N_{\text{r}}}$为某个长度尺度结构能够被解析的网格单元数量, 它与单位波长的波能够被重构所需的样本点数量一致.

对RANS模式, 由于其不能对湍流结构进行直接解析, 因此${L_{{\text{tur}}}} \leqslant {h_{\max }}$. 由(30)式可得, 此时${r_{\text{k}}} \geqslant 1$. 对LES模式, Larsson等^[39]指出, 湍流大尺度脉动场信息能够被准确捕捉的条件是${N_{\text{r}}}$至少为2, 即${L_{{\text{tur}}}} \geqslant 2{h_{\max }}$, 将其代入(30)式, 可得

$${r_{\text{k}}} \leqslant {\bar r_{{\text{k}}1}},$$

(32)

其中, ${\bar r_{{\text{k1}}}} = {\left( {{1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}} \right)^{{2 \mathord{\left/ {\vphantom {2 3}} \right. } 3}}}=0.63$. (32)式表明, 当${N_{\text{r}}} = 2$时, 湍动能解析度指标${\bar r_{{\text{k}}1}} = 0.63$. 这意味着, 一旦LES模式被激活, 则将至少有37%的湍流大尺度脉动信息能够被直接解析. 这同时意味着, 当${r_{\text{k}}} > {\bar r_{{\text{k1}}}}$时, LES将不能被真正激活.

有鉴于此, 设${r_{{\text{k1}}}} \leqslant {\bar r_{{\text{k1}}}}$. 本文将当地网格满足条件${r_{\text{k}}} > {r_{{\text{k}}1}}$的计算区域定义为WM-HRL模型的完全RANS区域, 即

$${D}_{\text{RANS}}=\left\{P|{r}_{\text{k|}P} > {r}_{\text{k1}},P\in {D}_{\text{copm}}\right\},$$

(33)

其中, $P$为计算区域${D_{{\text{comp}}}}$中的网格单元, ${r_{{\text{k|}}P}}$为网格单元$P$上的湍动能解析度指标值.

对WM-HRL模型来讲, 其第2个关键是如何合理地确定RANS/LES混合转化区域${D_{{\text{hyb}}}}$. 为此, 设${r_{{\text{k2}}}}$为WM-HRL模型进入完全LES模式时人们期望的大尺度湍流之解析度指标值. 据此, 可以将WM-HRL的RANS/LES混合转换区域${D_{{\text{hyb}}}}$定义为

$${D}_{\text{hyb}}=\left\{P|{r}_{\text{k2}} < {r}_{\text{k|}P} < {r}_{\text{k}1},P\in {D}_{\text{copm}}\right\}.$$

(34)

对PANS模型, ${r_{\text{k}}}$也是衡量其对钝体绕流场大尺度湍流解析能力的关键性控制参数. Pereira等^[2]对Re = 3900下圆柱绕流问题进行了系列数值模拟, 研究发现, 当${r_{\text{k}}} > 0.5$时, PANS模型尚不足以充分解析绕流场中大尺度涡结构的信息. 同时, 研究表明: 只有当${r_{\text{k}}} \leqslant 0.5$时, LES才具有较好地解析大尺度涡结构信息的能力. 因此, 把${r_{{\text{k2}}}}$取为小于0.5的某个值将是一种合理的选择.

当${r_{{\text{k2}}}}=0.5$时, 由(30)式可知, ${L_{{\text{tur}}}} \geqslant 2.83{h_{\max }}$. 再结合(31)式, 可取${N_{\text{r}}} = 3$. 此时, ${L_{{\text{tur}}}} \geqslant 3{h_{\max }}$, 将其代入(30)式可得${r_{\text{k}}} \leqslant 0.48$. 由此可见, 可将${\bar r_{{\text{k2}}}}$的取值修正为${\bar r_{{\text{k}}2}}=0.48$. 这意味着, 当WM-HRL模型进入完全LES模式后, 其在区域${D_{{\text{hyb}}}}$中对大尺度涡结构的解析能力至少可达52%.

至此, 对本文将构造的一种新的WM-HRL模型中如何合理地确定两个关键区域${D_{{\text{RANS}}}}$和${D_{{\text{hyb}}}}$进行阐述, 分别给出确定WM-HRL模型进行完全RANS模式的区域${D_{{\text{RANS}}}}$之${r_{{\text{k1}}}}$准则, 以及确定WM-HRL模型进行RANS/LES混合转换模式的区域${D_{{\text{hyb}}}}$之${r_{{\text{k}}2}}$准则. 其中, ${r_{{\text{k1}}}} \leqslant {\bar r_{{\text{k1}}}}$且${r_{{\text{k2}}}} \leqslant {\bar r_{{\text{k2}}}}$. 后文将其分别称为${r_{{\text{k1-Q}}}}$准则和${r_{{\text{k2-Q}}}}$准则.

对如上所述的${r_{{\text{k1-Q}}}}$和${r_{{\text{k2-Q}}}}$准则, 需要利用(30)式进行计算, 其中RANS含能尺度${L_{{\text{tur}}}}$在进行CFD计算之前属于未知量. 因此, 这两个准则无法用于在进行CFD计算网格设置时来具体确定两个关键区域${D_{{\text{RANS}}}}$和${D_{{\text{hyb}}}}$的边界位置.

为此, 下面继续构造${r_{\text{k}}}$的一种仅依赖于当地网格尺寸的计算公式. Pope^[40]指出, 在对数律层区, ${L_{{\text{tur}}}}$与${d_{\text{w}}}$成正比关系, 即${L_{{\text{tur}}}} = {C_{\text{w}}}{d_{\text{w}}}$. 在高雷诺数的情况, Han等^[41]指出, ${C_{\text{w}}}$可近似取为${C_{\text{w}}} \approx 2.5$. 由此可见, (30)式可以改写为:

$${r_{\text{k}}} = {\left( {0.4\frac{{{h_{\max }}}}{{{d_{\text{w}}}}}} \right)^{{2 \mathord{\left/ {\vphantom {2 3}} \right. } 3}}}.$$

(35)

对(35)式定义的${r_{\text{k}}}$, 其值有可能会出现大于1的情况. 为避免这种情况, 将 (35) 式修改为如下形式:

$${r_{\text{k}}} = \min \left[ {1.0,{{\left( {0.4\frac{{{h_{\max }}}}{{{d_{\text{w}}}}}} \right)}^{{2 \mathord{\left/ {\vphantom {2 3}} \right. } 3}}}} \right] ,$$

(36)

在(36)式的形式下, 所定义的湍动能解析度指标参数${r_{\text{k}}}$将始终不会超过1.0, 即${r_{\text{k}}} \leqslant 1.0$. 根据(36)式可得如下结论.

首先, 当${r_{\text{k}}} \geqslant {r_{{\text{k}}1}}$时, ${{{d_{\text{w}}}} \mathord{\left/ {\vphantom {{{d_{\text{w}}}} {{h_{\max }}}}} \right. } {{h_{\max }}}} \leqslant 0.4{\left( {{r_{{\text{k}}1}}} \right)^{ - {3 \mathord{\left/ {\vphantom {3 2}} \right. } 2}}}$. 在${r_{\text{k}}} \geqslant {\bar r_{{\text{k}}1}}$时, 可得${{{d_{\text{w}}}} \mathord{\left/ {\vphantom {{{d_{\text{w}}}} {{h_{\max }}}}} \right. } {{h_{\max }}}} \leqslant 0.8$, 此时LES将不能被真正激活. 一般地, 可将WM-HRL模型为完全RANS的区域定义为

$$\begin{split}{D}_{\text{RANS}}=\;&\{P|\left({d}_{\text{w}}/{h}_{{\mathrm{max}}}\right)|{}_{P}\leqslant 0.4{\left({r}_{\text{k}1}\right)}^{-3/2},\\ &P\in {D}_{\text{copm}}\}. \end{split}$$

(37)

其次, 当${r_{{\text{k}}2}} < {r_{\text{k}}} < {r_{{\text{k1}}}}$时, $0.4{{( {{r_{{\text{k1}}}}} )}^{ - {3 /2}}} < {d_{\text{w}}}/ {{h_{\max }}} < 0.4{\left( {{r_{{\text{k2}}}}} \right)^{ - {3 /2}}}$. 在${r_{{\text{k2}}}} \leqslant {\bar r_{{\text{k2}}}}$时, 可得${{{d_{\text{w}}}} \mathord{\left/ {\vphantom {{{d_{\text{w}}}} {{h_{\max }}}}} \right. } {{h_{\max }}}} \geqslant 1.2$, 此时LES正好被激活. 一般地, 可将WM-HRL模型为RANS/LES混合转换的区域定义为

$$\begin{split}{D}_{\text{hyb}}=\;&\{P|0.4{\left({r}_{\text{k1}}\right)}^{-3/2} < {d}_{\text{w}}/{h}_{{\mathrm{max}}} < 0.4{\left({r}_{\text{k}2}\right)}^{-3/2},\\ &P\in {D}_{\text{copm}}\}.\\[-1pt] \end{split}$$

(38)

将由(37)式确定的WM-HRL模型为完全RANS的区域称为$n{r_{{\text{k1-Q}}}}$准则, 而将由(38)式确定之WM-HRL模型为RANS/LES混合转换的区域称为$n{r_{{\text{k2-Q}}}}$准则. 由此可知, 与前面所述之原${r_{{\text{k1-Q}}}}$和${r_{{\text{k2-Q}}}}$准则相比, 新的$n{r_{{\text{k1-Q}}}}$和$n{r_{{\text{k2-Q}}}}$准则仅与当地网格尺寸相关, 因此可以很容易地根据这两个准则来进行WM-HRL模型中两个关键区域${D_{{\text{RANS}}}}$和${D_{{\text{hyb}}}}$的确定及其相应的网格设置.

根据如上所述两个$n{r_{{\text{k1-Q}}}}$和$n{r_{{\text{k2-Q}}}}$准则, 可构造一种新的混合函数${f_{{\text{nd}}}}$如下:

$${f_{{\text{nd}}}} = {1} - {\text{sign}}\left[ {{r_{{\text{k1}}}} - {m_{\text{k}}}} \right]{f_{\text{s}}},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {m_{\text{k}}} = {\text{min}}\left[ {{r_{{\text{k1}}}},{r_{\text{k}}}} \right],$$

(39)

其中, ${f_{\text{s}}}$为待定函数. 同时, 将混合长度修改为

$${\tilde l_{{\text{hyb}}}} = {f_{{\text{nd}}}}{\kern 1pt} {\kern 1pt} {l_{{\text{RANS}}}} + ({1} - {f_{{\text{nd}}}}){\kern 1pt} {\kern 1pt} {l_{{\text{LES}}}}.$$

(40)

由(39)式可知, 在当地网格满足$n{r_{{\text{k1-Q}}}}$准则的情况, 即${r_{\text{k}}} \geqslant {r_{{\text{k}}1}}$此时${m_{\text{k}}} = {r_{{\text{k1}}}}$, 可得${f_{{\text{nd}}}} = 1$. 再结合(40)式可知, 此时WM-HRL模型为完全RANS模式. 另一方面, 当${r_{\text{k}}} < {r_{{\text{k1}}}}$时, 由(39)式可知, ${m_{\text{k}}}= {r_{\text{k}}} < r {r_{{\text{k1}}}}$, 此时${f_{{\text{ns}}}} = 1 - {f_{\text{s}}}$.

根据$n{r_{{\text{k2-Q}}}}$准则, 函数${f_{\text{s}}}$需要满足如下条件: 第一, 当${r_{{\text{k2}}}} < {m_{\text{k}}} < {r_{{\text{k1}}}}$时, ${f_{\text{s}}}$需要满足$0 < {f_{\text{s}}} < 1$; 第二, 当${m_{\text{k}}} \leqslant {r_{{\text{k2}}}}$时, ${f_{\text{s}}}$需要满足${f_{\text{s}}} = 1$; 第三, 当${m_{\text{k}}} = {r_{{\text{k1}}}}$时, ${f_{\text{s}}}$ 需要满足${f_{\text{s}}} = 0$.

能够同时满足上述3个条件的函数可用一个关于${m_{\text{k}}}$的指数函数表示, 具体如下:

$${f_{\text{s}}}({m_{\text{k}}}) = {\text{min}}\left[ {{\text{2exp}}\left( { - 11{\alpha ^{2}}} \right),{1}{.0}} \right],$$

(41)

$$\alpha =\text{min}\left(-0\text{.25,}\frac{1}{2({r}_{\text{k2}}-{r}_{\text{k1}})}{m}_{k}+\frac{{r}_{\text{k1}}-3{r}_{\text{k2}}}{4({r}_{\text{k2}}-{r}_{\text{k1}})}\right),$$

(42)

其中, ${r_{{\text{k1}}}} > {r_{{\text{k2}}}}$为自定义参数, 且取值在0和1之间.

下面证明, 在(42)式的条件下, 由(41)式构造的函数${f_{\text{s}}}$, 满足其所需的3条要求: 首先, 由(42)式可知, 当${m_{\text{k}}} = {r_{{\text{k1}}}}$时, $\alpha = - 0.75$, 由(41)式可知, 此时${f_{\text{s}}} = 0$; 其次, 当${r_{{\text{k2}}}} < {m_{\text{k}}} < {r_{{\text{k1}}}}$时, $- 0.75 < \alpha < - 0.25$, 由(41)式可知, 此时$0 < {f_{\text{s}}} < 1$; 第三, 当${m_{\text{k}}} \leqslant {r_{{\text{k2}}}}$时, $\alpha = - 0.25$, 由(41)式可知, 此时${f_{\text{s}}}=1$.

目前, 常用的一类WMLES模型为所谓的代数壁面模化大涡模拟(Alg WMLES) ^[34], 其基本思想是对Smagorinsky型LES的亚格子涡黏系数乘以一个阻尼系数$Fr$, 即:

$$\begin{split}&{v_{{\text{sgs}}}} = Fr{\left( {{C_{{\text{SMAG}}}}\varDelta } \right)^2}\left| {\left\langle S \right\rangle } \right|,\\ &Fr = 1 - \exp \left[ { - {{\left( {{{{y^ + }} \mathord{\left/ {\vphantom {{{y^ + }} A}} \right. } A}} \right)}^3}} \right],\end{split}$$

(43)

其中, ${v_{{\text{sgs}}}}$为亚格子涡黏系数; C_SMAG为模型常数, 取值在0.1—0.18之间; $A$为常数, 一般取为25; y⁺为无量纲壁面距离.

由(43)式可知, 当${y^ + }$约大于60时, 阻尼函数$Fr$趋于1, 此时Alg WMLES即为完全Smagorinsky型亚格子涡黏模型. 对${y^ + } \leqslant 60$的这个近壁区域, 正好为黏性底层和过渡层, 由于黏性底层很薄, 其范围约为${y^ + } \leqslant 10.0$, 此时Alg WLMES的亚格子涡黏系数${v_{{\text{sgs}}}} \approx 0.0$, 这相当于DNS. 在过渡层内, $Fr$从0开始增大到1, 此时Alg WLMES相当于准DNS. 由此可见, 对高雷诺数钝体绕流问题, 为了捕捉黏性底层及过渡层这两个区域内足够精细的湍流信息, 此时Alg WMLES所需的网格量几乎与DNS的网格量相当.

另一类常用的WMLES模型为所谓的壁面应力模化大涡模拟(WRMLES) ^[39], 该模型根据湍流边界层速度剖面的对数律来计算壁面剪切力, 并输入到LES边界网格作为边界条件. 在WRMLES中, 网格间距($\varDelta$)与局部边界层厚度($\delta$)通常需要满足条件$\delta /\varDelta \approx 20—30$. 这种较粗的近壁网格有可能会导致缺乏湍流应力的现象, 同时基于最近邻LES速度对壁面应力进行建模的壁面应力边界条件会增大边界层内部区域的总应力, 因此可能会致壁面应力被低估或高估, 进而发生LLM问题^[40].

WM-HRL模型与Alg WMLES和WRMLES模型均不相同, 该模型在黏性底层内一般为完全RANS解析模式, 在过渡层内的某个区域内为RANS/LES混合解析模式, 在对数律层区域则为完全LES模式. 由于RANS模式和RAN/LES混合解析模式所需网格数量远小于DNS和准DNS模式的网格量, 而且在混合函数(39)式—(42)式的双重保护下, 不仅可以避免MSD及GIS的问题, 同时也可避免LLM的问题. 因此, 该WM-HRL模型有望成为高雷诺数壁湍流三维时空瞬态发展湍流场高置信度解析的一种实用化的CFD工具.

本文利用作者团队自研CFD软件NUWA: FLOW_UV系统, 对如上所述WM-HRL模型开发了相应的植入式程序. 该自研CFD软件系统, 采用有限体积法(FVM)求解RANS方程. 其中, 压力速度求解采用PISO算法, 对流项及扩散项的空间离散采用二阶中心差分格式, 时间离散格式为二阶隐式格式.

3.   计算模型与设置

本文对雷诺数Re = 3900下圆柱绕流场问题, 采用WM-HRL进行数值模拟与分析. 其中, $Re = {{{U_\infty }D} \mathord{\left/ {\vphantom {{{U_\infty }D} v}} \right. } v}$, ${U_\infty }$为来流速度, D为圆柱直径. 计算区域为矩形区域, 如图1所示. 具体设置如下: 圆柱底面中心位于坐标原点$\left( {0, 0, 0} \right)$, 入口边界位于${{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} D}} \right. } D} = - 10$; 出口边界位于${{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} D}} \right. } D} = 15$; 前后两个边界分别位于${{{x_2}} \mathord{\left/ {\vphantom {{{x_2}} D}} \right. } D} = \pm 4$; 展向高度为${{{L_3}} \mathord{\left/ {\vphantom {{{L_3}} D}} \right. } D} = {\text{π }}$.

边界条件设置如下: 在入口边界, 速度入口设置为自由来流条件, 即$\left( {{{\bar u}_1}, {{\bar u}_2}, {{\bar u}_3}} \right) = \left( {{U_\infty }, 0, 0} \right)$; 湍动能按湍流强度$I = 5{\text{%}}$设定, 其中$k = {{3{{\left( {{U_\infty }I} \right)}^2}} / 2}$; 比耗率按$\omega = {k \mathord{\left/ {\vphantom {k {{v_{\text{t}}}}}} \right. } {{v_{\text{t}}}}}$设定, 其中${{{v_{\text{t}}}} \mathord{\left/ {\vphantom {{{v_{\text{t}}}} v}} \right. } v} = 10.0$. 在出口边界, 设置为零压梯度出口, 即$\nabla \bar p = 0$. 在两个垂直侧面及上下边界, 均设置为对称边界条件. 在圆球表面边界, 设置为无滑移条件, 即${\bar u_1} = {\bar u_2} = {\bar u_3} = 0$, 湍动能设置为$k = 0$, 比耗率设置为$\omega = {{6 v}/ {{\beta ^*}d_{\text{w}}^{2}}}$.

采用ANSYS ICEM软件进行分块网格划分, 如图2所示. 单元网格为六面体, 圆柱体壁面处第一层网格位于${y^ + } = 0.66$. 沿圆柱体展向网格的尺寸为${{{\varDelta _3}} \mathord{\left/ {\vphantom {{{\varDelta _3}} D}} \right. } D} = {{\text{π }} \mathord{\left/ {\vphantom {{\text{π }} {64}}} \right. } {64}}$, 水平方向的网格平均尺寸为${{{\varDelta _1}({\varDelta _2})} \mathord{\left/ {\vphantom {{{\varDelta _1}({\varDelta _2})} D}} \right. } D} = 0.02$. 自第一层起, 各层网格在径向的增长率为1.07, 在边界层内总计布置了45层网格, 整个计算域的网格数量为143万.

对SST k-ω模型, 在钝体近壁区通常采用所谓的壁面函数模型^[42]. 该壁面函数模型对湍动能k、比耗率ω及壁面剪切应力等, 在黏性底层(${y^ + } < 5$)均给出了明确的边界条件. 在本文所构建的WM-HRL模型中, 在钝体近壁区采用的是SST k-ω模型, 为使所构建的WM-HRL模型在钝体近壁区发挥出实际的壁面模化作用, 第一层网格一般需要设置在黏性底层的${y^ + } < 1$内.

对雷诺数Re = 3900下的圆柱绕流问题, 本文通过系列数值模拟研究表明, 将第一层网格设置在${y^ + } = 0.66$处, 可同时兼顾计算精度及网格总量控制等要求. 另外, Lehmkuhl等^[10]采用DNS的研究表明, 在湍流核心区Kolmogorov尺度的平均值为${{\overline \eta } \mathord{\left/ {\vphantom {{\overline \eta } D}} \right. } D} = 0.02$, 为保证网格密度能够捕捉到最小尺度之圆柱绕流的相干结构, 网格尺寸需要满足: ${{\overline \varDelta } \mathord{\left/ {\vphantom {{\overline \varDelta } {\overline \eta }}} \right. } {\overline \eta }} = 0.9$. 有鉴于此, 本文将径向的增长率设置为1.07.

经统计, 对雷诺数Re = 3900的情况, 以圆柱体展向高度${{{L_3}} \mathord{\left/ {\vphantom {{{L_3}} D}} \right. } D} = {\text{π }}$为对象进行CFD计算的主要文献如表1所示. 其中, Pereira等^[2]采用的是${{{L_3}} \mathord{\left/ {\vphantom {{{L_3}} D}} \right. } D} = 3.0$. 表中还列出了展向网格分辨率${{{\varDelta _3}} \mathord{\left/ {\vphantom {{{\varDelta _3}} D}} \right. } D}$以及所用总网格量等信息, 并与本文所采用的相关网格信息进行比较.

由表1可知, Lehmkuhl等^[10]采用的展向网格分辨率最高, 所用网格量也最大. Tremblay^[8]和Breuer^[15]采用的展向网格分辨率一致, 但水平方向的网格分辨率不同. Luo等^[24]采用的展向网格分辨率略低于前三篇文献的情况, 但水平方向的网格分辨率高于Breuer^[15]. Pereira等^[2]和D'Alessandro等^[30]采用的展向网格分辨率最低, 两者的总网格量接近. 本文所采用的展向网格分辨率与Tremblay^[8]和Breuer^[15]的一致, 但水平方向的网格分辨率均要低于这两篇文献的情况. 总之, 本文所用计算网格数量均少于DNS^[10], LES^[8,15], PANS^[2,24]及DES类模型^[24,30]的计算网格数量.

在数值计算中, 无量纲化时间步长$\varDelta {t}^{*}(= \varDelta t{U}_{\infty }/D)$取值为$6.8 \times {10^{ - 3}}$, 库朗数CFL<1, 计算时长为70个泄涡周期, 并对后45个泄涡周期的数据做统计平均, 以获取圆柱绕流场统计量等信息, 并与文献中相关实验结果及CFD数值模拟结果进行比较分析, 如表2所示.

在表2中, ${\bar f_{{\text{vs}}}}$为无因次涡泄频率, ${\bar f_{{\text{kh}}}}$为无因次K-H不稳定性的频率, ${\phi _{\text{s}}}$为分离角, ${{{L_{\text{r}}}} \mathord{\left/ {\vphantom {{{L_{\text{r}}}} D}} \right. } D}$为无因次回流区长度, ${C_{\text{d}}}$为阻力系数, $- {C_{{\text{pb}}}}$为基线压力系数, 所谓“形状”, 指的是在站位${{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} D}} \right. } D} = 1.06$处平均流向速度的剖面形状. 其中, 无因次频率$\bar f$定义为$\bar f={{fD} \mathord{\left/ {\vphantom {{fD} {{U_\infty }}}} \right. } {{U_\infty }}}$, $f$为相应的有因次频率.

在表2中, 对PANS模型, 滤波控制函数f_k定义为f_k = k_u/k. Lehmkuhl等^[10]采用了两种DNS模式. 其中, “Mode H”表示“高能模态”, 利用该模态的DNS得到V形结果. 而“Mode L”表示“低能模态”, 利用该模态的DNS得到U形结果.

由表2可知, 对Pereira等^[2] (PANS) (f_k = 0.25, 0.5)的情况, 与表中所列实验结果^[18,27]相比, 分离角的计算结果均明显偏小. 对Luo等(PANS)^[24] (${f_{\text{k}}}$ = 0.5)的情况, 分离角${\phi _{\text{s}}}$、阻力系数${C_{\text{d}}}$及基线压力系数$- {C_{{\text{pb}}}}$的计算结果与表中所列实验结果结果^[27]相比均明显偏大, 相对误差分别达到9.176%, 37.76%和63.333%. 对D'Alessandro等(SA-DES) ^[30]的情况, 阻力系数${C_{\text{d}}}$的计算结果与表中所列实验结果^[27]相比均明显偏大, 相对误差达22.7%.

由表2可进一步发现, 在Parnaudeau等^[18]和Lourenco和Shih^[27]的实验文章中未给出剪切层小尺度K-H不稳定性的频谱特征及其无因次频率${\bar f_{{\text{kh}}}}$的实验结果, 而在Tremblay^[8], Breuer^[15], Luo等^[24]及D'Alessandro等^[30]的CFD计算中, 也未给出剪切层小尺度K-H不稳定性的频谱特征及其无因次频率${\bar f_{{\text{kh}}}}$的计算结果, 但Lehmkuhl等^[10]和Pereira等^[2]给出了剪切层小尺度K-H不稳定性的频谱特征及其无因次频率${\bar f_{{\text{kh}}}}$的计算结果.

Parnaudeau等^[18]的实验测得回流区长度为${{{L_{\text{r}}}} \mathord{\left/ {\vphantom {{{L_{\text{r}}}} D}} \right. } D}$ = 1.51, 并得到U形流向速度剖面. Lourenco和Shih^[27]的实验测得回流区长度为${{{L_{\text{r}}}} \mathord{\left/ {\vphantom {{{L_{\text{r}}}} D}} \right. } D}$ = 1.18, 并得到V形流向速度剖面. Lehmkuhl等^[10]利用DNS之Mode H计算得到的回流区长度${{{L_{\text{r}}}} \mathord{\left/ {\vphantom {{{L_{\text{r}}}} D}} \right. } D}$ = 1.26, 与Lourenco和Shih^[27]实验结果的相对误差为6.78%. 而Lehmkuhl等^[10]利用DNS之Mode L计算得到回流区长度${{{L_{\text{r}}}} \mathord{\left/ {\vphantom {{{L_{\text{r}}}} D}} \right. } D}$ = 1.55, 与Parnaudeau等^[18]实验结果的相对误差为2.65%.

Tremblay^[8]及Breuer ^[15]采用LES均得到V形剖面, 但是他们的CFD计算所得回流区长度与Lourenco和Shih^[27]实验测得的回流区长度有较大差异, 相对误差分别达到–11.86%和16.27%. Pereira等^[2]采用PANS在${f_{\text{k}}}$ = 0.25, 0.5下分别得到U形和V形剖面. 当${f_{\text{k}}}$ = 0.5时, CFD计算所得回流区长度与Lourenco和Shih^[27]实验测得的回流区长度接近, 相对误差为–5.08%. 当${f_{\text{k}}}$ = 0.25时, CFD计算所得回流区长度与Parnaudeau等^[18]实验测得的回流区长度相比偏大, 相对误差达14.56%. 此外, 对Pereira等^[2] (PANS) (${f_{\text{k}}}$ = 0.5)的情况, 与表中所列DNS^[10]结果相比, K-H不稳定性频率${\bar f_{{\text{kh}}}}$的计算结果明显偏大, 相对误差达15.67%.

Luo等^[24]采用PANS在${f_{\text{k}}}$ = 0.1和0.5下均得到V形剖面. 当${f_{\text{k}}}$ = 0.1时, CFD计算所得回流区长度与Lourenco和Shih^[27]实验测得的回流区长度的相对误差为7.63%. 当${f_{\text{k}}}$ = 0.5时, CFD计算所得回流区长度与Lourenco 等^[27]实验测得的回流区长度相比明显不符. Luo等^[24]采用SST-DES得到V形剖面, 但计算所得回流区长度与Lourenco和Shih^[27]实验测得的回流区长度相比也明显不符.

D'Alessandro等^[30]采用SA-DES也得到V形剖面, 但计算所得回流区长度与Lourenco和Shih^[27]实验测得的回流区长度相比明显不符. D'Alessandro等^[30]采用SA-IDDES及v²-f DES均得到U形剖面, 采用SA-IDDES计算得到的回流区长度与Parnaudeau等^[18]实验测得的回流区长度的相对误差为–5.50%. 但采用v²-f DES计算得到的回流区长度与Parnaudeau等^[18]实验测得的回流区长度相比偏大, 相对误差高达11.13%.

对钝体绕流问题, 其边界层包括黏性底层、过渡层及对数律层. 对本文WM-HRL模型来讲, 这三类边界层范围的信息将是至关重要的. 为此, 利用图2所示计算网格系统, 首先进行RANS数值模拟, 结果表明: 在雷诺数Re = 3900下圆柱绕流场的边界层厚度为${y^ + } \leqslant 180.0$, 黏性底层厚度为${y^ + } \leqslant 10.0$, 过渡层厚度为$10.0 < {y^ + } < 30.0$.

对雷诺数Re = 3900下的圆柱绕流场问题, 其相干结构主要包含两类结构. 其中, 第一类为与尾流中涡泄相关的大尺度不稳定性结构, 其无因次特征频率为${\bar f_{{\text{vs}}}}$, 而第二类为与分离剪切层脉动相关的小尺度不稳定性结构, 称为K-H不稳定性, 其无因次特征频率为${\bar f_{{\text{kh}}}}$. 对本文WM-HRL模型来讲, 圆柱绕流剪切层小尺度K-H不稳定性结构发生的位置信息将是至关重要的.

剪切层的速度梯度较大, 属于一个位置狭小的窄带, 且对监测点的分布位置敏感^[23]. 为获取圆柱绕流剪切层小尺度K-H不稳定性结构发生的准确边界位置, 结合相关文献中DNS^[10], PANS^[2]及LES^[23]等CFD数值模拟结果, 本文设置了如图3所示的18个监测点, 其坐标位置如表3所示. 其中, 监测点P1—P14位于剪切层区域, 监测点P15位于剪切层脱落区, 而监测点P16—P18位于尾流中心线上. 表3中所列的这些监测点均位于圆柱体底部位置. 同时, 对表3中所列的每个监测点, 沿圆柱体展向同时均匀布置了其他4个相应的监测点. 因此, 实际上总计布置了90个监测点.

在此基础上, 利用本文WM-HRL进行系列数值模拟, 对所得各监测点处的横向脉动速度进行Lomb频谱分析, 详细分析见4.3节. 在做Lomb谱分析时, 用的是表3中所列各监测点相对应的5个监测点处横向脉动速度做展向平均所得, 其做法与Lehmkuhl等^[10](DNS)的做法相同.

结果表明, 在各个监测点的Lomb频谱中, 均出现一个频率相对较低的谱峰, 该谱峰所对应频率与涡泄频率一致. 同时, 对各种工况的系列CFD数值模拟结果中, 在测点P4的Lomb频谱中, 均观察到一个频率相对较高的谱峰, 该谱峰所对应频率与剪切层小尺度K-H不稳定性结构的频率一致. 由此可见, 雷诺数Re = 3900下圆柱绕流场中剪切层小尺度K-H不稳定性结构发生在对数律层中${y^ + } \geqslant 89.4$的区域.

在图2所述网格系统设置下, 利用(36)式进行计算, 结果表明: 在对数律层的${y^ + } \geqslant 67.4$的区域中, ${r_{\text{k}}}$均不大于0.2, 即${r_{\text{k}}} \leqslant 0.2$. 由于圆柱绕流剪切层小尺度K-H不稳定性结构发生在${y^ + } \geqslant 89.4$的区域, 因此本文WM-HRL模型对该剪切层小尺度K-H不稳定结构为完全LES解析模式, 且其对湍动能的解析度至少为80%.

4.   数值结果与分析

当利用WM-HRL进行数值模拟时, 影响圆柱绕流场计算置信度的主要因素包括两个边界位置和3个区域的湍动能解析度. 其中, 两个边界分别为RANS结束边界(记为${\varGamma _{{\text{RANS}}}}$)和LES启动边界(记为${\varGamma _{{\text{LES}}}}$), 3个区域分别为RANS区、RANS/LES混合转换区和LES区. 为下文陈述简便计, 记$y_{{\text{RANS}}}^ +$为RANS结束边界${\varGamma _{{\text{RANS}}}}$的位置, 而记$y_{{\text{LES}}}^{+}$为LES启动边界${\varGamma _{{\text{LES}}}}$的位置.

4.1   WM-HRL模型参数影响规律

首先讨论LES启动边界${\varGamma _{{\text{LES}}}}$位于剪切层小尺度K-H不稳定性结构发生区域的情况, 并考虑两种情况: 第一种情况为$y_{{\text{LES}}}^{+}=105.8$, ${r_{{\text{k2}}}} = 0.{1556}$; 第二种情况为$y_{{\text{LES}}}^{+}=113.9$, ${r_{{\text{k2}}}} = 0.{1484}$. 在每种情况下, 均设置11种RANS结束边界${\varGamma _{{\text{RANS}}}}$的位置及其相应的${r_{{\text{k}}1}}$, 如表4所示. 其中, $y_{{\text{RANS}}}^{+}{\text{ = 7}}{.9}$位于黏性底层, $y_{{\text{RANS}}}^{+}{\text{ = 13}}{.3}, \;{14}{.9}, \;{18}{.4}, \;{20}{.4}$, 27.1位于过渡层, $y_{{\text{RANS}}}^ + = 38.4, \;41.7, \;49.2$, 72.7位于对数律层但在剪切层K-H不稳定区外, 而$y_{{\text{RANS}}}^ + =91.2$则位于剪切层K-H不稳定性结构发生区域内.

在各种$y_{{\text{RANS}}}^ +$及${r_{{\text{k1}}}}$和$y_{{\text{LES}}}^{+}$及${r_{{\text{k2}}}}$的组合下, 利用WM-HRL进行数值模拟所得相关流场统计量的结果如表4所示. 由表4可知, 只有当$y_{{\text{LES}}}^{+}$ = 105.8且$y_{{\text{RANS}}}^{+}{\text{ = 14}}{.9}$时, 才能同时准确获取剪切层小尺度K-H不稳定性结构的无因次频率${\bar f_{{\text{kh}}}}$及回流区长度${{{L_{\text{r}}}} \mathord{\left/ {\vphantom {{{L_{\text{r}}}} D}} \right. } D}$的值. 在此情况下, 在站位${x_1}/D = 1.06$处流向平均速度的剖面形状为V形, 与Lourenco和Shih^[27]的实验结果一致. ${{{L_{\text{r}}}} \mathord{\left/ {\vphantom {{{L_{\text{r}}}} D}} \right. } D}$的计算结果与Lourenco和Shih^[27]实验结果相比, 相对误差为0.85%. ${\bar f_{{\text{kh}}}}$的计算结果与Lehmkuhl等^[10]的DNS之Mode H计算结果相比, 相对误差为0.75%.

下面讨论LES启动边界${\varGamma _{{\text{LES}}}}$位于剪切层小尺度K-H不稳定性结构发生区域外, 但仍位于对数律层区的情况, 并考虑3种情况: 第一种情况为$y_{{\text{LES}}}^{+}{\text{ = 49}}{.2}$, ${r_{{\text{k}}2}} = 0.{2546}$; 第二种情况为$y_{{\text{LES}}}^{+}= 72.4$, ${r_{{\text{k2}}}} = 0.{1983}$; 第三种情况为$y_{{\text{LES}}}^ + {\text{ = 84}}{.6}$, ${r_{{\text{k2}}}} = 0.{1713}$. 在每种情况下, 均设置若干种RANS结束边界${\varGamma _{{\text{RANS}}}}$的位置及其相应的${r_{{\text{k1}}}}$, 其设置原则与表4一致, 如表5所示. 在各种$y_{{\text{RANS}}}^{+}$及${r_{{\text{k1}}}}$和$y_{{\text{LES}}}^{+}$及${r_{{\text{k2}}}}$的组合下, 利用WM-HRL进行数值模拟所得相关流场统计量的结果如表5所示.

由表5可知, 只有当RANS结束边界${\varGamma _{{\text{RANS}}}}$位于过渡层, 且${r_{{\text{k}}1}} \leqslant {\bar r_{{\text{k}}1}} = 0.63$时, 才能同时准确获得${\bar f_{{\text{kh}}}}$和${{{L_{\text{r}}}}/D}$的计算结果. 在此条件下, 当计算得到V形速度剖面时, ${{{L_{\text{r}}}} \mathord{\left/ {\vphantom {{{L_{\text{r}}}} D}} \right. } D}$的计算结果与Lourenco和Shih^[27]的实验结果相比, 相对误差最大为 –9.32%, 而${\bar f_{{\text{kh}}}}$的计算结果与Lehmkuhl等^[10]的DNS之Mode H计算结果相比, 相对误差最大为12.68%. 当计算得到U形速度剖面时, ${{{L_{\text{r}}}}/ D}$的计算结果与Parnaudeau等^[18]的实验结果相比, 相对误差最大为–12.58%, 而${\bar f_{{\text{kh}}}}$的计算结果与Lehmkuhl等^[10]的DNS之Mode H计算结果相比, 相对误差最大为–20.15%.

最后讨论LES启动边界${\varGamma _{{\text{LES}}}}$位于过渡层的情况, 并考虑5种情况, 分别为$y_{{\text{LES}}}^{+}{\text{ = 10}}{.4}$, $y_{{\text{LES}}}^{+}=13.3$, $y_{{\text{LES}}}^{+}{\text{ = 18}}.4$, $y_{{\text{LES}}}^{+}=20.4$和$y_{{\text{LES}}}^{+}{\text{ = 29}}{.6}$. 在每种情况下, 均设置若干种RANS结束边界${\varGamma _{{\text{RANS}}}}$的位置$y_{{\text{RANS}}}^ +$及其相应的${r_{{\text{k1}}}}$. 在各种$y_{{\text{RANS}}}^ +$及$y_{{\text{LES}}}^{+}$的组合下, 利用WM-HRL进行数值模拟所得相关流场统计量的结果如表6所示.

首先, 利用本文WM-HRL模型所得${\bar f_{{\text{vs}}}}$, ${\phi _{\text{s}}}$, ${C_{\text{d}}}$及$- {C_{{\text{pb}}}}$的计算结果之置信度进行分析. 由表2可知, 文献[18, 27]中实验及CFD计算^{[2,8,10,15,18,24,27,30]}所得${\bar f_{{\text{vs}}}}$之值在2.01—2.22范围内, 所得${\phi _{\text{s}}}$之值在85—89.28范围内, 所得${C_{\text{d}}}$之值在0.927—1.14范围内, 而所得$- {C_{{\text{pb}}}}$之值在0.829—1.077范围内.

需要指出的是, 在上文统计中, 对${\phi _{\text{s}}}$之CFD计算值, 已去掉了Pereira等^[2] (PANS) (${f_{\text{k}}}$ = 0.25和0.5)的异常计算结果, 以及Luo等^[24](PANS)(${f_{\text{k}}}$ = 0.5)的异常计算结果. 对${C_{\text{d}}}$之CFD计算值, 已去掉Luo等^[24](PANS)(${f_{\text{k}}}$ = 0.5)的异常计算结果, 以及D'Alessandro等^[30](SA-DES)的异常计算结果. 对$- {C_{{\text{pb}}}}$之CFD计算值, 已去掉Luo等^[24](PANS)(${f_{\text{k}}}$ = 0.5)的异常计算结果.

综合表4—表6可知, 利用本文WM-HRL模型, 所得${\bar f_{{\text{vs}}}}$的计算值在0.217—0.226范围内, 所得${\phi _{\text{s}}}$之计算值在86.4^o—88.4^o范围内, 所得${C_{\text{d}}}$之CFD计算值在1.03—1.17范围内, 而所得$- {C_{{\text{pb}}}}$之计算值在0.95—1.15范围内. 由此可见, 本文利用WM-HRL模型所得${\bar f_{{\text{vs}}}}$, ${\phi _{\text{s}}}$, ${C_{\text{d}}}$及$- {C_{{\text{pb}}}}$之计算结果, 与文献[18, 27]中实验及CFD计算^{[2,8,10,15,18,24,27,30]}所得相应结果具有良好的一致性.

其次, 就利用本文WM-HRL模型所得${\bar f_{{\text{kh}}}}$及${{{L_{\text{r}}}}/D}$的计算结果之置信度进行分析. 由表6可知, 利用本文WM-HRL模型所得${\bar f_{{\text{kh}}}}$的计算值与Lehmkuhl等^[10]利用Mode H之DNS计算结果相比, 相对误差最小为–1.49%, 最大为10.44%. 对利用本文WM-HRL模型所得${{{L_{\text{r}}}}/D}$的计算结果, 当计算得到V形速度剖面时, ${{{L_{\text{r}}}}/D}$之计算值与Lourenco和Shih^[27]的实验结果相比, 相对误差最小为–8.47%, 最大为10.0%. 当计算得到U形速度剖面时, ${{{L_{\text{r}}}}/D}$之计算值与Parnaudeau等^[18]的实验结果相比, 相对误差最小为–5.96%, 最大为–9.93%.

由于剪切层小尺度K-H不稳定性发生的位置是未知的, 但位于对数律层区中. 系列数值模拟结果表明, 对本文所构造的WM-HRL模型, 为了能够准确解析并捕捉剪切层小尺度K-H不稳定性的结构特征及其频谱特性, 可将其LES模式的启动边界均设置在过渡层内, 而RANS模式的结束边界则既可设置在过渡层内也可设置在黏性底层内, 并且使${r_{{\text{k2}}}} \leqslant 0.2$, 即在对数律层区的网格具有至少80%的湍动能解析度.

在圆柱体展向长度及其网格系统已经确定的条件下, 对DNS^[6-12]及LES^{[6–8,13-23]}来讲, 通过CFD计算只能得到回流区长度及流向速度剖面形状的其中一种分支结构. 对DES类模型^[24,30]来讲, 不同RANS湍流模型的类型及RANS/LES之间不同的混合转换方式等, 均可能会影响其对回流区长度及流向速度剖面形状分支结构的计算结果.

对PANS模型^[2,24–26]而言, 其核心思想是通过引入一个${f_{\text{k}}}$参数来调控流场可解/不可解湍流量的比例而实现. 对基于SST k-ω的PANS模型, ${f_{\text{k}}}$参数可调控其ω方程中耗散项之值的变化. 一般地, ${f_{\text{k}}}$值越小, 耗散项的值也越小, 从而使求解得到的比耗率ω增大, 进而使湍流涡黏系数${v_{\text{t}}}$减小, 可解尺度就释放的越多^[25].

对PANS模型, Pereira等^[2]指出: 为准确获取圆柱绕流回流区结构及其长度信息, ${f_{\text{k}}}$的值不能大于0.5, 即${f_{\text{k}}} \leqslant 0.5$. 同时, 当${f_{\text{k}}}$从0.5减小到某个值(文中为0.25)后, 圆柱尾部近壁面处流向速度剖面的形状从V形转变为U形. 但Luo等^[24]在${f_{\text{k}}}=0.5$时通过CFD计算并没有得到准确的回流区长度, 而且在${f_{\text{k}}}=0.1$时, 得到的是V形剖面. 导致这种相矛盾结果的原因之一, 可能与两位作者所使用网格结构及其空间分辨率分布不同有关(具体见表2).

本文采用WM-HRL模型的系列数值模拟结果表明, 在同一套网格系统下, 通过改变WM-HRL模型中两个边界位置及3个区域的湍动能解析度的组合, 既可以通过数值模拟获得与Parnaudeau等^[18]实验结果一致的回流区长度及圆柱近壁面处平均流向速度的U形剖面, 也可以通过数值模拟获得与Lourenco和Shih^[27]实验结果一致的回流区长度及圆柱近壁面处平均流向速度的V形剖面.

这表明对本文所建立的WM-HRL模型, 可以在同一套网格系统下通过变化湍动能解析度指标参数$n{r_{{\text{k1-Q}}}}$和$n{r_{{\text{k2-Q}}}}$准则的组合条件, 精细地解析圆柱绕流场中两类不同回流区长度结构特征, 及其对应的圆柱尾部近壁面处U形和V形两个流向速度剖面的分支结构.

4.2   一阶和二阶统计量特性

选取表6中相关数值模拟结果案列, 讨论雷诺数Re = 3900下圆柱绕流场的一阶和二阶统计量特性的计算结果, 并与Parnaudeau等^[18]及Lourenco和Shih^[27]的相关实验结果进行比较分析. 对U形流向速度剖面的情况, 选取4个组合工况如下. Case AU: $y_{{\text{RANS}}}^ +$ = 10.4, $y_{{\text{LES}}}^{+}$ = 18.4; Case BU: $y_{{\text{RANS}}}^ +$ = 14.9, $y_{{\text{LES}}}^{+}$ = 18.4; Case CU: $y_{{\text{RANS}}}^ +$ = 7.9, $y_{{\text{LES}}}^{+}$ = 20.4; Case DU: $y_{{\text{RANS}}}^ +$ = 27.1, $y_{{\text{LES}}}^{+}$ = 29.6. 其中, 对Case AU工况, 其RANS结束边界位于黏性底层和过渡层交界处, 而LES启动边界位于过渡层内; 对Case BU工况, 其RANS结束边界和LES启动边界均位于过度层内; 对Case CU工况, 其RANS结束边界位于黏性底层内, 而LES启动边界位于过渡层内; 对Case DU工况, 其RANS结束边界位于过渡层内, 而LES启动边界位于过渡层与对数律层交界处.

对V形流向速度剖面的情况, 选取4个组合工况如下. Case AV: $y_{{\text{RANS}}}^ +$ = 7.9, $y_{{\text{LES}}}^{+}$ = 13.3; Case BV: $y_{{\text{RANS}}}^ +$ = 10.4, $y_{{\text{LES}}}^{+}$ = 13.3; Case CV: $y_{{\text{RANS}}}^ +$ = 14.9, $y_{{\text{LES}}}^{+}$ = 20.4; Case DV: $y_{{\text{RANS}}}^ +$ = 20.4, $y_{{\text{LES}}}^{+}$ = 29.6. 其中, 对Case AV工况, 其RANS结束边界位于黏性底层内, 而LES启动边界位于黏性底层与过渡层交界处; 对Case BV工况, 其RANS结束边界位于黏性底层和过渡层交界处, 而LES启动边界位于过渡层内; 对Case CV工况, 其RANS结束边界和LES启动边界均位于过渡层内; 对Case DV工况, 其RANS结束边界位于过渡层内, 而LES启动边界位于过渡层与对数律层交界处.

在针对一阶和二阶统计量的分析中, 对Case AU—Case DU这4种工况, 除了将其与Parnaudeau等^[18]的实验结果进行比较外, 同时还与Lehmkuhl等^[10]利用Mode H之DNS计算结果进行了比较. 对Case AV—Case DV这4种工况, 除了将其与Lourenco和Shih^[27]的实验进行比较外, 同时还与Lehmkuhl等^[10]利用Mode L之DNS计算结果进行了比较.

4.2.1   一阶统计量特性

在图4中, 分别给出了上述8种工况下圆柱表面周向系数${C_{\text{p}}}$分布特性的数值结果. 由于Parnaudeau等^[18]和Lourenco和Shih^[27]的实验没有给出圆柱表面压力系数的数据, 此处采用Norberg^[43]的实验结果进行比较分析.

由图4可知, 无论是对Case AU—Case DU这4种工况, 还是对Case AV—Case DV这4种工况, 利用本文WM-HRL模型所得沿圆柱表面周向压力系数${C_{\text{p}}}$分布的数值结果之间的差异均很小. 对图4(a)的情况, 本文计算结果与Norberg^[43]实验结果相比的相对误差最大为10.3%, 而Lehmkuhl等^[10]利用DNS(Mode L)计算结果与Norberg^[43]实验结果相比的相对误差最大为3.1%. 对图4(b)的情况, 本文计算结果与Norberg^[43]实验结果相比的相对误差最大为14.0%, 而Lehmkuhl等^[10]利用DNS(Mode H)计算结果Norberg^[43]实验结果相比的相对误差最大为0.6%.

在图5中, 分别给出了前述8种工况下沿尾流中心线平均流向速度剖面特性的数值结果. 由图5可知, 在圆柱表面正后方中心线上的点P(0.5D, 0)处, 平均流向速度的值为0, 在达到一个负的最小值后开始增大, 并在中心线上的点Q(${L_{\text{r}}}{+}0.5 D$, 0)处又变为0. 其中, ${{{L_{\text{r}}}} / D}$为回流区长度. Parnaudeau等^[18]的实验测量获得${{{L_{\text{r}}}} /D} = 1.51$ , Lehmkuhl等^[10]利用Mode L之DNS计算得到${{{L_{\text{r}}}} /D}=1.55$, 如图5(a)所示. 而Lourenco和Shih^[27]的实验测量获得${{{L_{\text{r}}}} \mathord{\left/ {\vphantom {{{L_{\text{r}}}} D}} \right. } D}$ = 1.18, Lehmkuhl等^[10]利用Mode H之DNS计算得到${{{L_{\text{r}}}} \mathord{\left/ {\vphantom {{{L_{\text{r}}}} D}} \right. } D}$ = 1.26, 如图5(b)所示. 此即在雷诺数Re = 3900下圆柱绕流两个实验及DNS计算中出现两类不同回流区分支结构的情况.

由图5进一步可发现, 无论是对Case AU—Case DU这4种工况, 还是对Case AV—Case DV这4种工况, 利用本文WM-HRL模型所得沿尾流中心线平均流向速度分布的数值结果之间的差异均很小. 对图5(a)的情况, 在回流区外, 本文数值结果与Parnaudeau等^[18]的实验结果之间的相对误差最大为5.5%, 而与Lehmkuhl等^[10]利用Mode L之DNS的计算结果之间的相对误差最大为3.2%.

对图5(b)的情况, 本文数值结果与Lehmkuhl等^[10]利用Mode H之DNS的计算均吻合良好, 两者之间的相对误差最大为4.2%. 但无论是本文计算结果还是Lehmkuhl等^[10]利用Mode H之DNS的计算结果, 在${{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} D}} \right. } D} \in$[2.36, 3.26]的范围内, 它们与Lourenco和Shih^[27]的实验结果均有一定的差异.

图6中, 分别给出了前述8种工况下在3个不同站位(${{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} D}} \right. } D}$ = 1.06, 1.54, 2.02)处平均流向速度剖面特性的数值结果. 站位${{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} D}} \right. } D}$ = 1.06位于剪切层K-H不稳定性结构发生区域, 同时也位于回流区内. Parnaudeau等^[18]的实验测量获得U形速度剖面, 而Lehmkuhl等^[10]利用Mode L之DNS计算也得到U形剖面, 如图6(a)所示. 而Lourenco和Shih^[27]的实验测量获得V形速度剖面, 而Lehmkuhl等^[10]利用Mode H之DNS计算也得到V形剖面, 如图6(b)所示. 此即在雷诺数Re = 3900下, 圆柱绕流两个试验中所测量得到的平均流向速度剖面出现U形和V形两个流动分支结构的情况.

由图6可知, 无论是对Case AU—Case DU这4种工况, 还是对Case AV—Case DV这4种工况, 利用本文WM-HRL模型所得该站位处平均流向速度分布的数值结果之间的差异均很小. 在位于剪切层K-H不稳定性结构发生区域的站位${{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} D}} \right. } D}$ = 1.06处, 对Case AU—Case DU这4种工况, 本文计算所得平均流向速度剖面均为U形, 与Parnaudeau等^[18]的实验及Lehmkuhl等^[10]利用Mode L之DNS计算所得速度剖面形状均一致. 对Case AV—Case DV这4种工况, 本文计算所得平均流向速度剖面均为V形, 与Lourenco和Shih^[27]的实验及Lehmkuhl等^[10]利用Mode H之DNS计算所得速度剖面形状均一致.

站位${{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} D}} \right. } D}$ = 1.54和2.02位于回流区中. 由图6可见, Parnaudeau等^[18]和Lourenco和Shih^[27]的实验结果均获得V形流向速度剖面, 而Lehmkuhl等^[10]利用Mode L和H之DNS计算也均获得V形流向速度剖面. 同时, 由图6还可观察到, 对Case AU-Case DU这4种工况, 利用本文WM-HRL模型计算所得平均流向速度剖面均为V形, 与Parnaudeau等^[18]的实验及Lehmkuhl等^[10]利用Mode L之DNS计算所得速度剖面形状一致. 对Case AV—Case DV这4种工况, 利用本文WM-HRL模型计算所得平均流向速度剖面也均为V形, 与Lourenco和Shih^[27]的实验及Lehmkuhl等^[10]利用Mode H之DNS计算所得速度剖面形状均一致.

图7分别给出了前述8种工况下在3个不同站位(${{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} D}} \right. } D}$ = 1.06, 1.54, 2.02)处平均横向速度剖面特性的数值结果.

由图7可见, 无论是对Case AU—Case DU这4种工况, 还是对Case AV—Case DV这4种工况, 利用本文WM-HRL模型所得圆柱后方不同站位处平均横向速度分布的数值结果差异均很小.

对图7(a)的情况, 在${{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} D}} \right. } D}$ = 1.06站位处, 本文计算结果与Parnaudeau等^[18]实验结果之间的相对误差最大为3.5%, 而与Lehmkuhl等^[10]利用Mode L的DNS计算结果之间的相对误差最大为3.8%. 在${{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} D}} \right. } D}$ = 1.54站位处, 本文计算结果与Parnaudeau等^[18]实验结果之间的相对误差最大为9.9%, 而与Lehmkuhl等^[10]利用Mode L的DNS计算结果之间的相对误差最大为10.3%. 在${{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} D}} \right. } D}$ = 2.02站位处, 本文计算结果与Parnaudeau等^[18]实验结果之间的相对误差最大为15%, 而与Lehmkuhl等^[10]利用Mode L之DNS计算结果之间的相对误差最大为16.1%.

对图7(b)的情况, 本文计算结果与Lehmkuhl等^[10]利用Mode H之DNS计算结果一致, 但两者与Lourenco和Shih^[27]的实验结果均有一定的差异. 在${{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} D}} \right. } D}$ = 1.06站位处, 本文计算结果与Lehmkuhl等^[10]利用Mode H的DNS计算结果之间的相对误差最大为7.8%. 在${{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} D}} \right. } D}$ = 1.54站位处, 本文计算结果与Lehmkuhl等^[10]利用Mode H之DNS计算结果之间的相对误差最大为5.4%. 在${{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} D}} \right. } D}$ = 2.02站位处, 本文计算结果与Lehmkuhl等^[10]利用Mode H之DNS计算结果之间的相对误差最大为10.0%.

Tremblay^[8]采用DNS和LES均复现了Lourenco和Shih^[27]对雷诺数Re = 3900下, 圆柱绕流在其后方不同站位处平均流向及横向速度剖面特性的实验结果. 其中, LES亚格子模型分别采用了经典Smagorinsky模型(Smag)及动态模型(Dyn). 结果表明, DNS的结果要优于LES的结果, 而Dyn亚格子模型的结果则要优于Smag亚格子模型的结果.

对PANS模型, 其对圆柱后方各站位处平均流向及横向速度剖面的计算精度依赖于滤波控制参数${f_{\text{k}}}$的取值方式. 对DES类模型, 其对圆柱后方各站位处平均流向及横向速度剖面的计算精度则依赖于所用RANS模型的类型及RANS/LES的混合转换方法等. 总体上, 在网格空间分辨率分布及相关的模型参数设置合理的条件下, 这两类模型对圆柱后方各站位处平均流向及横向速度剖面的计算精度, 一般地可以达到与LES相当的计算精度^[2,24-26,30].

本文采用WM-HRL模型的系列数值模拟结果表明, 在同一套网格系统下, 通过改变WM-HRL模型中两个边界位置及3个区域的湍动能解析度的组合, 既可以通过数值模拟获得与Parnaudeau等^[18]实验所得一致的圆柱后方各站位处的平均流向及横向速度剖面, 也可以通过数值模拟获得与Lourenco和Shih^[27]实验所得一致的圆柱后方各站位处的平均流向及横向速度剖面.

4.2.2   二阶统计量特性

图8分别给出了前述8种工况下在3个不同站位(${{{x_1}} /D}$ = 1.06, 1.54, 2.02)处各向同性流向雷诺应力剖面特性的数值结果.

由图8可见, 对Case AU—Case DU的情况, 在位于剪切层小尺度K-H不稳定性结构发生区域的站位${{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} D}} \right. } D}$ = 1.06处, 本文计算结果与Parnaudeau等^[18]的实验结果及Lehmkuhl等^[10]利用Mode L之DNS计算结果相比, 主要差异在两个“峰”处. 在位于回流区的站位${{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} D}} \right. } D}$ = 1.54, 2.02处, 本文计算结果与Parnaudeau等^[18]的实验结果及Lehmkuhl等^[10]利用Mode L的DNS计算结果相比, 在两个“峰”处的值小于相应的实验值及DNS计算值.

对Case AV—Case DV的情况, 本文计算结果均与Lourenco和Shih^[27]实验结果及Lehmkuhl等^[10]利用Mode H之DNS计算结果具有良好的一致性. 在位于剪切层小尺度K-H不稳定区的站位${{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} D}} \right. } D}$ = 1.06处, 本文对BV和DV工况的计算结果与Lourenco和Shih^[27]实验结果之间的最大相对误差分别为1.7%和5%, 而与Lehmkuhl等^[10]利用Mode H之DNS计算结果之间的最大误差分别为0.3%和3.7%. 在${{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} D}} \right. } D}$ = 1.54站位处, 本文计算结果与Lourenco和Shih^[27]实验结果之间的最大相对误差分别为16.1%, 3.4%, 3.8%, 6.4%. 在${{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} D}} \right. } D}$ = 2.02站位处, 本文计算结果与Lourenco和Shih^[27]实验结果之间的最大相对误差分别为4.6%, 11%, 3.5%, 6.7%.

图9分别给出了前述8种工况下在3个不同站位(${{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} D}} \right. } D}$ = 1.06, 1.54, 2.02)处各向同性横向雷诺应力剖面特性的数值结果.

由图9可见, 对Case AU—Case DU的情况, 在位于剪切层小尺度K-H不稳定区的站位${{{x_1}} /D}$ = 1.06处, 本文计算结果与Parnaudeau等^[18]实验结果相比的最大相对误差分别为23.6%, 7.7%, 15.7%, 9.1%, 而Lehmkuhl等^[10]利用Mode L之DNS计算结果与Parnaudeau等^[18]实验结果相比的最大相对误差为14.5%.

在位于回流区的站位${{{x_1}}/D}$ = 1.54处, 本文计算结果与Parnaudeau等^[18]实验结果相比的最大相对误差分别为13.7%, 25.3%, 6.4%, 13.6%, 而Lehmkuhl等^[10]利用Mode L之DNS计算结果与Parnaudeau等^[18]实验结果相比的最大相对误差为15%.

在位于回流区的站位${{{x_1}}/D}$ = 2.02处, 本文计算结果与Parnaudeau等^[18]实验结果相比的最大相对误差分别为3.3%, 7.2%, 12.9%, 0.5%, 而Lehmkuhl等^[10]利用Mode L的DNS计算结果与Parnaudeau等^[18]实验结果相比的最大相对误差为7%.

对Case AV—Case DV的情况, 无论是本文计算结果还是Lehmkuhl等^[10]利用Mode H的DNS计算结果一致, 两者均与Lourenco和Shih^[27]的实验结果有一定的差异. 在位于剪切层小尺度K-H不稳定区的站位${{{x_1}} /D}$ = 1.06处, 本文对Case BV和Case DV工况的计算结果与Lehmkuhl等^[10]利用Mode H的DNS计算结果之间的相对误差较小, 分别为6.8%和3.7%, 而对Case AV和CaseCV工况, 本文计算结果与Lehmkuhl等^[10]利用Mode H之DNS计算结果的相对误差较大, 分别为29.5%和29.3%.

在位于回流区的站位${{{x_1}} /D}$ = 1.54处, 对Case AV—Case DV工况, 本文计算结果与Lehmkuhl等^[10]利用Mode H的DNS计算结果的相对误差分别为7.75%, 15.2%, 14%和1%. 在位于回流区的站位${{{x_1}}/D}$ = 2.02处, 对Case AV—Case DV工况, 本文计算结果与Lehmkuhl等^[10]利用Mode H之DNS计算结果的相对误差分别为1.7%, 1.2%, 0.5%, 2.2%.

图10分别给出了上述8种工况下在3个不同站位(${{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} D}} \right. } D}$ = 1.06, 1.54, 2.02)处各向异性雷诺应力剖面特性的数值结果. 对此情况, Lehmkuhl等^[10]没有给出相关的DNS之计算结果.

由图10可见, 对图10(a)的情况, 在位于剪切层小尺度K-H不稳定性发生区域的站位${{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} D}} \right. } D}$ = 1.06处, 本文在Case AU—Case DU这4种工况下的计算结果均与Parnaudeau等^[18]的实验结果具有良好的一致性. 在位于回流区的站位${{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} D}} \right. } D}$ = 1.54处, 本文对Case AU和Case BU的计算结果与Parnaudeau等^[18]的实验结果也具有良好的一致性, 但在$\left| {{{{x_2}} /D}} \right| < 0.48$的范围内本文对Case CU和Case DU的计算结果, 与Parnaudeau等^[18]的实验结果有一定的差异, 最大相对误差分别为38.7%和14.8%. 在位于回流区的站位${{{x_1}} /D}$ = 2.02处, 本文对Case BU—Case DU的计算结果与Parnaudeau等^[18]的实验结果具有良好的一致性, 但在$\left| {{{{x_2}} / D}} \right| < 0.48$的范围内本文对Case AU的计算结果, 与Parnaudeau等^[18]的实验结果有一定的差异, 最大相对误差为30.8%.

对图10(b)的情况, 在位于剪切层小尺度K-H不稳定性发生区域的站位${{{x_1}}/D}$ = 1.06处, 以及在位于回流区的站位${{{x_1}} /D}$ = 1.54处, 本文对Case AV—Case DV这4种工况下的计算结果均与Lourenco和Shih^[27]的实验结果具有良好的一致性. 在位于回流区的站位${{{x_1}} /D}= 2.02$ 处, 本文对Case AV—Case DV的计算结果与Lourenco和Shih^[27]的实验结果相比, 在“峰”和“谷”处有一定的差异. 其中, 在“峰”处的相对误差分别为30.2%, 29.5%, 38%, 43.8%, 而在“谷”处的相对误差分别为32.7%, 31.8%, 41.1%, 40.8%.

Tremblay^[8]采用DNS和LES同时也均复现了Lourenco和Shih^[27]对雷诺数Re = 3900下, 圆柱绕流后方不同站位处雷诺应力剖面特性的实验结果. 其中, LES亚格子模型分别包括经典Smagorinsky模型(Smago)及动态模型(Dyn). 结果表明, DNS的结果要优于LES的结果, 而Dyn亚格子模型的结果要优于Smago亚格子模型的结果.

对PANS模型, 其对圆柱后方各站位处雷诺应力剖面的计算精度同样地依赖于滤波控制参数${f_{\text{k}}}$的取值方式. 对DES类模型, 其对圆柱后方各站位处雷诺应力剖面的计算精度则依赖于所选用RANS模型的类型及RANS/LES之间的混合转换方法等.

由于这两类模型均采用了相关的RANS模型, 而RANS模型对湍流尺度的解析能力是基于对涡黏系数的调控而实现, 但这种调控与计算网格的空间分布特性对湍流尺度的系统性分辨率直接相关. 通常情况下, 虽然可以通过加密网格来调控这两类模型对湍流尺度的解析能力, 但在近壁面网格较密的条件下, 会出现对网格空间分布特性过于敏感的缺陷, 导致近壁面RANS区域被破坏, 而网格空间分布尺度还没有小到足够进行大涡模拟的程度^[25].

本文采用WM-HRL模型的系列数值模拟结果表明, 在同一套网格系统下, 通过改变WM-HRL模型中“两个边界位置”及“三个区域的湍动能解析度”的组合, 既可通过数值模拟获得与Parnaudeau等^[18]实验结果一致的圆柱后方各站位处的雷诺应力剖面, 也可以通过数值模拟获得与Lourenco和Shih^[27]实验结果一致的圆柱后方各站位处的雷诺应力剖面.

研究发现, 对本文所构建的WM-HRL模型, 在RANS/LES混合转换区边界位置及其两个指标参数$n{r_{{\text{k1-Q}}}}$和$n{r_{{\text{k2-Q}}}}$准则取值的某些组合下, 计算所得二阶统计量与相关文献的实验结果相比, 存在较大的误差, 其主要原因如下.

第一, 在本文所构建的WM-HRL模型中, 其RANS模型采用的是SST $k{\text{-}}\omega$模型. 该湍流模型在边界层内为标准的Wilcox $k{\text{-}}\omega$两方程模型(简称SKW模型). 对WM-HRL模型, 其RANS模型必须直接作用于圆柱近壁面处, 但对SKW模型, 会存在如下缺陷: 包括可能会过大地预测回流区的长度, 还可能会过大地预测壁面剪切应力的值及过小地预测壁面湍动能的值^[44].

第二, 在本文所构建的WM-HRL模型中, LES模型采用的是经典的Smargorinsky型亚格子模型. 该亚格子模型虽然概念简单且数值求解稳定性高, 但不仅会在近壁面产生过大的数值耗散, 而且定常的Smagorinsky常数难以描述复杂时空变化下的湍流场.

为克服上述两个缺陷, 对RANS模型可改用Peng等^[44]提出的低雷诺数Wilcox $k{\text{-}}\omega$两方程模型的修正版本. Peng等^[44]的研究表明, 该修正版本可有效地提高其对近壁面流的解析能力, 包括可减小近壁面处$\omega$的值及增大近壁面处k的值等, 并具有更好地模拟分离、回流及再附等复杂流动现象的能力.

另一方面, 对LES模型可改用动态(Dyn)亚格子模型. Tremblay^[8]的研究表明, 对雷诺数Re = 3900下的圆柱绕流场问题, 与基于经典Smagorinsky型亚格子涡黏系数的LES模型相比, 采用基于动态(Dny)亚格子涡黏系数的LES模型, 可以更好地获得圆柱绕流场相关统计量的计算结果. 这些正是本文作者团队将在下一阶段重点开展的相关研究工作.

4.3   相干结构频谱及流场特征

首先研究剪切层小尺度K-H不稳定性结构发生区域的范围问题. 为此, 选取Case CU和Case CV两种工况进行分析. 研究对象为这两种工况下, 在图3中P1—P12这12个监测点处横向脉动速度的Lomb谱特征, 结果如图11和图12所示.

由图11和图12可见, 在12个监测点P1—P12处横向脉动速度的Lomb谱中, 均出现一个频率相对较低的谱峰, 该峰值所对应频率与涡泄频率一致, 其无因次频率${\bar f_{{\text{vs}}}} \approx 0.22$. 同时, 在12个监测点P1—P12处横向脉动速度的Lomb谱中, 还会出现一个与该频率峰值相对应的倍频($2{\bar f_{{\text{vs}}}}$)成分. 前者为圆柱体升力脉动的主频成分, 而后者为圆柱体阻力脉动的主频成分.

对Case CU工况, 在测点P4—P10处横向脉动速度的Lomb谱中(见图11), 还可以观察到一个频率更高的谱峰, 该峰值所对应频率与剪切层小尺度K-H不稳定性的频率一致, 其无因次频率${\bar f_{{\text{kh}}}}$在1.33—1.47之间. 进一步的研究表明, 在此工况下圆柱绕流场中剪切层小尺度K-H不稳定性结构发生在对数律层中$89.4 \leqslant {y^ + } \leqslant 253.5$的区域.

对Case CV工况, 在测点P4—P7处横向脉动速度的Lomb谱中(见图12), 也可以观察到一个频率更高的谱峰, 该峰值所对应频率也与剪切层小尺度K-H不稳定性的频率一致, 其无因次频率${\bar f_{{\text{kh}}}}$也在1.43—1.44之间. 进一步的研究表明, 在此工况下圆柱绕流场中剪切层小尺度K-H不稳定性结构发生在对数律层中$89.4 \leqslant {y^ + } \leqslant 167.4$的区域.

讨论其他6个监测点P13—P18处流向和横向脉动速度的Lomb谱特征, 计算工况选取Case CU和Case CV两种情况, 结果如图13所示. 其中, 监测点P13和P14分别位于剪切层偏上及偏下的位置, 监测点P15位于剪切层脱落区, 监测点P16位于回流区内的尾流中线上, 而监测点P17和P18位于回流区外的尾流中线上.

由图13可见, 在Case CU和Case CV两种工况下, 对P13—P18这6个监测点, 无论是在其流向脉动速度的Lomb谱中, 还是在其横向脉动速度的Lomb谱中, 均清晰地出现两个频率较低的谱峰, 一个谱峰对应涡泄频率(${\bar f_{{\text{vs}}}}$), 另一个谱峰对应涡泄频率的倍频($2{\bar f_{{\text{vs}}}}$).

对位于剪切层脱落区的监测点P15的情况, 在Case CU和Case CV两种工况下, 无论是在其流向脉动速度的Lomb谱中, 还是在其横向脉动速度的Lomb谱中, 除了均清晰地出现单频频率${\bar f_{{\text{vs}}}}$及其倍频$2{\bar f_{{\text{vs}}}}$的谱峰外, 还可观察到一个三倍频率$3{\bar f_{{\text{vs}}}}$的谱峰.

对位于其他区域的监测点P16—P18的情况, 在Case CU和Case CV两种工况下, 从其流向脉动速度的Lomb谱中只能清晰地看到单频频率${\bar f_{{\text{vs}}}}$及其倍频$2{\bar f_{{\text{vs}}}}$的谱峰, 而从其横向瞬态速度的Lomb频谱图中则只能清晰地看到单频频率${\bar f_{{\text{vs}}}}$及其三倍频率$3{\bar f_{{\text{vs}}}}$的谱峰, 但没有出现明显的二倍频率$2{\bar f_{{\text{vs}}}}$的谱峰.

由图13进一步可见, 在Case CU和Case CV两种工况下, 对位于剪切层偏上方的监测点P13, 无论是在其流向脉动速度的Lomb谱中, 还是在其横向脉动速度的Lomb谱中, 除了均清晰地出现与涡泄相关的单频频率${\bar f_{{\text{vs}}}}$及倍频频率$2{\bar f_{{\text{vs}}}}$的谱峰外, 还可清晰地观察到一个频率更高的谱峰, 其对应频率与剪切层小尺度K-H不稳定性结构的频率一致, 其无因次频率${\bar f_{{\text{kh}}}}$在1.43—1.44之间.

在Case CU和Case CV两种工况下, 对位于剪切层偏下方的监测点P14, 在其流向脉动速度的Lomb谱图中, 只清晰地出现了与涡泄相关的单频频率${\bar f_{{\text{vs}}}}$及倍频频率$2{\bar f_{{\text{vs}}}}$的谱峰. 然而, 在其横向脉动速度的Lomb谱图中, 除了均清晰地出现与涡泄相关的单频频率${\bar f_{{\text{vs}}}}$及倍频频率$2{\bar f_{{\text{vs}}}}$的谱峰外, 还可清晰地观察到一个与剪切层小尺度K-H不稳定性结构频率一致的谱峰, 其无因次频率${\bar f_{{\text{kh}}}}$也在1.43—1.44之间.

最后, 从图13可见, 在Case CU和Case CV两种工况下, 对监测点P13和P14, 无论是在其流向脉动速度的Lomb频谱图中, 还是在其横向脉动速度的Lomb谱中, 均没有出现明显的三倍频率$3{\bar f_{{\text{vs}}}}$的谱峰.

对雷诺数Re = 3900下圆柱绕流场中出现三倍频率的现象, 在Parnaudeau等^[18]的DNS及Pereira^[2]等的PANS计算中都有发现, 而在Lehmkuhl等^[10]采用DNS及Tremblay^[8]和Breuer^[15]采用LES的计算中则没有发现. Ong和Wallace^[45]指出, 流体行为中的峰值现象可以通过线性稳定性理论进行预测, 用于研究各种不稳定性问题. 然而, 这个理论是否能够适用于本文所观察到的这种复杂现象, 或者本文算例所出现的这类复杂频谱现象是否还存在其他影响因素, 值得在后续工作中做进一步的研究.

最后讨论圆柱绕流场中两类相干结构的流场特征, 计算工况选取Case AU—DU及Case AV—DV两类情况, 结果如图14和图15所示. 其中, 图14的左列为${{{x_3}} \mathord{\left/ {\vphantom {{{x_3}} D}} \right. } D} = {{\text{π }} \mathord{\left/ {\vphantom {{\text{π }} 2}} \right. } 2}$断面上展向涡量${\omega _3}= {{\partial {u_2}} \mathord{\left/ {\vphantom {{\partial {u_2}} {\partial {x_1}}}} \right. } {\partial {x_1}}} - {{\partial {u_1}} \mathord{\left/ {\vphantom {{\partial {u_1}} {\partial {x_2}}}} \right. } {\partial {x_2}}}$的云图, 右列为流向速度云图, 且图中白色垂直虚线为回流区末端位置(即${{{u_1}} \mathord{\left/ {\vphantom {{{u_1}} {{U_\infty }}}} \right. } {{U_\infty }}} = 0$位置). 图15为沿圆柱体长度的三维展向涡量云图.

在Case AU—DU及Case AV—DV两类工况下, 从图14 (左)均清晰地可见附着于圆柱上下侧壁的两个层流分离剪切层结构, 表现为两条位置狭小的窄带结构. 同时, 从图14 (右)则可清晰地观察到圆柱尾部因流动分离而形成的两个回流区分支结构. 具体来说, 第1个分支结构发生在${{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} D}} \right. } D}$ = 1.06处平均流向速度剖面为U形的情况, 回流区长度较长; 而第2个分支结构发生在${{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} D}} \right. } D}$ = 1.06处平均流向速度剖面为V形的情况, 回流区长度较短.

结合图11和图12的结果进一步可知, 在Case AU—DU及Case AV—DV两类工况下, 两个剪切层中小尺度K-H不稳定性结构的触发位置都基本相同, 均约位于${y^ + } = 89.4$处, 但两个剪切层中小尺度K-H不稳定性结构的消逝位置不同. 在Case AU—DU工况下, 其消逝位置约为${y^ + } = 253.5$. 而在Case AV—DV工况下, 其消逝位置约为${y^ + } = 167.4$.

这意味着, 当${{{x_1}}/D}$ = 1.06处平均流向速度剖面为U形时, 剪切层中小尺度K-H不稳定性结构的长度较大, 而且相应的回流区长度也较长. 而当${{{x_1}}/D}$ = 1.06处平均流向速度剖面为V形时, 剪切层中小尺度K-H不稳定性结构的长度较小, 而且相应的回流区长度也较短. 因此, 对雷诺数Re = 3900下圆柱绕流中出现两类回流区分支结构, 以及在圆柱尾部近壁面处平均流向速度出现U形和V形两类剖面结构形状的形成机理, 可能与WM-HRL模型中两个边界位置及3个区域的湍动能解析度发生变化时, 该湍流模型对圆柱绕流剪切层中小尺度层流分离流动的解析能力不同有关.

对DNS模型, 在同一套网格系统下, 其对圆柱绕流剪切层中小尺度层流分离流动的解析能力是确定的, 因此只能获得某一种回流区分支结构及与之对应的圆柱尾部近壁面平均流向速度剖面的某一种形状. 但对不同的网格系统, 由于其对圆柱绕流剪切层中小尺度层流分离流动的解析能力相应地会不同, 进而可能会出现在不同网格系统下获得不同的回流区分支结构及与之对应的圆柱尾部近壁面流向速度剖面的U形或V形结构.

对PANS及DES类模型, 它们可选用不同的RANS湍流模型, 也可以改变湍流模型中的相关模型参数. 对DES类模型, 其执行RANS/LES之间转换过渡的方式还可以不同. 另一方面, 对LES模型, 其可以选用经典Smagorinsky及Dyn等^[9]不同的亚格子涡黏模型. 这些复杂的因素都可能会导致即使在同一套网格系统下, 使得这些湍流模型对圆柱绕流剪切层中小尺度层流分离流动的解析能力不同, 进而可能会获得不同的回流区分支结构及与之对应的圆柱尾部近壁面平均流向速度剖面的U形或V形结构.

亚临界雷诺数区圆柱绕流中主要有两类相干结构. 其中, 一类为剪切层小尺度K-H不稳定性结构, 而另一类为大尺度的Karman涡街结构. 为进一步阐述它们的形成机理及其特征, 在图15中, 给出了在Case AU—DU和Case AV—DV两类工况下, 圆柱绕流展向三维涡量云图的数值结果.

由图15并结合图14可见, 由于剪切层两侧的速度差诱发K-H不稳定性, 使剪切层失稳并出现滚卷现象. 一方面, 这种现象会在剪切层中导致大量小尺度的旋涡产生于此区域, 在横向瞬态速度的Lomb谱中表现为宽频信号的凸起(见图11和图12). 另一方面, 这种现象还会在圆柱尾迹区形成另一类大尺度的旋涡结构, 即Karman涡街结构, 在横向瞬态速度的Lomb谱中表现为频率相对较低的窄频信号的凸起(见图11—图13).

由图15还可以清晰地观察到, 附着于圆柱壁面的剪切层伴随着其失稳、滚卷等复杂现象在圆柱尾流区形成流向与展向相嵌套的三维涡体结构. 由于本文所构建的WM-HRL模型具有能够精细解析湍流谱结构的能力(见图11—图13), 进而不仅能够精细地解析各种小尺度涡结构(如K-H不稳定性结构等), 又能够精细地解析各种大尺度涡结构(如Karman涡街等). 因此, 该WM-HRL模型表现出预期的相当于LES甚至DNS对各种小尺度及大尺度运动解析的能力.

5.   结　论

本文通过修改SST-$k{\text{-}}\omega$湍流模型中k方程的色散项, 构造了一种新的具有壁面模化能力的钝体绕流混合RANS/LES模型, 即WM-HRL模型. 该模型在钝体近壁面的某个计算区域内采用RANS模式, 在经历一个RANS/LES混合转换过渡区后, 在其余计算区域则采用LES模式, 且该LES模式等效为经典Smagorinsky亚格子模型.

通过构造一个仅与当地网格空间分布尺寸相关的湍动能解析度指标参数${r_{\text{k}}}$, 提出了确定该WM-HRL模型两个边界位置和3个区域湍动能解析度的 $n{r_{{\text{k1-Q}}}}$和$n{r_{{\text{k2-Q}}}}$准则. 其中, $n{r_{{\text{k1-Q}}}}$准则用于确定WM-HRL模型中RANS模式结束边界位置, 及其在RANS区湍动能的最大解析度(即${r_{{\text{k}}1}}$), 而$n{r_{{\text{k2-Q}}}}$准则用于确定WM-HRL模型中LES模式启动边界位置及其在LES区对湍动能的最小解析度(即${r_{{\text{k2}}}}$).

在此基础上, 本文构造了一个新的具有双重保护功能的混合函数${f_{{\text{nd}}}}$. 第一重保护通过$n{r_{{\text{k1-Q}}}}$准则实现, 保证在${r_{\text{k}}} \geqslant {r_{{\text{k}}1}}$时WM-HRL为完全RANS模式. 通过第一重保护层的设置, 可使WM-HRL既具有延迟脱体涡模拟(DDES)的能力, 又具有壁面模化(WM)的能力. 在RANS区最大湍动能解析度指标${r_{{\text{k}}1}}$的合理设置下, 可避免MSD和GIS的问题.

第二重保护通过$n{r_{{\text{k2-Q}}}}$准则实现, 可保证在${r_{\text{k}}} \leqslant {r_{{\text{k}}2}}$时WM-HRL为完全LES模式. 在此第二重保护层设置的条件下, 可通过${r_{{\text{k}}2}}$值的合理设置, 一方面在进行计算网格设置时即可实现对LES启动边界位置的预先设定, 另一方面还可自定义LES区中WM-HRL对钝体绕流各种复杂湍流场的解析能力.

在合理设置$n{r_{{\text{k1-Q}}}}$和$n{r_{{\text{k2-Q}}}}$准则的条件下, 既可保证WM-LES模型在RANS/LES混合转换区即可激活LES模式, 同时还可保证在完全LES区即可完全激活LES模式, 进而避免LLM的问题.

基于该WM-HRL模型, 本文以雷诺数Re = 3900下圆柱绕流为对象, 开展了系列数值模拟分析与评估, 得到如下主要结论.

1)即使当LES的启动边界位于圆柱绕流剪切层小尺度K-H不稳定性区内时, 在合理设置的$n{r_{{\text{k2-Q}}}}$准则(一般需要${r_{{\text{k}}2}} \leqslant 0.5$)的条件下, 该WM-HRL模型也能获取准确的涡泄频率、分离角、阻力系数及基准线压力系数等统计量信息.

2)为使WM-HRL模型能够准确获取圆柱绕流剪切层小尺度K-H不稳定性结构特征及其频谱特性, 以及回流区两个分支结构的长度信息, 可将LES模式的启动边界位置均设置在过渡层内, 而RANS模式的结束边界位置则既可在黏性底层也可在过渡层内, 并且使${r_{{\text{k}}2}} \leqslant 0.2$, 即在对数律层区LES对湍动能具有至少具有80%的解析能力.

3)在上述2)的条件下, WM-HRL模型能够获取与相关实验精度相当的圆柱绕流一阶和二阶统计量的信息, 以及剪切层小尺度K-H不稳定性结构的无因次频率${\bar f_{{\text{kh}}}}$及回流区两个分支结构的长度信息, 并且当回流区长度的数值结果与Parnaudeau等^[18]的实验结果一致, 在${{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} D}} \right. } D} = 1.06$处的平均流向速度剖面为U形, 而当回流区长度的数值结果与Lourenco和Shih^[27]的实验结果一致, 则为V形.

特别地, 在本文的系列数值模拟中, 通过两个边界位置及3个区域的湍动能解析度的各种不同组合, 在仅用一套网格系统的条件下, 利用新构造的WM-HRL模型, 针对雷诺数Re = 3900下圆柱绕流场问题, 获得了两类长度不同的剪切层小尺度K-H不稳定性结构特性, 具体如下.

第一类剪切层小尺度K-H不稳定性结构的长度较长, 此时回流区长度也较长, 与之相应的圆柱尾部近壁面处的平均流向速度剖面为U形. 第二类剪切层小尺度K-H不稳定性结构的长度较短, 此时回流区长度也较短, 与之相应的圆柱尾部近壁面处的平均流向速度剖面为V形.

关于雷诺数Re = 3900下圆柱绕流场中回流区的两个分支结构, 以及与之相应的圆柱尾部近壁面处U形和V形两类平均流向速度剖面形状等问题, 已经成为困扰学界三十余年的难题. 本项研究为后续更深层次认识这一难题提供了一种新的手段.